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, const Bool_t bUseMC)
243 // Process comparison information
246 AliDebug(AliLog::kError, "esdEvent not available");
249 AliHeader* header = 0;
250 AliGenEventHeader* genHeader = 0;
257 AliDebug(AliLog::kError, "mcEvent not available");
260 // get MC event header
261 header = mcEvent->Header();
263 AliDebug(AliLog::kError, "Header not available");
267 stack = mcEvent->Stack();
269 AliDebug(AliLog::kError, "Stack not available");
273 genHeader = header->GenEventHeader();
275 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
278 genHeader->PrimaryVertex(vtxMC);
283 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
285 AliESDtrack *track = esdEvent->GetTrack(iTrack);
288 if(GetAnalysisMode() == 0) ProcessTPC(stack,track);
289 else if(GetAnalysisMode() == 1) ProcessTPCITS(stack,track);
290 else if(GetAnalysisMode() == 2) ProcessConstrained(stack,track);
292 printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
298 //_____________________________________________________________________________
299 void AliPerformanceDCA::Analyse()
302 // Analyse comparison information and store output histograms
303 // in the analysis folder "folderDCA"
306 TH1::AddDirectory(kFALSE);
310 TObjArray *aFolderObj = new TObjArray;
313 // set pt measurable range
314 fDCAHisto->GetAxis(3)->SetRangeUser(0.10,10.);
317 h2D = fDCAHisto->Projection(0,1); // inverse projection convention
318 h2D->SetName("dca_r_vs_dca_z");
319 h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(1)->GetTitle());
320 h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(0)->GetTitle());
321 sprintf(title,"%s vs %s",fDCAHisto->GetAxis(0)->GetTitle(),fDCAHisto->GetAxis(1)->GetTitle());
322 h2D->SetTitle(title);
323 aFolderObj->Add(h2D);
326 h2D = fDCAHisto->Projection(0,2);
327 h2D->SetName("dca_r_vs_eta");
328 h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
329 h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(0)->GetTitle());
330 sprintf(title,"%s vs %s",fDCAHisto->GetAxis(2)->GetTitle(),fDCAHisto->GetAxis(0)->GetTitle());
331 h2D->SetTitle(title);
332 aFolderObj->Add(h2D);
334 h1D = MakeStat1D(h2D,0,0);
335 h1D->SetName("mean_dca_r_vs_eta");
336 h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
337 h1D->GetYaxis()->SetTitle("mean_dca_r (cm)");
338 sprintf(title," mean_dca_r (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
339 h1D->SetTitle(title);
340 aFolderObj->Add(h1D);
342 h1D = MakeStat1D(h2D,0,1);
343 h1D->SetName("rms_dca_r_vs_eta");
344 h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
345 h1D->GetYaxis()->SetTitle("rms_dca_r (cm)");
346 sprintf(title," rms_dca_r (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
347 h1D->SetTitle(title);
348 aFolderObj->Add(h1D);
351 h2D = fDCAHisto->Projection(0,3);
352 h2D->SetName("dca_r_vs_pt");
353 h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
354 h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(0)->GetTitle());
355 sprintf(title,"%s vs %s",fDCAHisto->GetAxis(0)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
356 h2D->SetTitle(title);
357 h2D->SetBit(TH1::kLogX);
358 aFolderObj->Add(h2D);
360 h1D = MakeStat1D(h2D,0,0);
361 h1D->SetName("mean_dca_r_vs_pt");
362 h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
363 h1D->GetYaxis()->SetTitle("mean_dca_r (cm)");
364 sprintf(title,"mean_dca_r (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
365 h1D->SetTitle(title);
366 h1D->SetBit(TH1::kLogX);
367 aFolderObj->Add(h1D);
369 h1D = MakeStat1D(h2D,0,1);
370 h1D->SetName("rms_dca_r_vs_pt");
371 h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
372 h1D->GetYaxis()->SetTitle("rms_dca_r (cm)");
373 sprintf(title,"rms_dca_r (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
374 h1D->SetTitle(title);
375 h1D->SetBit(TH1::kLogX);
376 aFolderObj->Add(h1D);
379 h2D = fDCAHisto->Projection(1,2);
380 h2D->SetName("dca_z_vs_eta");
381 h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
382 h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(1)->GetTitle());
383 sprintf(title,"%s vs %s",fDCAHisto->GetAxis(1)->GetTitle(),fDCAHisto->GetAxis(2)->GetTitle());
384 h2D->SetTitle(title);
385 aFolderObj->Add(h2D);
387 h1D = MakeStat1D(h2D,0,0);
388 h1D->SetName("mean_dca_z_vs_eta");
389 h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
390 h1D->GetYaxis()->SetTitle("mean_dca_z (cm)");
391 sprintf(title,"mean_dca_z (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
392 h1D->SetTitle(title);
393 aFolderObj->Add(h1D);
395 h1D = MakeStat1D(h2D,0,1);
396 h1D->SetName("rms_dca_z_vs_eta");
397 h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
398 h1D->GetYaxis()->SetTitle("rms_dca_z (cm)");
399 sprintf(title,"rms_dca_z (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
400 h1D->SetTitle(title);
401 aFolderObj->Add(h1D);
404 h2D = fDCAHisto->Projection(1,3);
405 h2D->SetName("dca_z_vs_pt");
406 h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
407 h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(1)->GetTitle());
408 sprintf(title,"%s vs %s",fDCAHisto->GetAxis(1)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
409 h2D->SetTitle(title);
410 h2D->SetBit(TH1::kLogX);
411 aFolderObj->Add(h2D);
413 h1D = MakeStat1D(h2D,0,0);
414 h1D->SetName("mean_dca_z_vs_pt");
415 h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
416 h1D->GetYaxis()->SetTitle("mean_dca_z (cm)");
417 sprintf(title,"mean_dca_z (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
418 h1D->SetTitle(title);
419 h1D->SetBit(TH1::kLogX);
420 aFolderObj->Add(h1D);
422 h1D = MakeStat1D(h2D,0,1);
423 h1D->SetName("rms_dca_z_vs_pt");
424 h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
425 h1D->GetYaxis()->SetTitle("rms_dca_z (cm)");
426 sprintf(title,"rms_dca_z (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
427 h1D->SetTitle(title);
428 h1D->SetBit(TH1::kLogX);
429 aFolderObj->Add(h1D);
432 h3D = fDCAHisto->Projection(2,3,0); // normal 3D projection convention
433 h3D->SetName("dca_r_vs_eta_vs_pt");
435 h2D = MakeStat2D(h3D,0,0,0);
436 h2D->SetName("mean_dca_r_vs_eta_vs_pt");
437 h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
438 h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
439 h2D->GetZaxis()->SetTitle("mean_dca_r (cm)");
440 sprintf(title,"mean_dca_r (cm) vs %s vs %s",fDCAHisto->GetAxis(2)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
441 h2D->SetTitle(title);
442 aFolderObj->Add(h2D);
444 h2D = MakeStat2D(h3D,0,0,1);
445 h2D->SetName("rms_dca_r_vs_eta_vs_pt");
446 h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
447 h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
448 h2D->GetZaxis()->SetTitle("rms_dca_r (cm)");
449 sprintf(title,"rms_dca_r (cm) vs %s vs %s",fDCAHisto->GetAxis(2)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
450 h2D->SetTitle(title);
451 aFolderObj->Add(h2D);
454 h3D = fDCAHisto->Projection(2,3,1);
455 h3D->SetName("dca_z_vs_eta_vs_pt");
457 h2D = MakeStat2D(h3D,0,0,0);
458 h2D->SetName("mean_dca_z_vs_eta_vs_pt");
459 h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
460 h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
461 h2D->GetZaxis()->SetTitle("mean_dca_z (cm)");
462 sprintf(title,"mean_dca_z (cm) vs %s vs %s",fDCAHisto->GetAxis(2)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
463 h2D->SetTitle(title);
464 aFolderObj->Add(h2D);
466 h2D = MakeStat2D(h3D,0,0,1);
467 h2D->SetName("rms_dca_z_vs_eta_vs_pt");
468 h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
469 h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
470 h2D->GetZaxis()->SetTitle("rms_dca_z (cm)");
471 sprintf(title,"rms_dca_z (cm) vs %s vs %s",fDCAHisto->GetAxis(2)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
472 h2D->SetTitle(title);
473 aFolderObj->Add(h2D);
482 // export objects to analysis folder
483 fAnalysisFolder = ExportToFolder(aFolderObj);
485 // delete only TObjArray
486 if(aFolderObj) delete aFolderObj;
489 //_____________________________________________________________________________
490 TH1F* AliPerformanceDCA::MakeStat1D(TH2 *hist, Int_t delta0, Int_t type)
492 // Return TH1F histogram
493 // delta - number of bins to integrate
494 // with mean (type == 0) or RMS (type==1)
497 const char* suffix = "_stat1d";
498 sprintf(hname,"%s%s",hist->GetName(),suffix);
499 TAxis* xaxis = hist->GetXaxis();
500 Int_t nbinx = xaxis->GetNbins();
502 TH1F *hnew = (TH1F*)hist->ProjectionX()->Clone();
503 hnew->SetName(hname);
506 for (Int_t ix=0; ix<=nbinx;ix++) {
507 sprintf(name,"%s_%d",hist->GetName(),ix);
508 TH1 *projection = hist->ProjectionY(name,ix-delta0,ix+delta0);
510 Float_t stat= 0., stat_err =0.;
511 if (type==0) { stat = projection->GetMean(); stat_err = projection->GetMeanError(); }
512 if (type==1) { stat = projection->GetRMS(); stat_err = projection->GetRMSError(); }
514 hnew->SetBinContent(ix, stat);
515 hnew->SetBinError(ix, stat_err);
521 //_____________________________________________________________________________
522 TH2F* AliPerformanceDCA::MakeStat2D(TH3 *hist, Int_t delta0, Int_t delta1, Int_t type)
524 // Return TH1F histogram
525 // delta0 - number of bins to integrate in x
526 // delta1 - number of bins to integrate in y
527 // with mean (type==0) or RMS (type==1)
530 const char* suffix = "_stat2d";
531 sprintf(hname,"%s%s",hist->GetName(),suffix);
533 TAxis* xaxis = hist->GetXaxis();
534 Int_t nbinx = xaxis->GetNbins();
536 TH2F *hnew = (TH2F*)hist->Project3D("yx")->Clone();
537 hnew->SetName(hname);
539 TAxis* yaxis = hist->GetYaxis();
540 Int_t nbiny = yaxis->GetNbins();
543 for (Int_t ix=0; ix<=nbinx;ix++) {
544 for (Int_t iy=0; iy<=nbiny;iy++) {
545 sprintf(name,"%s_%d_%d",hist->GetName(),ix,iy);
546 TH1 *projection = hist->ProjectionZ(name,ix-delta0,ix+delta0,iy-delta1,iy+delta1);
548 Float_t stat= 0., stat_err =0.;
549 if (type==0) { stat = projection->GetMean(); stat_err = projection->GetMeanError(); }
550 if (type==1) { stat = projection->GetRMS(); stat_err = projection->GetRMSError(); }
552 hnew->SetBinContent(ix,iy,stat);
553 hnew->SetBinError(ix,iy,stat_err);
560 //_____________________________________________________________________________
561 TFolder* AliPerformanceDCA::ExportToFolder(TObjArray * array)
563 // recreate folder avery time and export objects to new one
565 AliPerformanceDCA * comp=this;
566 TFolder *folder = comp->GetAnalysisFolder();
569 TFolder *newFolder = 0;
571 Int_t size = array->GetSize();
574 // get name and title from old folder
575 name = folder->GetName();
576 title = folder->GetTitle();
582 newFolder = CreateFolder(name.Data(),title.Data());
583 newFolder->SetOwner();
585 // add objects to folder
587 newFolder->Add(array->At(i));
595 //_____________________________________________________________________________
596 TFolder* AliPerformanceDCA::CreateFolder(TString name,TString title) {
597 // create folder for analysed histograms
599 folder = new TFolder(name.Data(),title.Data());