#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 dependency graph for getShowerStart.cpp:
Go to the source code of this file.
Functions | |
| int | main (int argc, char **argv) |
Variables | |
| const int | maxEvents = 3000000 |
| const int | maxTrees = 100 |
| const int | maxHits = 200 |
| int main | ( | int | argc, | |
| char ** | argv | |||
| ) |
Definition at line 46 of file getShowerStart.cpp.
00046 { 00047 00048 if ( argc !=2 ) { 00049 FILE_LOG(logERROR) << "usage: example rootFile " << endl; 00050 exit(1); 00051 } 00052 00053 int search = 1; 00054 00055 00056 TH1F * hDensity = new TH1F("hDensity","Hits Density ",200,0,200.); 00057 TH2F * hxy = new TH2F("hxy","Hit Distribution",100,0,100.,100,0,100.); 00058 TH1I * hnhit = new TH1I("hnhit","Number of all hits",200,0,1000); 00059 TH1I * hShowerStart = new TH1I("hShowerStart","Start of the Shower",50, 0, 50); 00060 TH1I * hShowerStartDa = new TH1I("hShowerStartDa","Start of the Shower",50, 0, 50); 00061 TH1I * hz = new TH1I("hz","Nb of Hits per Layer",50,0,50); 00062 TProfile * tp_nhit = new TProfile("tp_nhit","Average Nb of Hits / Layer",50,0,50,0,200,""); 00063 TProfile * hLongprof = new TProfile("hLongprof","Long Profile reconstructed from 4 layers",50,0,50,0,150,""); 00064 TProfile * hLongprof1 = new TProfile("hLongprof1","Long Profile from 1 layer",50,0,50,0,150,""); 00065 TProfile * hLongprofL1 = new TProfile("hLongprofL1","Long Profile when Shower Starts Layer 1",50,0,50,0,150,""); 00066 TH1I * hnhit_roi = new TH1I("hnhit_roi","Number of hits in tight DeltaT",200,0,2000.); 00067 00068 TString name,title,name1,title1 ; 00069 int hbin=96; 00070 // Book Layer Hit maps 00071 TH2I * hxyMap[50]; 00072 TH2I * hxyMapTot[50]; 00073 for (int i=0;i<50;i++) 00074 { 00075 name = Form("hxyMap_%i",i); 00076 title = Form("Hit distribution in Layer %i",i); 00077 hxyMap[i]= new TH2I(name,title,96,0,hbin,96,0,hbin); 00078 name1 = Form("hxyMapTot_%i",i); 00079 title1 = Form("All events Hit distribution in Layer %i",i); 00080 hxyMapTot[i]= new TH2I(name1,title1,96,0,hbin,96,0,hbin); 00081 } 00082 TH2I * hxyTot = new TH2I("hxyTot","All Layers of 1 Event",96,0,hbin,96,0,hbin); 00083 00084 int NbHitLayer[100]; 00085 int NbHitLayer77[100]; 00086 float RG1G3 = -100.; 00087 00088 string rootName; 00089 rootName.assign(argv[1]); 00090 int nbHit = 0; 00091 int nbEvt = 0; 00092 TString RootFiles = "/lapp_data/LC/Detecteurs/Sid/karyotak/DataTB2012/rootfiles/" ; 00093 TString OutFilename = RootFiles+"ShowerStart"+rootName+".root"; 00094 TFile tf(OutFilename,"RECREATE"); 00095 int NbTrees = 0; 00096 int RunId=1; 00097 00098 TChain chain("T"); 00099 for (int i =0; i<5 ; i++){ 00100 TString SubFile = Form("_I%i",i); 00101 TString fname = Form("/lapp_data/LC/data/TB2012/SPS_H2_may_2012/Production/CaloHits/CaloHit_DHCAL_"+rootName+SubFile+"_0.slcio.root"); 00102 chain.Add(fname); 00103 cout << fname << endl; 00104 } 00105 00106 TObjArray *fileElements=chain.GetListOfFiles(); 00107 TIter next(fileElements); 00108 TChainElement *chEl=0; 00109 int NchainFiles=0; 00110 while (( chEl=(TChainElement*)next() )) 00111 { 00112 TFile f(chEl->GetTitle()); 00113 // TFile f(rootName.c_str()); 00114 TIter nextkey(f.GetListOfKeys()); 00115 TKey *key; 00116 NchainFiles++; 00117 cout << " Reading a new file in the chain " << NchainFiles << endl; 00118 while (key = (TKey*)nextkey()) 00119 { 00120 00121 TString keyname = key->GetName(); 00122 if ((keyname == "hResidualX") || (keyname == "hResidualY")){} 00123 else 00124 { 00125 00126 TTree *tree = (TTree*)key->ReadObj(); 00127 cout << "---------------------------- NEW TTREE -----------------"<< endl; 00128 nbHit = 0; 00129 NbTrees++; 00130 if( NbTrees > maxTrees ) continue; 00131 CaloRun *run=NULL; 00132 00133 try { 00134 CaloRun *run=dynamic_cast<CaloRun *>(tree->GetUserInfo()->First()); 00135 if( run != NULL ){cout << "Run[" << run->GetRunId() << "], temperature[" << run->GetTemperature() << "] pression[" << run->GetPressure() << "] " << endl ; } 00136 00137 } 00138 catch ( ... ) {} 00139 CaloEvent *evt = new CaloEvent(); 00140 TBranch *branch= tree->GetBranch("CaloEvent"); 00141 branch->SetAddress(&evt); 00142 00143 CaloHit* hit =NULL; 00144 int nbEntries = tree->GetEntries(); 00145 cout << " Nb of events to read " << nbEntries << endl; 00146 TString name = Form("deltaTRun_%i",RunId); 00147 TString title = Form("delta distribution Run %i",RunId); 00148 TH1I * deltaTRun = new TH1I(name,title,20,-10,10); 00149 00150 for ( int evtNum = 0; evtNum < nbEntries ; evtNum++) 00151 { 00152 tree->GetEntry(evtNum); 00153 nbEvt++; 00154 for (int i=0; i<50; i++){hxyMap[i]->Reset();} 00155 if(nbEvt>maxEvents) break; 00156 00157 int nbHits = evt->GetNHits(); 00158 hnhit->Fill(nbHits); 00159 00160 //number of hits in region of interest 00161 float xpos = 0; 00162 float ypos = 0; 00163 float zpos = 0; 00164 int ncol = 0; 00165 int nrow = 0; 00166 int zSlot = 0; 00167 int deltat=0; 00168 int nhit_roi = 0; 00169 float Rad1lay=10000.; 00170 int maxDeltat = evt->GetDeltaTmax(); 00171 int XOrg=1000; int YOrg=1000; 00172 for (int layer=0 ; layer <100 ; layer++){NbHitLayer[layer]=0;NbHitLayer77[layer]=0;} 00173 00174 for(int i=0 ; i < nbHits ; i++) 00175 { 00176 hit = (CaloHit*)evt->GetHits()->UncheckedAt(i); 00177 deltat=hit->GetDeltaT(); 00178 deltaTRun->Fill(maxDeltat-deltat); 00179 00180 xpos = hit->GetX()/10000. ; // unit = mm 00181 ypos = hit->GetY()/10000. ; 00182 zSlot = hit->GetLayer() ; 00183 zpos = hit->GetZ()/10000.; 00184 ncol = hit->GetColInChamber(); 00185 nrow = hit->GetRowInChamber(); 00186 00187 if(zSlot <48 && abs(deltat-maxDeltat)>2) continue; // This is for RPCs 00188 if(zSlot >=48 && abs(deltat-maxDeltat)>5) continue; // For Micromegas 00189 00190 nhit_roi++; // Nb of good hits 00191 00192 hz->Fill(zSlot); 00193 00194 hxy->Fill(ypos,xpos); 00195 hxyMap[zSlot]->Fill(ncol,nrow,1); 00196 hxyMapTot[zSlot]->Fill(ncol,nrow,1); 00197 hxyTot->Fill(ncol,nrow); 00198 00199 } // end hit 00200 hnhit_roi->Fill(nhit_roi); 00201 00202 if( nhit_roi<150 ) continue; // Reject muons already 00203 00204 int lG1=4 ; int lG2=15 ; 00205 int NbHitsG1 = 0; int NbHitsG2 = 0; int NbHitsG3 = 0; 00206 YOrg=hxyTot->ProjectionX()->GetMaximumBin()-1; // We get the central pad where beam goes for this event 00207 XOrg=hxyTot->ProjectionY()->GetMaximumBin()-1; 00208 // cout << " YOrg " << YOrg << " XOrg " << XOrg << endl; 00209 00210 for (int kl=0 ; kl<=49 ; kl++){ 00211 NbHitLayer[kl] = NbHitLayer[kl] + (int)hxyMap[kl]->GetEntries(); 00212 NbHitLayer77[kl] = hxyMap[kl]->Integral(YOrg-10,XOrg-10,YOrg+10, XOrg+10); 00213 if(kl <=lG1 ) { NbHitsG1 = NbHitsG1 + (int)hxyMap[kl]->GetEntries();} 00214 if(kl > lG1 && kl< lG2 ) { NbHitsG2 = NbHitsG2+ (int)hxyMap[kl]->GetEntries();} 00215 if(kl > lG2 ) { NbHitsG3 = NbHitsG3 + (int)hxyMap[kl]->GetEntries();} 00216 } 00217 RG1G3 = (float)NbHitsG1 / NbHitsG3; 00218 00219 // Compute Hit Density to find MIP events 00220 int NbLayers=0; 00221 for (int layer=0 ; layer < 50 ; layer++){ 00222 if(NbHitLayer[layer]>0) 00223 { 00224 NbLayers++; 00225 } 00226 } 00227 00228 float HitDensity = 1. * nhit_roi/NbLayers; 00229 hDensity->Fill(HitDensity); 00230 00231 //Select Pion events 00232 if( nhit_roi<150 || HitDensity<4 || RG1G3 > 1.0 ) continue; 00233 00234 // Find where the shower starts 00235 00236 NbLayers=0; 00237 int ShowerStart=-10; 00238 int ShowerStartDa=-10; 00239 00240 for (int layer=1 ; layer < 48 ; layer++){ 00241 if(NbHitLayer[layer]>0) 00242 { 00243 NbLayers++; 00244 tp_nhit->Fill(layer,NbHitLayer[layer]); 00245 if(ShowerStartDa<0 && NbHitLayer77[layer-1]<=2 && NbHitLayer77[layer]>=2 && NbHitLayer77[layer+1]>=4) {ShowerStartDa=layer; } 00246 if(ShowerStart<0 && NbHitLayer77[layer-1]<=3 && NbHitLayer77[layer+1]-NbHitLayer77[layer]>2 && NbHitLayer77[layer+2]-NbHitLayer77[layer+1]>2) { 00247 //if(ShowerStart<0 && NbHitLayer77[layer-1]<=2 && NbHitLayer77[layer+1]>NbHitLayer77[layer] && NbHitLayer77[layer+2]>NbHitLayer77[layer+1]) { 00248 // if(ShowerStart<0 && NbHitLayer77[layer-1]<=2 && NbHitLayer77[layer]>4) { 00249 ShowerStart = layer; // cout << ShowerStart << endl ; 00250 00251 } 00252 } 00253 } 00254 hShowerStart->Fill(ShowerStart); 00255 hShowerStartDa->Fill(ShowerStartDa); 00256 00257 int uMgaslayer=45; 00258 if(ShowerStart<=uMgaslayer){hLongprof->Fill(uMgaslayer-ShowerStart,NbHitLayer[uMgaslayer]);} 00259 if(ShowerStart<=uMgaslayer){hLongprof1->Fill(uMgaslayer-ShowerStart,NbHitLayer[uMgaslayer]);} 00260 uMgaslayer=9; 00261 if(ShowerStart<=uMgaslayer){hLongprof->Fill(uMgaslayer-ShowerStart,NbHitLayer[uMgaslayer]);} 00262 uMgaslayer=20; 00263 if(ShowerStart<=uMgaslayer){hLongprof->Fill(uMgaslayer-ShowerStart,NbHitLayer[uMgaslayer]);} 00264 uMgaslayer=32; 00265 if(ShowerStart<=uMgaslayer){hLongprof->Fill(uMgaslayer-ShowerStart,NbHitLayer[uMgaslayer]);} 00266 00267 if(ShowerStart==1){ 00268 for (int layer=0 ; layer < 48 ; layer++){ 00269 hLongprofL1->Fill(layer,NbHitLayer[layer]);} 00270 } 00271 00272 if (nbEvt%10000 == 0) { 00273 cout << "Processing event: " << nbEvt 00274 << endl; 00275 } 00276 tf.cd(0); 00277 } // end event 00278 cout << " Finished with events in this TREE " << endl; 00279 deltaTRun->Write(); 00280 delete deltaTRun; 00281 } // else 00282 } // end of tree or run 00283 cout << " Ending file in chain " << endl; 00284 } // end of chain 00285 cout << " Going to finsh the work after" << nbEvt << " events "<< endl; 00286 tf.cd(); 00287 hDensity->Write(); 00288 hxy->Write(); 00289 hz->Write(); 00290 hnhit->Write(); 00291 tp_nhit->Write(); 00292 hnhit_roi->Write(); 00293 hShowerStart->Write(); 00294 hShowerStartDa->Write(); 00295 hLongprof->Write(); 00296 hLongprof1->Write(); 00297 hLongprofL1->Write(); 00298 for (int i=0 ; i<50 ; i++){ 00299 hxyMapTot[i]->Write(); 00300 } 00301 00302 }
| const int maxEvents = 3000000 |
Definition at line 41 of file getShowerStart.cpp.
| const int maxTrees = 100 |
Definition at line 42 of file getShowerStart.cpp.
| const int maxHits = 200 |
Definition at line 43 of file getShowerStart.cpp.
1.4.7