#include <map>#include "TKey.h"#include <TH2F.h>#include <TProfile.h>#include <TROOT.h>#include <TRint.h>#include <TList.h>#include <TTree.h>#include <TFile.h>#include <TChain.h>#include "TChainElement.h"#include <TMath.h>#include <stdlib.h>#include <TSystem.h>#include <iostream>#include <sstream>#include "Log.hh"#include "MicroException.hh"#include <string>#include "mysql/Mysql.hh"#include "root/CaloRun.hh"#include "root/CaloEvent.hh"#include "root/CaloHit.hh"#include "root/MTRun.hh"#include "root/MTEvent.hh"#include "root/MTChannel.hh"#include "TGraphErrors.h"#include "TGraphAsymmErrors.h"#include <TF1.h>#include "tools/Arguments.hh"Include dependency graph for getNbOfHits.cpp:
Go to the source code of this file.
Functions | |
| int | main (int argc, char **argv) |
Variables | |
| const int | maxEvents = 200000 |
| const int | maxTrees = 100 |
| const int | maxHits = 5000 |
| const float | AverageMult = 1.7 |
| int main | ( | int | argc, | |
| char ** | argv | |||
| ) |
Definition at line 48 of file getNbOfHits.cpp.
00048 { 00049 Arguments arg(argc,argv," getNbOfHits run1 [run2] max [run20] -p particleType -r opt"); 00050 00051 if ( arg.getNbOfArgs() < 1) { 00052 arg.usage(); 00053 exit(1); 00054 } 00055 00056 // Analyse the command line to find what to do 00057 00058 int NbOfRuns = arg.getNbOfArgs(); 00059 00060 string rootName[20]; 00061 for (int nr=0 ; nr < NbOfRuns ; nr ++ ){rootName[nr].assign(arg.getArgument(nr+1)); } 00062 00063 string ParticleType = "Pions"; 00064 if ( arg.getOption("-p").size() != 0 ) 00065 { 00066 ParticleType = arg.getOption("-p"); 00067 } 00068 00069 string RootFileType = "EM"; // With Efficiency and Multiplicity Corrections 00070 if ( arg.getOption("-r").size() != 0 ) 00071 { 00072 RootFileType = arg.getOption("-r"); 00073 } 00074 00075 TString Period = "SPS_aug_2012"; 00076 if ( arg.getOption("-t").size() != 0 ) 00077 { 00078 Period = arg.getOption("-t"); 00079 } 00080 00081 cout << " Starting Analysis Run " << rootName[0] << endl; 00082 int search = 1; 00083 00084 int NbHitLayer[49]; 00085 int IsPenetratingMIP; 00086 00087 TH1I* hnhit_roi = new TH1I("hnhit_roi","Number of hits in tight DeltaT",200,0,2000.); 00088 TH1F* hDensity = new TH1F("hDensity","Hits Density ",200,0,50.); 00089 TH1I* hHitsTot48 = new TH1I("hHitsTot48", "Nb of Hits Thresh 1 Chamber 48", 200, 0, 200); 00090 TH1I* hHitsTot49 = new TH1I("hHitsTot49", "Nb of Hits Thresh 1 Chamber 49", 200, 0, 200); 00091 TH1I * hnHitForw = new TH1I("hnhitForw","MIP Number of hits Forw",50,0,50); 00092 TH2I* hxy = new TH2I("hxy","Hit Map on first layer",96, 0, 96, 96, 0, 96); 00093 TH2I* hxyOrg = new TH2I("hxyOrg","Hit Map on Central Pad",96, 0, 96, 96, 0, 96); 00094 TH2I* hPadCentr = new TH2I("hPadCentr","Central Pad",96, 0, 96, 96, 0, 96); 00095 TH2F* hXYExtr = new TH2F("hXYExtr","XY Extrapolation",96, 0, 96, 96, 0, 96); 00096 TH2F* hS3S9 = new TH2F("hS3S9","NB of Hits in 3x3/10x10 Nhits", 400,0.,2000., 100, 0., 1.); 00097 TProfile * hRegMult = new TProfile("hRegMult", "Nb of Hits in n+1xn+1 cumulative matrice", 50, 0, 50., 0.9, 1000.); 00098 TProfile * hRegMultC = new TProfile("hRegMultC", "Nb of Hits Total-center matri", 50, 0, 50., 0., 1000.); 00099 00100 TH2F* hRMSx = new TH2F("hRMSx","NB of Hits VS RMSX", 400,0.,2000.,100,0,25); 00101 TH2F* hRMSy = new TH2F("hRMSy","NB of Hits VS RMSY", 400,0.,2000.,100,0,25); 00102 TH2I* hNhitDeb = new TH2I("hNhitDeb", " Deb", 400,0.,2000.,400,0.,2000.); 00103 TH1F* hNhitDebrap = new TH1F("hNhitDebrap", " Debaa", 100, 0., 5.); 00104 TH2F* hNhitRapHits = new TH2F("hNhitRapHits", " Debaa", 400,0.,2000.,100, 0., 5.); 00105 TH1F* hRG1G2 = new TH1F("hRG1G2", "G1/G2", 100, 0., 10.); 00106 TH2F* hnHitsRG1G2 = new TH2F("hnHitsRG1G2", "Nhits vs G1/G3", 400,0.,2000.,100, 0., 10.); 00107 // if(rootName[0]=="715615" || rootName[0]=="715692" || rootName[0]=="715699" || rootName[0]=="715747"){xhlim = 700;} 00108 TH1I* hNbHitsPions = new TH1I("hhNbHitsPions","Number of hits for Pions",520,0,2600.); 00109 TH1I* hNbHitsPions1 = new TH1I("hhNbHitsPions1","Number of hits for Pions",520,0,2600.); 00110 TH1I* hNbHitsPions2 = new TH1I("hhNbHitsPions2","Number of hits for Pions",200,0,1000.); 00111 TH1I* hNbHitsPions3 = new TH1I("hhNbHitsPions3","Number of hits for Pions",100,0,500.); 00112 TH1I* hNbHitsPionsSC = new TH1I("hhNbHitsPionsSC","Number of hits for Pions",520,0,2600.); 00113 TProfile * hNoiseLay = new TProfile("hNoiseLay","Noise Per layer",50, 0., 50., 0., 100.,"S"); 00114 TH1F* hShCont1 = new TH1F("hShCont1","Rear to Forward Nb hits ratio", 100, 0., 1.); 00115 TH1F* hNbHitsRear = new TH1F("hNbHitsRear","Nb Hits Rear", 100, 0., 100.); 00116 TH1F* hNbHitsBoard = new TH1F("hNbHitsBoard","Nb Hits Board", 200, 0., 500.); 00117 00118 TH2I * hxyMap[50]; 00119 TH2I * hxyMap1[50]; 00120 TH2I * hxyMap2[50]; 00121 TH2I * hxyMap3[50]; 00122 00123 TString name,title; 00124 int NbinsEffx=12; 00125 int NbinsEffy=12; 00126 int hbin=96; 00127 for (int i=0;i<50;i++) 00128 { 00129 name = Form("hxyMap_%i",i); 00130 title = Form("Hit distribution in Layer %i",i); 00131 hxyMap[i]= new TH2I(name,title,96,0,hbin,96,0,hbin); 00132 hxyMap1[i]= new TH2I(name+" Th1",title,96,0,hbin,96,0,hbin); 00133 hxyMap2[i]= new TH2I(name+" Th2",title,96,0,hbin,96,0,hbin); 00134 hxyMap3[i]= new TH2I(name+" Th3",title,96,0,hbin,96,0,hbin); 00135 } 00136 TH2I * hxyTot = new TH2I("hxyTot","",96,0,hbin,96,0,hbin); 00137 00138 TString OutFilename = "/lapp_data/LC/Detecteurs/Sid/karyotak/DataTB2012/rootfiles/"+ParticleType+rootName[0]+RootFileType+".root"; 00139 TFile tf(OutFilename,"RECREATE"); 00140 00141 00142 cout << " Writting file : " << OutFilename << endl; 00143 00144 int nbEvtTot=0; 00145 int NbTrees = 0; 00146 int RunId=1; 00147 int nhwrt=0; 00148 int nhit_roi_old=1; 00149 int nhit_roi; 00150 float RG1G3 = -100.; 00151 int maxDeltatOld = -100; 00152 int maxDeltat = -100; 00153 UInt_t BcIdAbsOld =0; 00154 UInt_t BcIdAbs =0; 00155 string Ebeam; double BeamEnergy ; 00156 double nHitsTot; double nHitsTot1; double nHitsTot2; double nHitsTot3; 00157 double ShCont1; double NbHitsRear; double NbHitsBoard; int zstart; 00158 00159 TString PionsOutputFileName="/lapp_data/LC/Detecteurs/Sid/karyotak/DataTB2012/rootfiles/DSTPions"+rootName[0]+RootFileType+".root"; 00160 TFile* PionsRootFile = new TFile(PionsOutputFileName, "recreate"); 00161 TTree* Ptree = new TTree("Pions",""); 00162 Ptree->Branch("Ebeam", &BeamEnergy, "BeamEnergy/D"); 00163 Ptree->Branch("Nhtot", &nHitsTot, "nHitsTot/D"); 00164 Ptree->Branch("Nhtot1", &nHitsTot1, "nHitsTot/D"); 00165 Ptree->Branch("Nhtot2", &nHitsTot2, "nHitsTot/D"); 00166 Ptree->Branch("Nhtot3", &nHitsTot3, "nHitsTot/D"); 00167 Ptree->Branch("ShCont1", &ShCont1, "ShCont1/D"); 00168 Ptree->Branch("NbHitsRear", &NbHitsRear, "NbHitsRear/D"); 00169 Ptree->Branch("NbHitsBoard", &NbHitsBoard, "NbHitsBoard/D"); 00170 Ptree->Branch("Zstart", &zstart, "zstart/I"); 00171 00172 // Define Beam spot limits 00173 00174 Int_t YPads[10] = {52, 48, 63, 43, 63, 63, 44, 45, 45, 50}; 00175 Int_t XPads[10] = {58, 50, 57, 48, 58, 61, 45, 61, 55, 42}; 00176 Int_t YPadsL[10] = {10, 10, 10, 10, 25, 25, 10, 10, 10, 5}; 00177 Int_t YPadsH[10] = {10, 10, 10, 10, 10, 10, 10, 25, 10, 5}; 00178 Int_t XPadsL[10] = {10, 10, 10, 10, 25, 25, 10, 25, 10, 5}; 00179 Int_t XPadsH[10] = {10, 10, 10, 10, 10, 10, 10, 10, 10, 5}; 00180 int jBeam = 9; 00181 TString Runs[9]={"714489","714525","714531","714546", "714553","714556","714559","714562","714573"}; 00182 for (int r=0; r<9 && ParticleType=="Pions" ; r++){ 00183 cout << r << " " << Runs[r] << " " << " Check run for Beam spot " << rootName[0] << endl; 00184 if( rootName[0]==Runs[r] ){jBeam = r; cout << "JBeam " << jBeam << " Beam spot at " << YPads[r] << " " << XPads[r] << endl; break; } 00185 } 00186 cout << " Found Beam spot at " << YPads[jBeam] << " " << XPads[jBeam] << endl; 00187 00188 TChain chain("T"); 00189 int nbHit = 0; 00190 Double_t Zlayer[maxHits]; // Before it was Float 00191 Double_t Xlayer[maxHits]; 00192 Double_t Ylayer[maxHits]; 00193 Double_t ErPos[maxHits]; 00194 00195 cout << " Reading Files : " << endl; 00196 for ( int nr=0; nr< NbOfRuns; nr++){ 00197 for (int i =0; i<5 ; i++){ 00198 TString SubFile = Form("_I%i",i); 00199 TString fname = Form("/lapp_data/LC/data/TB2012/"+Period+"/Production/CaloHits/CaloHit_DHCAL_"+rootName[nr]+SubFile+"_0.slcio.root"); 00200 chain.Add(fname); 00201 cout << fname << endl; 00202 } 00203 } 00204 Float_t EffMap[49]; Float_t MultMap[49]; 00205 for (int m=0; m<=49 ; m++){EffMap[m]=1.; MultMap[m] = AverageMult;} 00206 if(RootFileType =="EM" || RootFileType == "E" ){ // Load files only if needed 00207 // Get efficiencies and Multiplicities 00208 00209 TString EffFilename="/gpfs/LAPP-DATA/LC/Detecteurs/Sid/karyotak/DataTB2012/rootfiles/Muons"+rootName[0]+"Eff.root"; // Gets the efficiency of the 1st run 00210 // TString EffFilename="/lapp_data/LC/Detecteurs/Sid/karyotak/DataTB2012/rootfiles/MuonsEff"+rootName[0]+".root"; // Gets the efficiency of the 1st run 00211 TFile ChEff(EffFilename); 00212 TH1F* hGlEffChamber = (TH1F*) ChEff.Get("hGlEffChamber"); 00213 TProfile * hLayerMult = (TProfile*) ChEff.Get("hLayerMult"); 00214 cout << "Reading efficiencies and Multiplicities " << endl; 00215 for (int l = 0; l<=49; l++){ 00216 if(RootFileType=="E" || RootFileType=="EM"){EffMap[l] = 0.01*hGlEffChamber->GetBinContent(l+1);} 00217 if(RootFileType=="M" || RootFileType=="EM")MultMap[l]=hLayerMult->GetBinContent(l+1); 00218 cout << "Layer " << l << " Eff " << EffMap[l] << " Multiplicity " << MultMap[l] << endl; 00219 } 00220 } 00221 00222 TObjArray *fileElements=chain.GetListOfFiles(); 00223 TIter next(fileElements); 00224 TChainElement *chEl=0; 00225 int NchainFiles=0; 00226 while (( chEl=(TChainElement*)next() )) 00227 { 00228 TFile f(chEl->GetTitle()); 00229 // TFile f(rootName.c_str()); 00230 TIter nextkey(f.GetListOfKeys()); 00231 TKey *key; 00232 NchainFiles++; 00233 00234 int nbEvt = 0; 00235 00236 cout << " Reading a new file in the chain " << NchainFiles << endl; 00237 while (key = (TKey*)nextkey()) 00238 { 00239 00240 TString keyname = key->GetName(); 00241 if ((keyname == "hResidualX") || (keyname == "hResidualY")){} 00242 else 00243 { 00244 00245 TTree *tree = (TTree*)key->ReadObj(); 00246 cout << "---------------------------- NEW TTREE -----------------"<< endl; 00247 nbHit = 0; 00248 NbTrees++; 00249 if( NbTrees > maxTrees ) break; 00250 CaloRun *run=NULL; 00251 00252 try { 00253 CaloRun *run=dynamic_cast<CaloRun *>(tree->GetUserInfo()->First()); 00254 if( run != NULL ){cout << "Run[" << run->GetRunId() << "], temperature[" << run->GetTemperature() << "] pression[" << run->GetPressure() << "] " << endl ; } 00255 00256 00257 // Get Beam energy form runId 00258 Mysql mysql; 00259 try { 00260 Ebeam = mysql.getRunEnergy(run->GetRunId()); 00261 } 00262 catch (MicroException e) {} 00263 int le = Ebeam.find(" "); string EbeamN = Ebeam.substr(0,le); BeamEnergy = strtod (Ebeam.c_str(),NULL); 00264 cout << "Beam[" << Ebeam << "] " << BeamEnergy << endl; 00265 } 00266 catch ( ... ) {} 00267 CaloEvent *evt = new CaloEvent(); 00268 TBranch *branch= tree->GetBranch("CaloEvent"); 00269 branch->SetAddress(&evt); 00270 00271 CaloHit* hit =NULL; 00272 int nbEntries = tree->GetEntries(); 00273 cout << " Nb of events to read " << nbEntries << endl; 00274 00275 for ( int evtNum = 0; evtNum < nbEntries ; evtNum++) 00276 { 00277 if(nbEvt>=maxEvents) break; 00278 nbEvt++; nbEvtTot++; 00279 00280 hxyTot->Reset(); 00281 for (int i=0; i<50; i++){hxyMap[i]->Reset();hxyMap1[i]->Reset();hxyMap2[i]->Reset();hxyMap3[i]->Reset();} 00282 tree->GetEntry(evtNum); 00283 00284 int nbHits = evt->GetNHits(); 00285 00286 //number of hits in region of interest 00287 nHitsTot=0.; nHitsTot1=0.; nHitsTot2=0.; nHitsTot3=0.; 00288 int nhit_fit = 0; 00289 float xpos = 0; 00290 float ypos = 0; 00291 float zpos = 0; 00292 int ncol = 0; 00293 int nrow = 0; 00294 int zSlot = 0; 00295 int deltat=0; 00296 int Thresh = 0; 00297 nhit_roi_old = nhit_roi; 00298 maxDeltatOld = maxDeltat; 00299 BcIdAbsOld = BcIdAbs; 00300 00301 nhit_roi = 0; 00302 float Rad1lay=10000.; 00303 int XOrg=1000; int YOrg=1000; 00304 maxDeltat = evt->GetDeltaTmax(); 00305 BcIdAbs = evt->GetBcIdAbs(); 00306 00307 for (int layer=0 ; layer <=49 ; layer++){ 00308 NbHitLayer[layer]=0; 00309 } 00310 for (int i=0 ; i<maxHits ;i++){ 00311 Zlayer[i]=0; 00312 Xlayer[i]=0; 00313 Ylayer[i]=0; 00314 double theta = (0.0136/100.)*sqrt((i+1)*2/1.757)*(1+0.038*log((i+1)*2/1.757)); 00315 ErPos[i] = sqrt(pow(10./sqrt(12.),2)+pow(2*28.*theta,2)); 00316 } 00317 00318 int NbHitLayer33=0; 00319 int NbHitLayer77=0; 00320 00321 for(int i=0 ; i < nbHits ; i++) 00322 { 00323 hit = (CaloHit*)evt->GetHits()->UncheckedAt(i); 00324 00325 zSlot = hit->GetLayer() ; 00326 deltat=hit->GetDeltaT(); 00327 xpos = hit->GetX()/1000. ; // unit = mm 00328 ypos = hit->GetY()/1000. ; 00329 zpos = hit->GetZ()/1000.; 00330 ncol = hit->GetColInChamber(); 00331 nrow = hit->GetRowInChamber(); 00332 Thresh = hit->GetThreshold(); 00333 // cout << Thresh << " " << endl; 00334 if(ParticleType=="Noise"){ 00335 if(zSlot <48 && abs(deltat-maxDeltat)<=2) continue; // This is for RPCs 00336 if(zSlot >=48 && abs(deltat-maxDeltat)<=5) continue; // For Micromegas 00337 } 00338 else{ 00339 if(zSlot <48 && abs(deltat-maxDeltat)>2) continue; // This is for RPCs 00340 if(zSlot >=48 && abs(deltat-maxDeltat)>5) continue; // For Micromegas 00341 } 00342 00343 // if(zSlot >=48 ) continue; // Ignore Micromegas 00344 00345 00346 nhit_roi++; // Nb of good hits 00347 hxyMap[zSlot]->Fill(ncol,nrow,1); 00348 if(Thresh == 1) {hxyMap1[zSlot]->Fill(ncol,nrow,1);} 00349 if(Thresh == 2) {hxyMap2[zSlot]->Fill(ncol,nrow,1);} 00350 if(Thresh == 3) {hxyMap3[zSlot]->Fill(ncol,nrow,1);} 00351 hxyTot->Fill(ncol,nrow); 00352 00353 if(zSlot==0){ 00354 // hxy->Fill(ncol,nrow); 00355 float rad = sqrt( (nrow-58.)*(nrow-58.)+(ncol-50.)*(ncol-50.) ); 00356 if(rad < Rad1lay) {XOrg=ncol; YOrg=nrow; Rad1lay=rad;} 00357 } 00358 // cout << XOrg << " Y org" << YOrg << endl; 00359 NbHitLayer[zSlot]++; 00360 if(abs(ncol-XOrg)<=2 && abs(nrow-YOrg)<=2){NbHitLayer33++;} 00361 if(abs(ncol-XOrg)<=10 && abs(nrow-YOrg)<=10){NbHitLayer77++;} 00362 00363 00364 // Prepare for fitting a muon 00365 // Choose layers to make the fit 00366 if(zSlot>=48) continue; // Do not include MicroMegas in the fit 00367 if(nhit_fit==0){ 00368 Zlayer[nhit_fit] = zpos; // Convert coordinates to mm 00369 Xlayer[nhit_fit] = xpos; 00370 Ylayer[nhit_fit] = ypos; 00371 nhit_fit++; 00372 } 00373 else{ 00374 if (nhit_fit < maxHits){ 00375 if(NbHitLayer[zSlot] >1 ) continue ; // keep only 1st hit in layers 00376 // Select hits that have almost the same slope as the previous one 00377 Double_t slopy = TMath::ATan2( (ypos-Ylayer[nhit_fit-1]) , (zpos-Zlayer[nhit_fit-1]) ); 00378 Double_t slopx = TMath::ATan2( (xpos-Xlayer[nhit_fit-1]) , (zpos-Zlayer[nhit_fit-1]) ); 00379 if ( TMath::Abs(slopx)> 5. || TMath::Abs(slopy)>5. ) continue; 00380 Zlayer[nhit_fit] = zpos; 00381 Xlayer[nhit_fit] = xpos; 00382 Ylayer[nhit_fit] = ypos; 00383 nhit_fit++; 00384 } 00385 } 00386 } // end hit 00387 00388 // Ca c'est une rustine provisoire pour eviter les doublons d'evenements 00389 hNhitDeb->Fill(nhit_roi,nhit_roi_old); hNhitDebrap->Fill((float)nhit_roi/nhit_roi_old); 00390 hNhitRapHits->Fill((float)nhit_roi,(float)nhit_roi/nhit_roi_old); 00391 // if(fabs( (float)nhit_roi/nhit_roi_old - 1.)<=0.12 ) continue; 00392 00393 YOrg=hxyTot->ProjectionX()->GetMaximumBin()-1; 00394 XOrg=hxyTot->ProjectionY()->GetMaximumBin()-1; 00395 hxyOrg->Fill(YOrg,XOrg); 00396 float RmsY = hxyTot->GetRMS(1); 00397 float RmsX = hxyTot->GetRMS(2); 00398 hRMSx->Fill(nhit_roi,RmsX); 00399 hRMSy->Fill(nhit_roi,RmsY); 00400 00401 hnhit_roi->Fill(nhit_roi); 00402 00403 // Compute Hit Density to find MIP events 00404 int NbLayers=0; 00405 for (int layer=0 ; layer < 50 ; layer++){ 00406 if(NbHitLayer[layer]>0) 00407 { 00408 NbLayers++; 00409 } 00410 } 00411 00412 float HitDensity = 1. * nhit_roi/NbLayers; 00413 hDensity->Fill(HitDensity); 00414 00415 //Look for penetrating MIPS----------------------------- 00416 00417 int NbLayersForw=0; 00418 int NbLayersBack=0; 00419 00420 IsPenetratingMIP = 0; 00421 if(HitDensity>1 && HitDensity < 3){ 00422 00423 for (int layer=0 ; layer < 10 ; layer++){ 00424 if(0 < NbHitLayer[layer] && NbHitLayer[layer] <= 3){NbLayersForw++;} 00425 } 00426 for (int layer=38 ; layer < 48 ; layer++){ 00427 if(0 < NbHitLayer[layer] && NbHitLayer[layer] <= 3){NbLayersBack++;} 00428 } 00429 00430 hnHitForw->Fill(NbLayersForw); 00431 hnHitForw->Fill(NbLayersBack+20); 00432 00433 //Define penetrating MIP 00434 if(NbLayersForw>6 && NbLayersBack > 6) { IsPenetratingMIP = 1;} 00435 } // end of definition of Penetrating MIPS 00436 00437 hS3S9->Fill( (float)nhit_roi, (float)NbHitLayer33/NbHitLayer77); 00438 // Looking for Noise !------------------------------ 00439 if(ParticleType=="Noise"){ 00440 for (int nly=0 ; nly<=47 ; nly++){ 00441 nHitsTot = nHitsTot + hxyMap[nly]->GetEntries() ; 00442 hNoiseLay->Fill(nly, hxyMap[nly]->GetEntries()); 00443 hNbHitsPions->Fill( (int) nHitsTot ); } 00444 } 00445 // Looking for Pions --------------------------------------------------- 00446 // if( nhit_roi>150 && HitDensity>3 && ParticleType == "Pions"){ 00447 if( nhit_roi>80 && HitDensity>4 && ParticleType == "Pions" ){ 00448 hxy->Fill(YOrg,XOrg); 00449 00450 // if(nhwrt<15) {tf.cd();nhwrt++; hxyTot->Write();f.cd();} 00451 int lG1=4 ; int lG2=15 ; int lG30=30; 00452 int NbHitsG1 = 0; int NbHitsG2 = 0; int NbHitsG3 = 0; double NbHits30=0; NbHitsRear=0; NbHitsBoard=0; 00453 int zstartFound=0; zstart=0; 00454 for (int kl=0 ; kl<=47 ; kl++){ 00455 if(kl <=lG1 ) { NbHitsG1 = NbHitsG1 + (int)hxyMap[kl]->GetEntries();} 00456 if(kl > lG1 && kl< lG2 ) { NbHitsG2 = NbHitsG2+ (int)hxyMap[kl]->GetEntries();} 00457 if(kl > lG2 ) { NbHitsG3 = NbHitsG3 + (int)hxyMap[kl]->GetEntries();} 00458 if(kl <= lG30 ) { NbHits30 = NbHits30 + hxyMap[kl]->GetEntries();} 00459 if(kl >= 41 && kl< 47 ) { NbHitsRear = NbHitsRear + hxyMap[kl]->GetEntries();} 00460 00461 /******** rising number of hits in 3 consecutive layer *********/ 00462 00463 if ( zstartFound==0 && ((int)hxyMap[kl+2]->GetEntries() > (int)hxyMap[kl+1]->GetEntries()) && 00464 ((int)hxyMap[kl+1]->GetEntries() > (int)hxyMap[kl]->GetEntries()) ){ 00465 for (int k=0 ; k<3 ; k++){ 00466 /******* shower begins when nhit >= 4 ***********/ 00467 if ((int)hxyMap[kl+k]->GetEntries() >=4){ 00468 zstartFound=1; 00469 zstart = kl+k; 00470 break; 00471 } 00472 } 00473 } 00474 } 00475 RG1G3 = (float)NbHitsG1 / NbHitsG3; 00476 hRG1G2->Fill( (float)NbHitsG1/NbHitsG3 ); 00477 hnHitsRG1G2->Fill(nhit_roi, (float)NbHitsG1/NbHitsG3 ); 00478 ShCont1 = NbHitsRear / NbHits30; 00479 00480 00481 //--------------Define Pions--------------- 00482 if(RG1G3 < 0.5 ) { 00483 if( // Fiducial cuts 00484 (XPads[jBeam]-XPadsL[jBeam] < XOrg && XOrg < XPads[jBeam]+XPadsH[jBeam]) && 00485 (YPads[jBeam]-YPadsL[jBeam] < YOrg && YOrg < YPads[jBeam]+YPadsH[jBeam]) ) { 00486 // cout << " Hidsto" <<ShCont1 << " " << NbHitsRear << endl; 00487 hShCont1->Fill(ShCont1); hNbHitsRear->Fill(NbHitsRear); 00488 int nHits48 = (int)hxyMap[48]->GetEntries(); 00489 int nHits49 = (int)hxyMap[49]->GetEntries(); 00490 hHitsTot48->Fill(nHits48); 00491 hHitsTot49->Fill(nHits49); 00492 for (int nly=0 ; nly<=47 ; nly++){ 00493 if( EffMap[nly] > 0) {nHitsTot = nHitsTot + (AverageMult)/MultMap[nly] * hxyMap[nly]->GetEntries() / EffMap[nly] ;} 00494 if( EffMap[nly] > 0) {nHitsTot1 = nHitsTot1 + (AverageMult)/MultMap[nly] * hxyMap1[nly]->GetEntries() / EffMap[nly] ;} 00495 if( EffMap[nly] > 0) {nHitsTot2 = nHitsTot2 + (AverageMult)/MultMap[nly] * hxyMap2[nly]->GetEntries() / EffMap[nly] ;} 00496 if( EffMap[nly] > 0) {nHitsTot3 = nHitsTot3 + (AverageMult)/MultMap[nly] * hxyMap3[nly]->GetEntries() / EffMap[nly] ;} 00497 NbHitsBoard = NbHitsBoard + (hxyMap[nly]->GetEntries() - hxyMap[nly]->Integral(9,88,9,88)); 00498 } 00499 // if(nHitsTot<180 && nhwrt<50) {tf.cd();nhwrt++; hxyTot->Write();f.cd();} 00500 hNbHitsBoard->Fill(NbHitsBoard); 00501 hNbHitsPions->Fill( (int) nHitsTot ); 00502 hNbHitsPions1->Fill( (int) nHitsTot1 ); 00503 hNbHitsPions2->Fill( (int) nHitsTot2 ); 00504 hNbHitsPions3->Fill( (int) nHitsTot3 ); 00505 if(ShCont1 <0.1 || NbHitsRear < 6) {hNbHitsPionsSC->Fill( (int) nHitsTot );} 00506 // cout << "Tree " << ShCont1 << " " << NbHitsRear << endl; 00507 Ptree->Fill(); 00508 }//end of fiducial cuts 00509 } // end of selection on g1/g3 00510 } // End of Pions 00511 00512 // Now I select penetrating Muons 00513 // cout << IsPenetratingMIP << " " << nhit_roi << endl; 00514 if(IsPenetratingMIP==1 && nhit_roi<150 && nhit_roi>50 && ParticleType == "Muons" ){ 00515 float zmax = 1372.; 00516 TGraphErrors *grx = new TGraphErrors(nhit_fit,Zlayer,Xlayer,0,ErPos); 00517 TF1 *f1 = new TF1("f1","pol1",0,zmax); 00518 grx->Fit("f1","Q"); 00519 double betaX = f1->GetParameter(0); 00520 double alfaX = f1->GetParameter(1); 00521 double chi2X = f1->GetChisquare(); 00522 delete f1; delete grx; 00523 TGraphErrors *gry = new TGraphErrors(nhit_fit,Zlayer,Ylayer,0,ErPos); 00524 TF1 *f2 = new TF1("f2","pol1",0,zmax); 00525 gry->Fit("f2","Q"); 00526 double betaY = f2->GetParameter(0); 00527 double alfaY = f2->GetParameter(1); 00528 double chi2Y = f2->GetChisquare(); 00529 delete f2; delete gry; 00530 // if(nhwrt<32 & nhit_roi>8 ) { // Now check for doublons 00531 // tf.cd(); 00532 // nhwrt++; 00533 // cout << " num Evt local " << nbEvt << " new/old Nb of hits " << nhit_roi << " / " << nhit_roi_old << " new/old delatMax " << maxDeltatOld << " " << maxDeltat << " BcId " << BcIdAbs << " / " << BcIdAbsOld << endl ; 00534 // hxyTot->Write(); 00535 // f.cd(); 00536 // } 00537 for(int i=48; i<=49 && ( chi2X/(nhit_fit-2) < 6. && chi2Y/(nhit_fit-2) < 6 ); i++){ 00538 float Xextrapol = alfaX*i*28. + betaX ; // X extrapolate to chamber i in mm 00539 float Yextrapol = alfaY*i*28. + betaY ; // Y extrapolate to chamber i in mm 00540 hXYExtr->Fill(Yextrapol/10.,Xextrapol/10.); 00541 // Find Central pad from computed x,y position 00542 Int_t xpadCentr = -1000; 00543 Int_t ypadCentr = -1000; 00544 if( i>47 && ( fabs(Yextrapol-322.5) < 20 || fabs(Yextrapol-650.) < 20 ) ) continue; 00545 if( i>47 && ( fabs(Xextrapol-485.0) < 20 ) ) continue; 00546 if(i<48){ 00547 xpadCentr = (int)(Xextrapol/10.4); 00548 ypadCentr = (int)(Yextrapol/10.4); 00549 } 00550 else{ 00551 if(Xextrapol<=480.) {xpadCentr = (int)(Xextrapol/10.);} 00552 else {xpadCentr = (int)( (Xextrapol-5)/10 );} 00553 if(Yextrapol<=320.) {ypadCentr = (int)( (Yextrapol/10.) );} 00554 if(Yextrapol>325 && Yextrapol<=645.) {ypadCentr = (int)( (Yextrapol-5.)/10.);} 00555 if(Yextrapol>=646) {ypadCentr = (int)((Yextrapol-10.)/10.);} 00556 } 00557 hPadCentr->Fill(ypadCentr,xpadCentr); 00558 if (xpadCentr >8 && ypadCentr >8 && xpadCentr <88 && ypadCentr < 88 ){ 00559 // if (xpadCentr >0 && ypadCentr >0 ){ 00560 int nHits48 = (int)hxyMap[48]->GetEntries(); 00561 int nHits49 = (int)hxyMap[49]->GetEntries(); 00562 if(i==48){hHitsTot48->Fill(nHits48); } // cout << " Nb hits 48 " << nHits48 << endl; 00563 if(i==49){hHitsTot49->Fill(nHits49);} 00564 00565 00566 00567 } // end checking the pad center 00568 } // end loop on all layers 00569 } // end penetrating muon 00570 00571 if (nbEvtTot%10000 == 0 ) { 00572 cout << "Processing event: " << nbEvtTot 00573 << " Local event " << nbEvt << endl; 00574 } 00575 tf.cd(0); 00576 } // end event 00577 cout << " Finished with events " << nbEvt << " in this TREE " << endl; 00578 } // else 00579 } // end of tree or run 00580 cout << " Ending file in chain " << endl; 00581 } // end of chain 00582 cout << " Going to finsh the work after" << nbEvtTot << " events "<< endl; 00583 tf.cd(); 00584 hnhit_roi->Write(); 00585 hDensity->Write(); 00586 hHitsTot48->Write(); 00587 hHitsTot49->Write(); 00588 Double_t b148=hHitsTot48->GetBinContent(1); Double_t ent48=hHitsTot48->GetEntries(); cout << " Innefficiency 48 " << b148/ent48 << " 1st bin " << b148 << " All " << ent48 << endl; 00589 Double_t b149=hHitsTot49->GetBinContent(1); Double_t ent49=hHitsTot49->GetEntries(); cout << " Innefficiency 49 " << b149/ent49 << endl; 00590 hnHitForw->Write(); 00591 hxy->Write(); 00592 hRegMult->Write(); 00593 hRegMultC->Write(); 00594 hPadCentr->Write(); 00595 hXYExtr->Write(); 00596 hS3S9->Write(); 00597 hRMSx->Write(); 00598 hRMSy->Write(); 00599 hNhitDeb->Write(); 00600 hNhitDebrap->Write();hNhitRapHits->Write(); 00601 hRG1G2->Write() ; hnHitsRG1G2->Write(); 00602 hNbHitsPions->Write(); hNbHitsPions1->Write(); hNbHitsPions2->Write(); hNbHitsPions3->Write(); 00603 hNbHitsPionsSC->Write(); hNoiseLay->Write(); 00604 hxyOrg->Write(); 00605 hShCont1->Write();hNbHitsRear->Write();hNbHitsBoard->Write(); 00606 00607 PionsRootFile->cd(); 00608 Ptree->Write(); 00609 00610 }
| const int maxEvents = 200000 |
Definition at line 43 of file getNbOfHits.cpp.
| const int maxTrees = 100 |
Definition at line 44 of file getNbOfHits.cpp.
| const int maxHits = 5000 |
Definition at line 45 of file getNbOfHits.cpp.
| const float AverageMult = 1.7 |
1.4.7