#include <map>#include "TKey.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 <TCanvas.h>#include "TGraphErrors.h"#include "TGraphAsymmErrors.h"#include <TH2F.h>Include dependency graph for ShowerProfile_CH.cpp:
Go to the source code of this file.
Functions | |
| int | main (int argc, char **argv) |
| int main | ( | int | argc, | |
| char ** | argv | |||
| ) |
Definition at line 42 of file ShowerProfile_CH.cpp.
00042 { 00043 00044 if ( argc !=2 ) { 00045 FILE_LOG(logERROR) << "usage: example rootFile " << endl; 00046 exit(1); 00047 } 00048 00049 //User's parameters 00050 string rootName; 00051 rootName.assign(argv[1]); 00052 00053 // example of data file 00054 // /lapp_data/LC/Detecteurs/MicroMegas/data/TB2012/SPS_H2_may_2012/Production/CaloHits/DHCAL_714405_I*_0.slcio.root 00055 00056 char nameofplot[1000]; 00057 sprintf(nameofplot,"./ShowerProfile_714405_DetLayer48.gif"); 00058 00059 //Plotting window 00060 TCanvas* c1 = new TCanvas("c1","Shower Profile", 200, 0, 1200, 1200); 00061 00062 //Histogramms 00063 TH2F* h02DProfx = new TH2F("02DProfx", "Basic Energy 2D Profile;Layer (z);X", 50, 0., 50.,100,0.,100.); 00064 TH2F* h02DProfy = new TH2F("02DProfy", "Basic Energy 2D Profile;Layer (z);Y", 50, 0., 50.,100,0.,100.); 00065 TH1F* h0EProf = new TH1F("0EProf", "Basic Energy Profile;Layer (z); Hits",50, 0., 50.); 00066 TH2F* h2DProfx = new TH2F("2DProfx", "Normal Energy 2D Profile;Layer (z);X", 50, 0., 50.,100,0.,100.); 00067 TH2F* h2DProfy = new TH2F("2DProfy", "Normal Energy 2D Profile;Layer (z);Y", 50, 0., 50.,100,0.,100.); 00068 TH1F* hEProf = new TH1F("EProf", "Normal Energy Profile;Layer (z); Hits",50, 0., 50.); 00069 TH2F* hKs2DProfx = new TH2F("Ks2DProfx", "Micromegas Energy 2D Profile (detector in layer 48);Reconsructed layer (z);X", 50, 0., 50.,100,0.,100.); 00070 TH2F* hKs2DProfy = new TH2F("Ks2DProfy", "Micromegas Energy 2D Profile (detector in layer 48);Reconsructed layer (z);Y", 50, 0., 50.,100,0.,100.); 00071 TH1F* hKsEProf = new TH1F("KsEProf", "Micromegas Energy Profile (detector in layer 48);Reconstructed layer (z); Hits",50, 0., 50.); 00072 00073 //Variables 00074 00075 int nbEvt = 0; 00076 int nbEvt_tot_used = 0; 00077 int nbEvt_ShowerStart[50]; 00078 for(int layer=0;layer<50;layer++) nbEvt_ShowerStart[layer]=0; 00079 int nbEvt_Micromegas = 0; 00080 00081 //first data reading structure 00082 int NbTrees = 0; 00083 TChain chain("T"); 00084 for (int i =0; i<5 ; i++){ 00085 TString SubFile = Form("_I%i",i); 00086 TString fname = Form("/lapp_data/LC/Detecteurs/MicroMegas/data/TB2012/SPS_H2_may_2012/Production/CaloHits/CaloHit_DHCAL_"+rootName+SubFile+"_0.slcio.root"); 00087 chain.Add(fname); 00088 cout << fname << endl; 00089 } 00090 00091 TObjArray *fileElements=chain.GetListOfFiles(); 00092 TIter next(fileElements); 00093 TChainElement *chEl=0; 00094 int NchainFiles=0; 00095 while (( chEl=(TChainElement*)next() )) 00096 { 00097 TFile f(chEl->GetTitle()); 00098 TIter nextkey(f.GetListOfKeys()); 00099 TKey *key; 00100 NchainFiles++; 00101 cout << " Reading a new file in the chain " << NchainFiles << endl; 00102 while (key = (TKey*)nextkey()) 00103 { 00104 TString keyname = key->GetName(); 00105 if ((keyname == "hResidualX") || (keyname == "hResidualY")){} 00106 else 00107 { 00108 TTree *tree = (TTree*)key->ReadObj(); 00109 cout << "---------------------------- NEW TTREE -----------------"<< endl; 00110 NbTrees++; 00111 CaloRun *run=NULL; 00112 00113 try { 00114 CaloRun *run=dynamic_cast<CaloRun *>(tree->GetUserInfo()->First()); 00115 if( run != NULL ){cout << "Run[" << run->GetRunId() << "], temperature[" << run->GetTemperature() << "] pression[" << run->GetPressure() << "] " << endl ; } 00116 00117 } 00118 catch ( ... ) {} 00119 CaloEvent *evt = new CaloEvent(); 00120 TBranch *branch= tree->GetBranch("CaloEvent"); 00121 branch->SetAddress(&evt); 00122 00123 CaloHit* hit =NULL; 00124 int nbEntries = tree->GetEntries(); 00125 cout << " Nb of events to read " << nbEntries << endl; 00126 00127 // A first loop over events to get the number of evts with a shower starting at each layer 00128 for (int evtNum = 0; evtNum < nbEntries; evtNum++){ 00129 00130 tree->GetEntry(evtNum); 00131 00132 //if(evt->IsSquare(400,int asu,int abs_time, int chamber=48)) continue; 00133 //if(evt->IsBorder(80,int asu,int abs_time, int chamber=48)) continue; 00134 //if(evt->Ischip(40,int chip,int abs_time, int chamber=48)) continue; 00135 00136 int nHits = evt->GetNHits(); 00137 double HitX[nHits]; 00138 double HitY[nHits]; 00139 int zSlot[nHits]; 00140 for(int i=0 ; i < nHits ; i++) 00141 { 00142 hit = (CaloHit*)evt->GetHits()->UncheckedAt(i); 00143 HitX[i] = hit->GetX() - 5000 ; 00144 HitY[i] = hit->GetY() - 5000 ; 00145 zSlot[i] = hit->GetLayer(); 00146 } 00147 00148 //Count the hits of the main shower (near the center) at each layer 00149 int MainShowerHitperLayer[50]; 00150 for (int layer = 0; layer < 50; layer++) MainShowerHitperLayer[layer] = 0; 00151 for (int hitNum = 0; hitNum < nHits; hitNum++) if(HitX[hitNum] > 400000 && HitX[hitNum] < 600000 && HitY[hitNum] > 400000 && HitY[hitNum] < 600000) MainShowerHitperLayer[zSlot[hitNum]]++; 00152 00153 //Find shower start of this event 00154 int ShowerStart = -1; 00155 for (int layer = 0 ; layer < 50 ; layer++){ 00156 if (MainShowerHitperLayer[layer] >= 2 && MainShowerHitperLayer[layer+1] >= 4 ){ 00157 ShowerStart=layer; 00158 break; 00159 } 00160 } 00161 00162 //Remember the number of evts with a shower starting at each layer 00163 for(int layer=0;layer<50;layer++) if(ShowerStart == layer) nbEvt_ShowerStart[layer]++; 00164 00165 //Remember the number of evts with a shower starting in the calo 00166 if(ShowerStart>0) nbEvt_tot_used++; 00167 00168 //Remember the number of evts with a shower starting in the "micromegas" interval 00169 if(ShowerStart<=48 && ShowerStart>0) nbEvt_Micromegas++; 00170 00171 }//end first loop over events 00172 }//end else 00173 }//end key 00174 }//end chain 00175 00176 cout<<endl; 00177 00178 //Main data reading structure 00179 NbTrees = 0; 00180 TObjArray *fileElement=chain.GetListOfFiles(); 00181 TIter nex(fileElement); 00182 TChainElement *chElt=0; 00183 NchainFiles=0; 00184 while (( chElt=(TChainElement*)nex() )) 00185 { 00186 TFile f(chElt->GetTitle()); 00187 TIter nextkey(f.GetListOfKeys()); 00188 TKey *key; 00189 NchainFiles++; 00190 while (key = (TKey*)nextkey()) 00191 { 00192 TString keyname = key->GetName(); 00193 if ((keyname == "hResidualX") || (keyname == "hResidualY")){} 00194 else 00195 { 00196 TTree *tree = (TTree*)key->ReadObj(); 00197 NbTrees++; 00198 CaloRun *run=NULL; 00199 00200 try { 00201 CaloRun *run=dynamic_cast<CaloRun *>(tree->GetUserInfo()->First()); 00202 00203 } 00204 catch ( ... ) {} 00205 CaloEvent *evt = new CaloEvent(); 00206 TBranch *branch= tree->GetBranch("CaloEvent"); 00207 branch->SetAddress(&evt); 00208 00209 CaloHit* hit =NULL; 00210 int nbEntries = tree->GetEntries(); 00211 00212 // Main loop over events 00213 int ev=0; 00214 for ( int evtNum = 0; evtNum < nbEntries ; evtNum++) 00215 { 00216 00217 nbEvt++; 00218 tree->GetEntry(evtNum); 00219 00220 //if(evt->IsSquare(800,int asu,int abs_time, int chamber=48)) continue; 00221 //if(evt->IsBorder(80,int asu,int abs_time, int chamber=48)) continue; 00222 //if(evt->Ischip(40,int chip,int abs_time, int chamber=48)) continue; 00223 00224 int nHits = evt->GetNHits(); 00225 double HitX[nHits]; 00226 double HitY[nHits]; 00227 int zSlot[nHits]; 00228 //double HitEng[nHits]; 00229 for(int i=0 ; i < nHits ; i++){ 00230 hit = (CaloHit*)evt->GetHits()->UncheckedAt(i); 00231 HitX[i] = hit->GetX() - 5000 ; 00232 HitY[i] = hit->GetY() - 5000 ; 00233 zSlot[i] = hit->GetLayer(); 00234 //HitEng[i] = hit->GetThreshold(); 00235 //cout<<"i "<<i<<"; X "<<HitX[i]<<"; Y "<<HitY[i]<<"; z "<<zSlot[i]<<"; Z "<<HitX[i]<<"; E "<<HitEng[i]<<endl; 00236 } 00237 00238 //Find the number of hits at each layer 00239 int MainShowerHitperLayer[50]; 00240 int TotalHitperLayer[50]; 00241 for (int layer = 0; layer < 50; layer++) { 00242 MainShowerHitperLayer[layer] = 0; 00243 TotalHitperLayer[layer] = 0; 00244 } 00245 00246 for (int hitNum = 0; hitNum < nHits; hitNum++) { 00247 //Basic 00248 TotalHitperLayer[zSlot[hitNum]]++; 00249 00250 //Count the hits of the main shower (near the center) at each layer 00251 if(HitX[hitNum] > 400000 && HitX[hitNum] < 600000 && HitY[hitNum] > 400000 && HitY[hitNum] < 600000) { 00252 MainShowerHitperLayer[zSlot[hitNum]]++; 00253 } 00254 } 00255 00256 //Find shower start of this event 00257 int ShowerStart = -1; 00258 for (int layer = 0 ; layer < 48 ; layer++){ 00259 if (MainShowerHitperLayer[layer] >= 2 && MainShowerHitperLayer[layer+1] >= 4 ){ 00260 ShowerStart=layer; 00261 break; 00262 } 00263 00264 } 00265 //cout<<"sh "<<ShowerStart<<endl; 00266 00267 // Loop over hits 00268 for (int hit = 0; hit < nHits; hit++) { 00269 int thelayer = zSlot[hit]; 00270 //float Edep = HitEng[hit]; 00271 double x = HitX[hit]/10000.; 00272 double y = HitY[hit]/10000.; 00273 00274 //if(evtNum<10 && hit<10) cout<<"thelayer "<<thelayer<<"; x "<<x<<"; y "<<y<<"; ShowerStart "<<ShowerStart<<"; nbEntries "<<nbEntries<<"; nbEvt_ShowerStart[1] "<<nbEvt_ShowerStart[1]<<"; nbEvt_ShowerStart[ShowerStart]"<<nbEvt_ShowerStart[ShowerStart]<<endl; //<<"; Edep "<<Edep<<endl; 00275 00276 //Basic energy profile: 00277 //All events cumulated 00278 if(ShowerStart>0){ 00279 h02DProfx->Fill(thelayer, x, 100000./nbEvt_tot_used); // I normalize al histos to 100 000 evts 00280 h02DProfy->Fill(thelayer, y, 100000./nbEvt_tot_used); 00281 h0EProf->Fill(thelayer, 100000./nbEvt_tot_used); 00282 } 00283 00284 //Normal method energy profile: 00285 //Events whith ShowerStart=1; 00286 //We move our detector at every layer; 00287 //->we plot in bin 48 the energy deposited in layer 48 for ShoweStart=1, 00288 //->we plot in bin 48-1 the energy deposited in layer 48-1 for ShowerStart=1, 00289 //=>we plot in bin i the energy deposited in layer i, for all events with ShowerStart=1, 00290 //to normalize we divide by the number of events that we have used. 00291 if (ShowerStart == 1){ 00292 h2DProfx->Fill(thelayer,x, 100000./nbEvt_ShowerStart[1]); 00293 h2DProfy->Fill(thelayer,y, 100000./nbEvt_ShowerStart[1]); 00294 hEProf->Fill(thelayer,100000./nbEvt_ShowerStart[1]); 00295 } 00296 00297 //MicroMegas method energy profile: 00298 //We are fixed in the 48th layer; 00299 //->we plot in bin 48 the energy deposited in layer 48 for ShowerStart=1, to normalize we divide by the number of events used: nbEvt_ShowerStart[layer1] 00300 //->we plot in bin 48-1 the energy deposited in layer 48 for ShowerStart=2, to normalize we divide by the number of events used: nbEvt_ShowerStart[layer2] 00301 //=>we plot in bin 48-ShowerStart+1 the energy deposited in layer 48, to normalize we divide by the number of events used: nbEvt_ShowerStart[ShowerStart] 00302 if(thelayer==48 && ShowerStart>0 && ShowerStart<=48){ 00303 hKs2DProfx->Fill(48-ShowerStart+1,x, 100000./nbEvt_ShowerStart[ShowerStart]); 00304 hKs2DProfy->Fill(48-ShowerStart+1,y, 100000./nbEvt_ShowerStart[ShowerStart]); 00305 hKsEProf->Fill(48-ShowerStart+1,100000./nbEvt_ShowerStart[ShowerStart]); 00306 } 00307 } // End of the loop over hits 00308 00309 cout << "Processing event number " 00310 << ev << "\t\r" 00311 << flush; 00312 ev++; 00313 00314 } // End of the loop over events 00315 cout << " Finished with events in this TREE " << endl; 00316 } // End else 00317 } // End Key 00318 } //End Chain 00319 00320 //1 Plot for all events 00321 c1->Divide(3,3); 00322 00323 //Basic histo, all events 00324 c1->cd(1); 00325 char Xtitle[100]; 00326 sprintf(Xtitle,"All "); 00327 char buf[10]; 00328 sprintf(buf,"%d",nbEvt_tot_used); 00329 strcat(Xtitle, buf); 00330 strcat(Xtitle," events. Layer (z)"); 00331 h02DProfx->GetXaxis()->SetTitle(Xtitle); 00332 h02DProfx->Draw("Zcol"); 00333 00334 //Normal histo, all events with ShowerStart=1 00335 c1->cd(2); 00336 //char Xtitle[100]; 00337 sprintf(Xtitle,"All events with showers start at layer 1; Sum of "); 00338 //char buf[10]; 00339 sprintf(buf,"%d",nbEvt_ShowerStart[1]); 00340 strcat(Xtitle, buf); 00341 strcat(Xtitle," events. Layer (z)"); 00342 h2DProfx->GetXaxis()->SetTitle(Xtitle); 00343 h2DProfx->Draw("Zcol"); 00344 00345 //Micromegas histo, all events with ShowerStart<48 00346 c1->cd(3); 00347 sprintf(Xtitle,"All events with ShowerStart<48; Sum of "); 00348 sprintf(buf,"%d",nbEvt_Micromegas); 00349 strcat(Xtitle, buf); 00350 strcat(Xtitle," events. Layer (z)"); 00351 hKs2DProfx->GetXaxis()->SetTitle(Xtitle); 00352 hKs2DProfx->Draw("Zcol"); 00353 00354 //Basic histo, all events 00355 c1->cd(4); 00356 sprintf(Xtitle,"All "); 00357 sprintf(buf,"%d",nbEvt_tot_used); 00358 strcat(Xtitle, buf); 00359 strcat(Xtitle," events. Layer (z)"); 00360 h02DProfy->GetXaxis()->SetTitle(Xtitle); 00361 h02DProfy->Draw("Zcol"); 00362 00363 //Normal histo, all events with ShowerStart=1 00364 c1->cd(5); 00365 sprintf(Xtitle,"All events with showers start at layer 1; Sum of "); 00366 sprintf(buf,"%d",nbEvt_ShowerStart[1]); 00367 strcat(Xtitle, buf); 00368 strcat(Xtitle," events. Layer (z)"); 00369 h2DProfy->GetXaxis()->SetTitle(Xtitle); 00370 h2DProfy->Draw("Zcol"); 00371 00372 //Micromegas histo, all events with ShowerStart<48 00373 c1->cd(6); 00374 sprintf(Xtitle,"All events with ShowerStart<48; Sum of "); 00375 sprintf(buf,"%d",nbEvt_Micromegas); 00376 strcat(Xtitle, buf); 00377 strcat(Xtitle," events. Layer (z)"); 00378 hKs2DProfy->GetXaxis()->SetTitle(Xtitle); 00379 hKs2DProfy->Draw("Zcol"); 00380 00381 //Basic histo, all events 00382 c1->cd(7); 00383 h0EProf->SetLineColor(1); 00384 sprintf(Xtitle,"All "); 00385 sprintf(buf,"%d",nbEvt_tot_used); 00386 strcat(Xtitle, buf); 00387 strcat(Xtitle," events. Layer (z)"); 00388 h0EProf->GetXaxis()->SetTitle(Xtitle); 00389 h0EProf->Draw(); 00390 00391 //Normal histo, all events with ShowerStart=1 00392 c1->cd(8); 00393 hEProf->SetLineColor(1); 00394 sprintf(Xtitle,"All events with showers start at layer 1; Sum of "); 00395 sprintf(buf,"%d",nbEvt_ShowerStart[1]); 00396 strcat(Xtitle, buf); 00397 strcat(Xtitle," events. Layer (z)"); 00398 hEProf->GetXaxis()->SetTitle(Xtitle); 00399 hEProf->Draw(); 00400 00401 //Micromegas histo, all events with ShowerStart<48 00402 c1->cd(9); 00403 hKsEProf->SetLineColor(1); 00404 sprintf(Xtitle,"All events with ShowerStart<48; Sum of "); 00405 sprintf(buf,"%d",nbEvt_Micromegas); 00406 strcat(Xtitle, buf); 00407 strcat(Xtitle," events. Layer (z)"); 00408 hKsEProf->GetXaxis()->SetTitle(Xtitle); 00409 hKsEProf->Draw(); 00410 00411 //save the plot 00412 c1->SaveAs(nameofplot); 00413 00414 cout << " Going to finsh the work after" << nbEvt << " events "<< endl; 00415 }
1.4.7