1 //------------------------------------------------------------------------------
2 // Implementation of AliPerformanceDCA class. It keeps information from
3 // comparison of reconstructed and MC particle tracks. In addtion,
4 // it keeps selection cuts used during comparison. The comparison
5 // information is stored in the ROOT histograms. Analysis of these
6 // histograms can be done by using Analyse() class function. The result of
7 // the analysis (histograms/graphs) are stored in the folder
8 // which is a data member of AliPerformanceDCA.
10 // Author: J.Otwinowski 04/02/2008
11 //------------------------------------------------------------------------------
15 // after running comparison task, read the file, and get component
16 gROOT->LoadMacro("$ALICE_ROOT/PWG1/Macros/LoadMyLibs.C");
18 TFile f("Output.root");
19 AliPerformanceDCA * compObj = (AliPerformanceDCA*)coutput->FindObject("AliPerformanceDCA");
21 // Analyse comparison data
24 // the output histograms/graphs will be stored in the folder "folderDCA"
25 compObj->GetAnalysisFolder()->ls("*");
27 // user can save whole comparison object (or only folder with anlysed histograms)
28 // in the seperate output file (e.g.)
29 TFile fout("Analysed_DCA.root","recreate");
30 compObj->Write(); // compObj->GetAnalysisFolder()->Write();
43 #include "AliPerformanceDCA.h"
44 #include "AliESDEvent.h"
45 #include "AliESDVertex.h"
47 #include "AliMathBase.h"
48 #include "AliRecInfoCuts.h"
49 #include "AliMCInfoCuts.h"
51 #include "AliMCEvent.h"
52 #include "AliTracker.h"
53 #include "AliHeader.h"
54 #include "AliGenEventHeader.h"
58 ClassImp(AliPerformanceDCA)
60 //_____________________________________________________________________________
61 AliPerformanceDCA::AliPerformanceDCA():
62 AliPerformanceObject("AliPerformanceDCA"),
74 // default constructor
78 //_____________________________________________________________________________
79 AliPerformanceDCA::AliPerformanceDCA(Char_t* name="AliPerformanceDCA", Char_t* title="AliPerformanceDCA",Int_t analysisMode=0, Bool_t hptGenerator=kFALSE):
80 AliPerformanceObject(name,title),
94 SetAnalysisMode(analysisMode);
95 SetHptGenerator(hptGenerator);
99 //_____________________________________________________________________________
100 AliPerformanceDCA::~AliPerformanceDCA()
103 if(fDCAHisto) delete fDCAHisto; fDCAHisto=0;
104 if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
107 //_____________________________________________________________________________
108 void AliPerformanceDCA::Init()
112 Double_t ptMin = 1.e-2, ptMax = 10.;
114 Double_t *binsPt = 0;
115 if (IsHptGenerator()) {
116 nPtBins = 100; ptMax = 100.;
117 binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
119 binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
124 Double_t binsPt[32] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.25,2.5,2.75,3.,3.5,4.,5.,6.,8.,10.};
125 Double_t ptMin = 0., ptMax = 10.;
127 if(IsHptGenerator() == kTRUE) {
129 ptMin = 0.; ptMax = 100.;
133 //dca_r, dca_z, eta, pt
134 Int_t binsQA[4] = {100,100,30,nPtBins};
135 Double_t xminQA[4] = {-10.,-10.,-1.5,ptMin};
136 Double_t xmaxQA[4] = {10.,10.,1.5,ptMax};
138 fDCAHisto = new THnSparseF("fDCAHisto","dca_r:dca_z:eta:pt",4,binsQA,xminQA,xmaxQA);
139 fDCAHisto->SetBinEdges(3,binsPt);
141 fDCAHisto->GetAxis(0)->SetTitle("dca_r (cm)");
142 fDCAHisto->GetAxis(1)->SetTitle("dca_z (cm)");
143 fDCAHisto->GetAxis(2)->SetTitle("#eta");
144 fDCAHisto->GetAxis(3)->SetTitle("p_{T} (GeV/c)");
149 AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
151 AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
154 fAnalysisFolder = CreateFolder("folderDCA","Analysis DCA Folder");
157 //_____________________________________________________________________________
158 void AliPerformanceDCA::ProcessTPC(AliStack* const stack, AliESDtrack *const esdTrack)
160 // Fill DCA comparison information
161 if(!esdTrack) return;
163 const AliExternalTrackParam *track = esdTrack->GetTPCInnerParam();
166 Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
167 esdTrack->GetImpactParametersTPC(dca,cov);
169 if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
171 Double_t vDCAHisto[4]={dca[0],dca[1],track->Eta(),track->Pt()};
172 fDCAHisto->Fill(vDCAHisto);
175 // Fill rec vs MC information
181 //_____________________________________________________________________________
182 void AliPerformanceDCA::ProcessTPCITS(AliStack* const stack, AliESDtrack *const esdTrack)
184 // Fill DCA comparison information
185 if(!esdTrack) return;
187 Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
188 esdTrack->GetImpactParameters(dca,cov);
190 if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
191 if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
192 Int_t clusterITS[200];
193 if(esdTrack->GetITSclusters(clusterITS)<fCutsRC->GetMinNClustersITS()) return; // min. nb. ITS clusters
195 Double_t vDCAHisto[4]={dca[0],dca[1],esdTrack->Eta(),esdTrack->Pt()};
196 fDCAHisto->Fill(vDCAHisto);
199 // Fill rec vs MC information
205 void AliPerformanceDCA::ProcessConstrained(AliStack* const /*stack*/, AliESDtrack *const /*esdTrack*/)
207 // Fill DCA comparison information
209 AliDebug(AliLog::kWarning, "Warning: Not implemented");
212 //_____________________________________________________________________________
213 Long64_t AliPerformanceDCA::Merge(TCollection* const list)
215 // Merge list of objects (needed by PROOF)
223 TIterator* iter = list->MakeIterator();
226 // collection of generated histograms
228 while((obj = iter->Next()) != 0)
230 AliPerformanceDCA* entry = dynamic_cast<AliPerformanceDCA*>(obj);
231 if (entry == 0) continue;
233 fDCAHisto->Add(entry->fDCAHisto);
240 //_____________________________________________________________________________
241 void AliPerformanceDCA::Exec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent, AliESDfriend *const esdFriend, const Bool_t bUseMC, const Bool_t bUseESDfriend)
243 // Process comparison information
247 Error("Exec","esdEvent not available");
250 AliHeader* header = 0;
251 AliGenEventHeader* genHeader = 0;
258 Error("Exec","mcEvent not available");
261 // get MC event header
262 header = mcEvent->Header();
264 Error("Exec","Header not available");
268 stack = mcEvent->Stack();
270 Error("Exec","Stack not available");
274 genHeader = header->GenEventHeader();
276 Error("Exec","Could not retrieve genHeader from Header");
279 genHeader->PrimaryVertex(vtxMC);
285 Error("Exec","esdFriend not available");
291 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
293 AliESDtrack *track = esdEvent->GetTrack(iTrack);
296 if(GetAnalysisMode() == 0) ProcessTPC(stack,track);
297 else if(GetAnalysisMode() == 1) ProcessTPCITS(stack,track);
298 else if(GetAnalysisMode() == 2) ProcessConstrained(stack,track);
300 printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
306 //_____________________________________________________________________________
307 void AliPerformanceDCA::Analyse()
310 // Analyse comparison information and store output histograms
311 // in the analysis folder "folderDCA"
314 TH1::AddDirectory(kFALSE);
318 TObjArray *aFolderObj = new TObjArray;
321 // set pt measurable range
322 fDCAHisto->GetAxis(3)->SetRangeUser(0.10,10.);
325 h2D = fDCAHisto->Projection(0,1); // inverse projection convention
326 h2D->SetName("dca_r_vs_dca_z");
327 h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(1)->GetTitle());
328 h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(0)->GetTitle());
329 sprintf(title,"%s vs %s",fDCAHisto->GetAxis(0)->GetTitle(),fDCAHisto->GetAxis(1)->GetTitle());
330 h2D->SetTitle(title);
331 aFolderObj->Add(h2D);
334 h2D = fDCAHisto->Projection(0,2);
335 h2D->SetName("dca_r_vs_eta");
336 h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
337 h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(0)->GetTitle());
338 sprintf(title,"%s vs %s",fDCAHisto->GetAxis(2)->GetTitle(),fDCAHisto->GetAxis(0)->GetTitle());
339 h2D->SetTitle(title);
340 aFolderObj->Add(h2D);
342 h1D = MakeStat1D(h2D,0,0);
343 h1D->SetName("mean_dca_r_vs_eta");
344 h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
345 h1D->GetYaxis()->SetTitle("mean_dca_r (cm)");
346 sprintf(title," mean_dca_r (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
347 h1D->SetTitle(title);
348 aFolderObj->Add(h1D);
350 h1D = MakeStat1D(h2D,0,1);
351 h1D->SetName("rms_dca_r_vs_eta");
352 h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
353 h1D->GetYaxis()->SetTitle("rms_dca_r (cm)");
354 sprintf(title," rms_dca_r (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
355 h1D->SetTitle(title);
356 aFolderObj->Add(h1D);
359 h2D = fDCAHisto->Projection(0,3);
360 h2D->SetName("dca_r_vs_pt");
361 h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
362 h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(0)->GetTitle());
363 sprintf(title,"%s vs %s",fDCAHisto->GetAxis(0)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
364 h2D->SetTitle(title);
365 h2D->SetBit(TH1::kLogX);
366 aFolderObj->Add(h2D);
368 h1D = MakeStat1D(h2D,0,0);
369 h1D->SetName("mean_dca_r_vs_pt");
370 h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
371 h1D->GetYaxis()->SetTitle("mean_dca_r (cm)");
372 sprintf(title,"mean_dca_r (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
373 h1D->SetTitle(title);
374 h1D->SetBit(TH1::kLogX);
375 aFolderObj->Add(h1D);
377 h1D = MakeStat1D(h2D,0,1);
378 h1D->SetName("rms_dca_r_vs_pt");
379 h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
380 h1D->GetYaxis()->SetTitle("rms_dca_r (cm)");
381 sprintf(title,"rms_dca_r (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
382 h1D->SetTitle(title);
383 h1D->SetBit(TH1::kLogX);
384 aFolderObj->Add(h1D);
387 h2D = fDCAHisto->Projection(1,2);
388 h2D->SetName("dca_z_vs_eta");
389 h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
390 h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(1)->GetTitle());
391 sprintf(title,"%s vs %s",fDCAHisto->GetAxis(1)->GetTitle(),fDCAHisto->GetAxis(2)->GetTitle());
392 h2D->SetTitle(title);
393 aFolderObj->Add(h2D);
395 h1D = MakeStat1D(h2D,0,0);
396 h1D->SetName("mean_dca_z_vs_eta");
397 h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
398 h1D->GetYaxis()->SetTitle("mean_dca_z (cm)");
399 sprintf(title,"mean_dca_z (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
400 h1D->SetTitle(title);
401 aFolderObj->Add(h1D);
403 h1D = MakeStat1D(h2D,0,1);
404 h1D->SetName("rms_dca_z_vs_eta");
405 h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
406 h1D->GetYaxis()->SetTitle("rms_dca_z (cm)");
407 sprintf(title,"rms_dca_z (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
408 h1D->SetTitle(title);
409 aFolderObj->Add(h1D);
412 h2D = fDCAHisto->Projection(1,3);
413 h2D->SetName("dca_z_vs_pt");
414 h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
415 h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(1)->GetTitle());
416 sprintf(title,"%s vs %s",fDCAHisto->GetAxis(1)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
417 h2D->SetTitle(title);
418 h2D->SetBit(TH1::kLogX);
419 aFolderObj->Add(h2D);
421 h1D = MakeStat1D(h2D,0,0);
422 h1D->SetName("mean_dca_z_vs_pt");
423 h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
424 h1D->GetYaxis()->SetTitle("mean_dca_z (cm)");
425 sprintf(title,"mean_dca_z (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
426 h1D->SetTitle(title);
427 h1D->SetBit(TH1::kLogX);
428 aFolderObj->Add(h1D);
430 h1D = MakeStat1D(h2D,0,1);
431 h1D->SetName("rms_dca_z_vs_pt");
432 h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
433 h1D->GetYaxis()->SetTitle("rms_dca_z (cm)");
434 sprintf(title,"rms_dca_z (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
435 h1D->SetTitle(title);
436 h1D->SetBit(TH1::kLogX);
437 aFolderObj->Add(h1D);
440 h3D = fDCAHisto->Projection(2,3,0); // normal 3D projection convention
441 h3D->SetName("dca_r_vs_eta_vs_pt");
443 h2D = MakeStat2D(h3D,0,0,0);
444 h2D->SetName("mean_dca_r_vs_eta_vs_pt");
445 h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
446 h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
447 h2D->GetZaxis()->SetTitle("mean_dca_r (cm)");
448 sprintf(title,"mean_dca_r (cm) vs %s vs %s",fDCAHisto->GetAxis(2)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
449 h2D->SetTitle(title);
450 aFolderObj->Add(h2D);
452 h2D = MakeStat2D(h3D,0,0,1);
453 h2D->SetName("rms_dca_r_vs_eta_vs_pt");
454 h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
455 h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
456 h2D->GetZaxis()->SetTitle("rms_dca_r (cm)");
457 sprintf(title,"rms_dca_r (cm) vs %s vs %s",fDCAHisto->GetAxis(2)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
458 h2D->SetTitle(title);
459 aFolderObj->Add(h2D);
462 h3D = fDCAHisto->Projection(2,3,1);
463 h3D->SetName("dca_z_vs_eta_vs_pt");
465 h2D = MakeStat2D(h3D,0,0,0);
466 h2D->SetName("mean_dca_z_vs_eta_vs_pt");
467 h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
468 h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
469 h2D->GetZaxis()->SetTitle("mean_dca_z (cm)");
470 sprintf(title,"mean_dca_z (cm) vs %s vs %s",fDCAHisto->GetAxis(2)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
471 h2D->SetTitle(title);
472 aFolderObj->Add(h2D);
474 h2D = MakeStat2D(h3D,0,0,1);
475 h2D->SetName("rms_dca_z_vs_eta_vs_pt");
476 h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
477 h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
478 h2D->GetZaxis()->SetTitle("rms_dca_z (cm)");
479 sprintf(title,"rms_dca_z (cm) vs %s vs %s",fDCAHisto->GetAxis(2)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
480 h2D->SetTitle(title);
481 aFolderObj->Add(h2D);
490 // export objects to analysis folder
491 fAnalysisFolder = ExportToFolder(aFolderObj);
493 // delete only TObjArray
494 if(aFolderObj) delete aFolderObj;
497 //_____________________________________________________________________________
498 TH1F* AliPerformanceDCA::MakeStat1D(TH2 *hist, Int_t delta0, Int_t type)
500 // Return TH1F histogram
501 // delta - number of bins to integrate
502 // with mean (type == 0) or RMS (type==1)
505 const char* suffix = "_stat1d";
506 sprintf(hname,"%s%s",hist->GetName(),suffix);
507 TAxis* xaxis = hist->GetXaxis();
508 Int_t nbinx = xaxis->GetNbins();
510 TH1F *hnew = (TH1F*)hist->ProjectionX()->Clone();
511 hnew->SetName(hname);
514 for (Int_t ix=0; ix<=nbinx;ix++) {
515 sprintf(name,"%s_%d",hist->GetName(),ix);
516 TH1 *projection = hist->ProjectionY(name,ix-delta0,ix+delta0);
518 Float_t stat= 0., stat_err =0.;
519 if (type==0) { stat = projection->GetMean(); stat_err = projection->GetMeanError(); }
520 if (type==1) { stat = projection->GetRMS(); stat_err = projection->GetRMSError(); }
522 hnew->SetBinContent(ix, stat);
523 hnew->SetBinError(ix, stat_err);
529 //_____________________________________________________________________________
530 TH2F* AliPerformanceDCA::MakeStat2D(TH3 *hist, Int_t delta0, Int_t delta1, Int_t type)
532 // Return TH1F histogram
533 // delta0 - number of bins to integrate in x
534 // delta1 - number of bins to integrate in y
535 // with mean (type==0) or RMS (type==1)
538 const char* suffix = "_stat2d";
539 sprintf(hname,"%s%s",hist->GetName(),suffix);
541 TAxis* xaxis = hist->GetXaxis();
542 Int_t nbinx = xaxis->GetNbins();
544 TH2F *hnew = (TH2F*)hist->Project3D("yx")->Clone();
545 hnew->SetName(hname);
547 TAxis* yaxis = hist->GetYaxis();
548 Int_t nbiny = yaxis->GetNbins();
551 for (Int_t ix=0; ix<=nbinx;ix++) {
552 for (Int_t iy=0; iy<=nbiny;iy++) {
553 sprintf(name,"%s_%d_%d",hist->GetName(),ix,iy);
554 TH1 *projection = hist->ProjectionZ(name,ix-delta0,ix+delta0,iy-delta1,iy+delta1);
556 Float_t stat= 0., stat_err =0.;
557 if (type==0) { stat = projection->GetMean(); stat_err = projection->GetMeanError(); }
558 if (type==1) { stat = projection->GetRMS(); stat_err = projection->GetRMSError(); }
560 hnew->SetBinContent(ix,iy,stat);
561 hnew->SetBinError(ix,iy,stat_err);
568 //_____________________________________________________________________________
569 TFolder* AliPerformanceDCA::ExportToFolder(TObjArray * array)
571 // recreate folder avery time and export objects to new one
573 AliPerformanceDCA * comp=this;
574 TFolder *folder = comp->GetAnalysisFolder();
577 TFolder *newFolder = 0;
579 Int_t size = array->GetSize();
582 // get name and title from old folder
583 name = folder->GetName();
584 title = folder->GetTitle();
590 newFolder = CreateFolder(name.Data(),title.Data());
591 newFolder->SetOwner();
593 // add objects to folder
595 newFolder->Add(array->At(i));
603 //_____________________________________________________________________________
604 TFolder* AliPerformanceDCA::CreateFolder(TString name,TString title) {
605 // create folder for analysed histograms
607 folder = new TFolder(name.Data(),title.Data());