#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 <TF1.h>Include dependency graph for residual.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 40 of file residual.cpp.
00040 { 00041 00042 if ( argc !=2 ) { 00043 FILE_LOG(logERROR) << "usage: example rootFile " << endl; 00044 exit(1); 00045 } 00046 00047 00048 int mip = 0; 00049 int penetrating_mip = 0; 00050 float ngoodchi2 = 0; 00051 float ntest48 = 0; 00052 float neff48 = 0; 00053 00054 TString name,title; 00055 00056 TF1 * fzx = new TF1("fzx","pol1",0,60); 00057 TF1 * fzy = new TF1("fzy","pol1",0,60); 00058 00059 TH1I * hx_evt = new TH1I("hx_evt","Horizontal distribution;x (pad number)",1000000,0,1000000); 00060 TH1I * hy_evt = new TH1I("hy_evt","Vertical distribution;y (pad number)",1000000,0,1000000); 00061 00062 TH1I * hx_all = new TH1I("hx_all","Horizontal distribution;x (pad number)",1000000,0,1000000); 00063 TH1I * hy_all = new TH1I("hy_all","Vertical distribution;y (pad number)",1000000,0,1000000); 00064 00065 TH1F * heff = new TH1F("heff","Efficiency per layer",60,0,60); 00066 TH1F * hmult = new TH1F("hmult","Multiplicity per layer",60,0,60); 00067 00068 TH1I * hnhitcalo = new TH1I("hnhitcalo","Number of hits - All calo",1000,0,1000); 00069 TH1I * hdtcalo = new TH1I("hdtcalo","Time to event time - All calo",41,-20,20); 00070 TH1F * hhit_dens = new TH1F("hhit_dens","Number of hit per layer",100,0,50); 00071 00072 TH1F * hzx_slope = new TH1F("hzx_slope","horizontal slope",101,-3,3); 00073 TH1F * hzy_slope = new TH1F("hzy_slope","vertical slope",101,-3,3); 00074 00075 TH1F * hzx_const = new TH1F("hzx_const","horizontal const",100,0,100); 00076 TH1F * hzy_const = new TH1F("hzy_const","vertical const",100,0,100); 00077 00078 TH1F * hzx_chi2 = new TH1F("hzx_chi2","horizontal chi2/ndf",1000,0,100); 00079 TH1F * hzy_chi2 = new TH1F("hzy_chi2","vertical chi2/ndf",1000,0,100); 00080 00081 TH2I * hxy[50]; 00082 TH1I * hdt[50]; 00083 TH1I * hnhit[50]; 00084 TH1F * hresx[50]; 00085 TH1F * hresy[50]; 00086 int nhit[50]; 00087 float neff[50]; 00088 00089 for (int i=0;i<50;i++) 00090 { 00091 nhit[i] = 0; 00092 00093 name = Form("hdt_%i",i); 00094 title = Form("Time to event time - chb %i",i); 00095 hdt[i]= new TH1I(name,title,41,-20,20); 00096 00097 name = Form("hresx_%i",i); 00098 title = Form("Horizontal residuals - chb %i",i); 00099 hresx[i]= new TH1F(name,title,101,-500000,500000); 00100 00101 name = Form("hresy_%i",i); 00102 title = Form("Vertical residuals - chb %i",i); 00103 hresy[i]= new TH1F(name,title,101,-500000,500000); 00104 00105 name = Form("hnhit_%i",i); 00106 title = Form("Number of CaloHit - chb %i",i); 00107 hnhit[i]= new TH1I(name,title,100,0,100); 00108 00109 name = Form("hxy_%i",i); 00110 title = Form("CaloHit position - chb %i",i); 00111 hxy[i] = new TH2I(name,title,100,0.,1000000.,100,0.,1000000.); 00112 00113 00114 } 00115 00116 00117 00118 00119 00120 00121 string rootName; 00122 rootName.assign(argv[1]); 00123 00124 TChain chain("T"); 00125 for (int i=0 ; i<6 ; i++) 00126 { 00127 if ( i != 5 ) continue; 00128 TString SubFile = Form("_I%i",i); 00129 TString fname = Form("/lapp_data/LC/data/TB2012/SPS_H2_may_2012/Production/CaloHits/CaloHit_DHCAL_"+rootName+SubFile+"_0.slcio.root"); 00130 //TString fname = Form("/lapp_data/LC/data/TB2012/SPS_H2_may_2012/Production/CaloHits/Calo_HCAL_714408"+SubFile+"_0.slcio.root"); 00131 chain.Add(fname); 00132 cout << fname << endl; 00133 } 00134 00135 TObjArray *fileElements=chain.GetListOfFiles(); 00136 TIter next(fileElements); 00137 TChainElement *chEl=0; 00138 int NchainFiles=0; 00139 00140 while (( chEl=(TChainElement*)next() )) 00141 { 00142 TFile f(chEl->GetTitle()); 00143 TIter nextkey(f.GetListOfKeys()); 00144 TKey *key; 00145 00146 int NbTrees = 0; 00147 00148 while (key = (TKey*)nextkey()) 00149 { 00150 00151 TString keyname = key->GetName(); 00152 if ((keyname == "hResidualX") || (keyname == "hResidualY")){} 00153 else{ 00154 00155 TTree *tree = (TTree*)key->ReadObj(); 00156 cout << "---------------------------- NEW TTREE -----------------"<< endl; 00157 int nbEvt =0; 00158 int nbHit = 0; 00159 NbTrees++; 00160 CaloRun *run=NULL; 00161 00162 try 00163 { 00164 CaloRun *run=dynamic_cast<CaloRun *>(tree->GetUserInfo()->First()); 00165 if( run != NULL ) 00166 { 00167 cout << "Run[" << run->GetRunId() << "], temperature[" << run->GetTemperature() << "] pression[" << run->GetPressure() << "] " << endl ; 00168 } 00169 } 00170 00171 catch ( ... ) {} 00172 00173 int RunId=1; 00174 CaloEvent *evt = new CaloEvent(); 00175 TBranch *branch= tree->GetBranch("CaloEvent"); 00176 branch->SetAddress(&evt); 00177 00178 CaloHit* hit =NULL; 00179 int nbEntries = tree->GetEntries(); 00180 00181 // nbEntries = 30; 00182 for ( int evtNum = 0; evtNum < nbEntries ; evtNum++) 00183 { 00184 TGraph * tgzx = new TGraph(); 00185 TGraph * tgzy = new TGraph(); 00186 00187 tree->GetEntry(evtNum); 00188 nbEvt++; 00189 00190 int nbHits = evt->GetNHits(); 00191 int maxDeltat = evt->GetDeltaTmax(); 00192 00193 cout<<"evt "<<evtNum<<" -> nhit = "<<nbHits<<" \r"<<flush; 00194 00195 hnhitcalo->Fill(nbHits); 00196 00197 for(int i=0 ; i < nbHits ; i++) 00198 { 00199 hit = (CaloHit*)evt->GetHits()->UncheckedAt(i); 00200 00201 int deltat = hit->GetDeltaT(); 00202 int xpos = hit->GetX();// * 1e-4; 00203 int ypos = hit->GetY();// * 1e-4; 00204 int zSlot = hit->GetLayer() ; 00205 00206 hdtcalo->Fill(maxDeltat-deltat); 00207 hdt[zSlot]->Fill(maxDeltat-deltat); 00208 00209 if (zSlot >= 48) 00210 { 00211 if (fabs(maxDeltat-deltat) <= 2) 00212 { 00213 nhit[zSlot]++; 00214 hx_evt->Fill(xpos); 00215 hy_evt->Fill(ypos); 00216 } 00217 } 00218 else 00219 { 00220 if (fabs(maxDeltat-deltat) <= 1) 00221 { 00222 nhit[zSlot]++; 00223 hx_evt->Fill(xpos); 00224 hy_evt->Fill(ypos); 00225 } 00226 } 00227 } 00228 00229 int bin_max_x = hx_evt->GetMaximumBin(); 00230 int bin_max_y = hy_evt->GetMaximumBin(); 00231 00232 00233 for(int i=0 ; i < nbHits ; i++) 00234 { 00235 hit = (CaloHit*)evt->GetHits()->UncheckedAt(i); 00236 00237 int deltat = hit->GetDeltaT(); 00238 int xpos = hit->GetX();// * 1e-4; 00239 int ypos = hit->GetY();// * 1e-4; 00240 int zSlot = hit->GetLayer() ; 00241 bool close = 0; 00242 int close_cut = 5; 00243 00244 if (zSlot >= 48) 00245 { 00246 //inverse DIF 00247 //if (xpos<=32){xpos+=64;} 00248 //else{if (xpos>64){xpos-=64;}} 00249 00250 if (fabs(maxDeltat-deltat) <= 2) 00251 { 00252 //hxy[zSlot]->Fill(xpos,ypos); 00253 close = (fabs(bin_max_x-1 - xpos) <= close_cut) && (fabs(bin_max_y-1 - ypos) <= close_cut); 00254 00255 if (nhit[zSlot] <= 4 && close) 00256 { 00257 tgzx->SetPoint(tgzx->GetN(),zSlot+1,xpos); 00258 tgzy->SetPoint(tgzy->GetN(),zSlot+1,ypos); 00259 } 00260 } 00261 } 00262 else 00263 { 00264 if (fabs(maxDeltat-deltat) <= 1) 00265 { 00266 //hxy[zSlot]->Fill(xpos,ypos); 00267 close = (fabs(bin_max_x-1 - xpos) <= close_cut) && (fabs(bin_max_y-1 - ypos) <= close_cut); 00268 00269 if (nhit[zSlot] <= 4 && close) 00270 { 00271 tgzx->SetPoint(tgzx->GetN(),zSlot+1,xpos); 00272 tgzy->SetPoint(tgzy->GetN(),zSlot+1,ypos); 00273 } 00274 } 00275 } 00276 } 00277 00278 int nhot_layer = 0; 00279 int nhit_tot = 0; 00280 int nfirst_ten = 0; 00281 int nlast_ten = 0; 00282 00283 for (int chb=0;chb<50;chb++) 00284 { 00285 hnhit[chb]->Fill(nhit[chb]); 00286 00287 if (nhit[chb]!=0) 00288 { 00289 if (chb<10){nfirst_ten++;} 00290 if (chb>=35 && chb<45){nlast_ten++;} 00291 00292 nhot_layer++; 00293 nhit_tot += nhit[chb]; 00294 } 00295 } 00296 00297 00298 float hit_density = 1. * nhit_tot / nhot_layer; 00299 hhit_dens->Fill(hit_density); 00300 00301 00302 if (hit_density <= 4) 00303 { 00304 mip++; 00305 00306 if (nfirst_ten>=7 && nlast_ten>=7) 00307 { 00308 penetrating_mip++; 00309 00310 fzx->SetParameters(bin_max_x,0); 00311 fzy->SetParameters(bin_max_y,0); 00312 00313 tgzx->Fit(fzx,"q"); 00314 tgzy->Fit(fzy,"q"); 00315 00316 float chi2zx = fzx->GetChisquare()/fzx->GetNDF(); 00317 float chi2zy = fzy->GetChisquare()/fzy->GetNDF(); 00318 00319 hzx_chi2->Fill(chi2zx); 00320 hzy_chi2->Fill(chi2zy); 00321 00322 hzx_slope->Fill(fzx->GetParameter(1)); 00323 hzy_slope->Fill(fzy->GetParameter(1)); 00324 00325 hzx_const->Fill(fzx->GetParameter(0)); 00326 hzy_const->Fill(fzy->GetParameter(0)); 00327 00328 if (chi2zx<=1 && chi2zy<=1) 00329 { 00330 ngoodchi2++; 00331 00332 for(int i=0 ; i < nbHits ; i++) 00333 { 00334 hit = (CaloHit*)evt->GetHits()->UncheckedAt(i); 00335 00336 int deltat = hit->GetDeltaT(); 00337 int xpos = hit->GetX();// * 1e-4; 00338 int ypos = hit->GetY();// * 1e-4; 00339 int zSlot = hit->GetLayer() ; 00340 00341 if (zSlot >= 48) 00342 { 00343 //inverse DIF 00344 //if (xpos<=32){xpos+=64;} 00345 //else{if (xpos>64){xpos-=64;}} 00346 00347 if (fabs(maxDeltat-deltat) <= 2) 00348 { 00349 if (nhit[zSlot]<=4) 00350 { 00351 hxy[zSlot]->Fill(xpos,ypos); 00352 00353 hresx[zSlot]->Fill(fzx->Eval(zSlot+1) - xpos); 00354 hresy[zSlot]->Fill(fzy->Eval(zSlot+1) - ypos); 00355 } 00356 } 00357 } 00358 else 00359 { 00360 if (fabs(maxDeltat-deltat) <= 1) 00361 { 00362 if (nhit[zSlot]<=4) 00363 { 00364 //cout << "zSlot[" << zSlot << "]" << endl; 00365 //cout << "xpos[" << xpos << "]" << endl; 00366 //cout << "ypos[" << ypos << "]" << endl; 00367 //cout << "hxy[zSlot]" << hxy[zSlot] << "]" << endl; 00368 hxy[zSlot]->Fill(xpos,ypos); 00369 hresx[zSlot]->Fill(fzx->Eval(zSlot+1) - xpos); 00370 hresy[zSlot]->Fill(fzy->Eval(zSlot+1) - ypos); 00371 } 00372 } 00373 } 00374 } 00375 00376 for (int chb=0;chb<50;chb++) 00377 { 00378 if (nhit[chb]!=0) 00379 { 00380 neff[chb]++; 00381 hmult->Fill(chb,nhit[chb]); 00382 } 00383 } 00384 00385 if (nhit[49]!=0) 00386 { 00387 ntest48++; 00388 if (nhit[48]!=0) 00389 { 00390 neff48++; 00391 } 00392 } 00393 00394 00395 } 00396 00397 name = Form("tgzx_%i",evtNum); 00398 title = Form("Hit in zx evt %i",evtNum); 00399 tgzx->SetName(name); 00400 tgzx->SetTitle(title); 00401 00402 name = Form("tgzy_%i",evtNum); 00403 title = Form("Hit in zy evt %i",evtNum); 00404 tgzy->SetName(name); 00405 tgzy->SetTitle(title); 00406 00407 tgzx->SetMarkerStyle(20); 00408 tgzy->SetMarkerStyle(20); 00409 00410 00411 00412 for (int b=0;b<1000000;b++) 00413 { 00414 int content = hx_evt->GetBinContent(b+1); 00415 hx_all->Fill(b+1 - bin_max_x,content); 00416 00417 content = hy_evt->GetBinContent(b+1); 00418 hy_all->Fill(b+1 - bin_max_y,content); 00419 } 00420 } 00421 } 00422 00423 for (int chb=0;chb<50;chb++) 00424 { 00425 nhit[chb] = 0; 00426 } 00427 00428 hx_evt->Reset(); 00429 hy_evt->Reset(); 00430 00431 tgzx->Delete(); 00432 tgzy->Delete(); 00433 00434 } 00435 } 00436 } 00437 } 00438 00439 cout<<endl; 00440 cout<<"MIP = "<<mip<<endl; 00441 cout<<"penetrating = "<<penetrating_mip<<endl; 00442 cout<<"good chi2 = "<<ngoodchi2<<endl; 00443 cout<<endl; 00444 00445 for (int chb=0;chb<50;chb++) 00446 { 00447 float eff = neff[chb]/(1.*ngoodchi2); 00448 float eff_err = eff*(1-eff)/sqrt(1.*ngoodchi2); 00449 00450 heff->SetBinContent(chb+1,eff); 00451 heff->SetBinError(chb+1,eff_err); 00452 00453 hmult->SetBinContent(chb+1,hmult->GetBinContent(chb+1)/neff[chb]); 00454 00455 cout<<"Chb "<<chb+1<<" eff = "<<eff<<" +/- "<<eff_err<<" (Nevt = "<<ngoodchi2<<")"<<endl; 00456 } 00457 00458 float eff48 = neff48/ntest48; 00459 float eff48_err = eff48 * (1-eff48) / sqrt(ntest48); 00460 00461 cout<<endl; 00462 cout<<"Chamber 49 (uM#2) = "<<neff48<<" / "<<ntest48<<" = "<<eff48<<" +/- "<<eff48_err<<endl; 00463 cout<<endl; 00464 00465 00466 00467 00468 TFile tf("Muons100GeV_deadzone.root","RECREATE"); 00469 00470 hx_all->Write(); 00471 hy_all->Write(); 00472 00473 heff->Write(); 00474 hmult->Write(); 00475 00476 00477 hzy_slope->Write(); 00478 00479 hzx_const->Write(); 00480 hzy_const->Write(); 00481 00482 hzx_chi2->Write(); 00483 hzy_chi2->Write(); 00484 00485 hhit_dens->Write(); 00486 hdtcalo->Write(); 00487 hnhitcalo->Write(); 00488 00489 for (int i=0;i<50;i++) 00490 { 00491 hdt[i]->Write(); 00492 hxy[i]->Write(); 00493 hnhit[i]->Write(); 00494 hresx[i]->Write(); 00495 hresy[i]->Write(); 00496 } 00497 }
1.4.7