1 //------------------------------------------------------------------------------
2 // Implementation of AliComparisonDCA 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 AliComparisonDCA.
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 //AliComparisonDCA * compObj = (AliComparisonDCA*)f.Get("AliComparisonDCA");
20 AliComparisonDCA * compObj = (AliComparisonDCA*)cOutput->FindObject("AliComparisonDCA");
22 // Analyse comparison data
25 // the output histograms/graphs will be stored in the folder "folderDCA"
26 compObj->GetAnalysisFolder()->ls("*");
28 // user can save whole comparison object (or only folder with anlysed histograms)
29 // in the seperate output file (e.g.)
30 TFile fout("Analysed_DCA.root","recreate");
31 compObj->Write(); // compObj->GetAnalysisFolder()->Write();
42 #include "AliComparisonDCA.h"
43 #include "AliESDEvent.h"
44 #include "AliESDRecInfo.h"
45 #include "AliESDVertex.h"
47 #include "AliMCInfo.h"
48 #include "AliMCInfoCuts.h"
49 #include "AliMathBase.h"
50 #include "AliRecInfoCuts.h"
51 #include "AliTracker.h"
55 ClassImp(AliComparisonDCA)
57 //_____________________________________________________________________________
58 AliComparisonDCA::AliComparisonDCA():
59 AliComparisonObject("AliComparisonDCA"),
79 // default constructor
82 //_____________________________________________________________________________
83 AliComparisonDCA::AliComparisonDCA(Char_t* name="AliComparisonDCA", Char_t* title="AliComparisonDCA",Int_t analysisMode=0, Bool_t hptGenerator=kFALSE):
84 AliComparisonObject(name,title),
106 SetAnalysisMode(analysisMode);
107 SetHptGenerator(hptGenerator);
112 //_____________________________________________________________________________
113 AliComparisonDCA::~AliComparisonDCA()
116 if(fDCAHisto) delete fDCAHisto; fDCAHisto=0;
118 if(fD0TanSPtTPCITS) delete fD0TanSPtTPCITS; fD0TanSPtTPCITS=0;
119 if(fD1TanSPtTPCITS) delete fD1TanSPtTPCITS; fD1TanSPtTPCITS=0;
120 if(fD0TanSPt) delete fD0TanSPt; fD0TanSPt=0;
121 if(fD1TanSPt) delete fD1TanSPt; fD1TanSPt=0;
122 if(fD0TanSPtTPC) delete fD0TanSPtTPC; fD0TanSPtTPC=0;
123 if(fD1TanSPtTPC) delete fD1TanSPtTPC; fD1TanSPtTPC=0;
125 if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
129 //_____________________________________________________________________________
130 void AliComparisonDCA::Init()
135 Double_t binsP[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.};
136 Double_t pMin = 0., pMax = 10.;
138 if(IsHptGenerator() == kTRUE) {
140 pMin = 0.; pMax = 100.;
143 //dca_r, dca_z, eta, pt
144 Int_t binsQA[4] = {100,100,20,nPBins};
145 Double_t xminQA[4] = {-10.,-10.,-1., pMin};
146 Double_t xmaxQA[4] = {10.,10.,1., pMax};
148 fDCAHisto = new THnSparseF("fDCAHisto","dca_r:dca_z:eta:pt",4,binsQA,xminQA,xmaxQA);
149 if(!IsHptGenerator()) fDCAHisto->SetBinEdges(3,binsP);
151 fDCAHisto->GetAxis(0)->SetTitle("dca_r (cm)");
152 fDCAHisto->GetAxis(1)->SetTitle("dca_z (cm)");
153 fDCAHisto->GetAxis(2)->SetTitle("eta");
154 fDCAHisto->GetAxis(3)->SetTitle("pt (GeV/c)");
158 fD0TanSPtTPCITS = new TH3F("DCAyTanSPtTPCITS","DCAyTanSPt",40,-2,2, 10,0.3,3, 100,-1,1);
159 fD0TanSPtTPCITS->SetXTitle("tan(#theta)");
160 fD0TanSPtTPCITS->SetYTitle("#sqrt{p_{t}(GeV)}");
161 fD0TanSPtTPCITS->SetZTitle("DCA_{xy}");
163 fD1TanSPtTPCITS = new TH3F("DCAzTanSPtTPCITS","DCAzTanSPt",40,-2,2, 10,0.3,3, 100,-1,1);
164 fD1TanSPtTPCITS->SetXTitle("tan(#theta)");
165 fD1TanSPtTPCITS->SetYTitle("#sqrt(p_{t}(GeV))");
166 fD1TanSPtTPCITS->SetZTitle("DCA_{z}");
168 fD0TanSPt = new TH3F("DCAyTanSPt","DCAyTanSPt",40,-2,2, 10,0.3,3, 100,-1,1);
169 fD0TanSPt->SetXTitle("tan(#theta)");
170 fD0TanSPt->SetYTitle("#sqrt{p_{t}(GeV)}");
171 fD0TanSPt->SetZTitle("DCA_{xy}");
173 fD1TanSPt = new TH3F("DCAzTanSPt","DCAzTanSPt",40,-2,2, 10,0.3,3, 100, -1,1);
174 fD1TanSPt->SetXTitle("tan(#theta)");
175 fD1TanSPt->SetYTitle("#sqrt{p_{t}(GeV)}");
176 fD1TanSPt->SetZTitle("DCA_{z}");
178 fD0TanSPtTPC = new TH3F("DCAyTanSPtTPC","DCAyTanSPt",40,-2,2, 10,0.3,3, 100,-1,1);
179 fD0TanSPtTPC->SetXTitle("tan(#theta)");
180 fD0TanSPtTPC->SetYTitle("#sqrt{p_{t}(GeV)}");
181 fD0TanSPtTPC->SetZTitle("DCA_{xy}");
183 fD1TanSPtTPC = new TH3F("DCAzTanSPtTPC","DCAzTanSPt",40,-2,2, 10,0.3,3, 100, -1,1);
184 fD1TanSPtTPC->SetXTitle("tan(#theta)");
185 fD1TanSPtTPC->SetYTitle("#sqrt{p_{t}(GeV)}");
186 fD1TanSPtTPC->SetZTitle("DCA_{z}");
191 AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
193 AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
196 fAnalysisFolder = CreateFolder("folderDCA","Analysis DCA Folder");
199 //_____________________________________________________________________________
200 void AliComparisonDCA::ProcessTPC(AliMCInfo* const infoMC, AliESDRecInfo * const infoRC)
202 // Fill DCA comparison information
203 AliExternalTrackParam *track = 0;
204 //Double_t field = AliTracker::GetBz(); // nominal Bz field [kG]
205 Double_t kMaxD = 123456.0; // max distance
207 Double_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
209 Float_t mcpt = infoMC->GetParticle().Pt();
210 Float_t mceta = infoMC->GetParticle().Eta();
211 //Float_t spt = TMath::Sqrt(mcpt);
213 // distance to Prim. vertex
214 const Double_t* dv = infoMC->GetVDist();
216 Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();
218 // Check selection cuts
219 if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return;
221 if (infoRC->GetStatus(1)!=3) return;
222 if (!infoRC->GetESDtrack()) return;
223 if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
224 //if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;
226 // calculate and set prim. vertex
227 AliESDVertex vertexMC;
228 vertexMC.SetXv( infoMC->GetParticle().Vx() - dv[0] );
229 vertexMC.SetYv( infoMC->GetParticle().Vy() - dv[1] );
230 vertexMC.SetZv( infoMC->GetParticle().Vz() - dv[2] );
232 // calculate track parameters at vertex
233 if (infoRC->GetESDtrack()->GetTPCInnerParam())
235 if ((track = new AliExternalTrackParam(*infoRC->GetESDtrack()->GetTPCInnerParam())) != 0 )
238 //Bool_t bDCAStatus = track->PropagateToDCA(&vertexMC,field,kMaxD,dca,cov);
239 Double_t field[3]; track->GetBxByBz(field);
240 Bool_t bDCAStatus = track->PropagateToDCABxByBz(&vertexMC,field,kMaxD,dca,cov);
243 Double_t vDCAHisto[4]={dca[0],dca[1],mceta,mcpt};
244 fDCAHisto->Fill(vDCAHisto);
251 //_____________________________________________________________________________
252 void AliComparisonDCA::ProcessTPCITS(AliMCInfo* const infoMC, AliESDRecInfo * const infoRC)
254 // Fill DCA comparison information
255 Int_t clusterITS[200];
256 Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
258 Float_t mcpt = infoMC->GetParticle().Pt();
259 Float_t mceta = infoMC->GetParticle().Eta();
260 //Float_t spt = TMath::Sqrt(mcpt);
262 // distance to Prim. vertex
263 const Double_t* dv = infoMC->GetVDist();
264 Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();
266 // Check selection cuts
267 if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return;
269 if (infoRC->GetStatus(1)!=3) return;
270 if (!infoRC->GetESDtrack()) return;
271 if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
272 //if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;
274 infoRC->GetESDtrack()->GetImpactParameters(dca,cov);
277 if(infoRC->GetESDtrack()->GetITSclusters(clusterITS)>fCutsRC->GetMinNClustersITS())
279 Double_t vDCAHisto[4]={dca[0],dca[1],mceta,mcpt};
280 fDCAHisto->Fill(vDCAHisto);
284 void AliComparisonDCA::ProcessConstrained(AliMCInfo* const /*infoMC*/, AliESDRecInfo * const /*infoRC*/)
286 // Fill DCA comparison information
288 AliDebug(AliLog::kWarning, "Warning: Not implemented");
291 //_____________________________________________________________________________
292 Long64_t AliComparisonDCA::Merge(TCollection* const list)
294 // Merge list of objects (needed by PROOF)
302 TIterator* iter = list->MakeIterator();
305 // collection of generated histograms
307 while((obj = iter->Next()) != 0)
309 AliComparisonDCA* entry = dynamic_cast<AliComparisonDCA*>(obj);
310 if (entry == 0) continue;
312 fDCAHisto->Add(entry->fDCAHisto);
314 fD0TanSPtTPCITS->Add(entry->fD0TanSPtTPCITS);
315 fD1TanSPtTPCITS->Add(entry->fD1TanSPtTPCITS);
316 fD0TanSPt->Add(entry->fD0TanSPt);
317 fD1TanSPt->Add(entry->fD1TanSPt);
318 fD0TanSPtTPC->Add(entry->fD0TanSPtTPC);
319 fD1TanSPtTPC->Add(entry->fD1TanSPtTPC);
328 //_____________________________________________________________________________
329 void AliComparisonDCA::Exec(AliMCInfo* const infoMC, AliESDRecInfo * const infoRC)
331 // Process comparison information
332 if(GetAnalysisMode() == 0) ProcessTPC(infoMC,infoRC);
333 else if(GetAnalysisMode() == 1) ProcessTPCITS(infoMC,infoRC);
334 else if(GetAnalysisMode() == 2) ProcessConstrained(infoMC,infoRC);
336 printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
341 //_____________________________________________________________________________
342 void AliComparisonDCA::Analyse()
345 // Analyse comparison information and store output histograms
346 // in the analysis folder "folderDCA"
349 TH1::AddDirectory(kFALSE);
350 TObjArray *aFolderObj = new TObjArray;
353 TGraph * gr[8]= { 0,0,0,0,0,0,0,0 };
354 TGraph2D *gr2[8]= { 0,0,0,0,0,0,0,0};
355 AliComparisonDCA * comp=this;
357 // write results in the folder
358 // Canvas to draw analysed histograms
359 TCanvas * c = new TCanvas("canDCA","DCA resolution");
365 gr[0] = AliMathBase::MakeStat1D(comp->fD0TanSPtTPC,2,5);
366 gr[0]->GetXaxis()->SetTitle("Tan(#theta)");
367 gr[0]->GetYaxis()->SetTitle("#sigmaDCA_xy (cm)");
368 gr[0]->SetName("DCAXYResolTanTPC");
369 gr[0]->SetTitle("resol. DCA_xy (TPC only)");
372 aFolderObj->Add(gr[0]);
375 gr[1] = AliMathBase::MakeStat1D(comp->fD1TanSPtTPC,2,5);
376 gr[1]->GetXaxis()->SetTitle("Tan(#theta)");
377 gr[1]->GetYaxis()->SetTitle("#sigmaDCA_z (cm)");
378 gr[1]->SetName("DCAZResolTanTPC");
379 gr[1]->SetTitle("resol. DCA_z (TPC only)");
382 aFolderObj->Add(gr[1]);
385 gr[2] = AliMathBase::MakeStat1D(comp->fD0TanSPtTPCITS,2,5);
386 gr[2]->GetXaxis()->SetTitle("Tan(#theta)");
387 gr[2]->GetYaxis()->SetTitle("#sigmaDCA_xy (cm)");
388 gr[2]->SetName("DCAXYResolTanTPCITS");
389 gr[2]->SetTitle("resol. DCA_xy (TPC+ITS)");
392 aFolderObj->Add(gr[2]);
395 gr[3] = AliMathBase::MakeStat1D(comp->fD1TanSPtTPCITS,2,5);
396 gr[3]->GetXaxis()->SetTitle("Tan(#theta)");
397 gr[3]->GetYaxis()->SetTitle("#sigmaDCA_z (cm)");
398 gr[3]->SetName("DCAZResolTanTPCITS");
399 gr[3]->SetTitle("resol. DCA_z (TPC+ITS)");
402 aFolderObj->Add(gr[3]);
408 gr[4] = AliMathBase::MakeStat1D(comp->fD0TanSPtTPC,2,4);
409 gr[4]->GetXaxis()->SetTitle("Tan(#theta)");
410 gr[4]->GetYaxis()->SetTitle("mean DCA_xy (cm)");
411 gr[4]->SetName("DCAXYMeanTanTPC");
412 gr[4]->SetTitle("mean DCA_xy (TPC only)");
415 aFolderObj->Add(gr[4]);
418 gr[5] = AliMathBase::MakeStat1D(comp->fD1TanSPtTPC,2,4);
419 gr[5]->GetXaxis()->SetTitle("Tan(#theta)");
420 gr[5]->GetYaxis()->SetTitle("mean DCA_z (cm)");
421 gr[5]->SetName("DCAZMeanTanTPC");
422 gr[5]->SetTitle("mean DCA_z (TPC only)");
425 aFolderObj->Add(gr[5]);
428 gr[6] = AliMathBase::MakeStat1D(comp->fD0TanSPtTPCITS,2,4);
429 gr[6]->GetXaxis()->SetTitle("Tan(#theta)");
430 gr[6]->GetYaxis()->SetTitle("mean DCA_xy (cm)");
431 gr[6]->SetName("DCAXYMeanTanTPCITS");
432 gr[6]->SetTitle("mean DCA_xy (TPC+ITS)");
435 aFolderObj->Add(gr[6]);
438 gr[7] = AliMathBase::MakeStat1D(comp->fD1TanSPtTPCITS,2,4);
439 gr[7]->GetXaxis()->SetTitle("Tan(#theta)");
440 gr[7]->GetYaxis()->SetTitle("mean DCA_z (cm)");
441 gr[7]->SetName("DCAZMeanTanTPCITS");
442 gr[7]->SetTitle("mean DCA_z (TPC+ITS)");
445 aFolderObj->Add(gr[7]);
448 TCanvas * c1 = new TCanvas("canDCA1","2D DCA resolution");
453 gr2[0] = AliMathBase::MakeStat2D(comp->fD0TanSPtTPC,4,2,5);
454 gr2[0]->GetXaxis()->SetTitle("Tan(#theta)");
455 gr2[0]->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");
456 gr2[0]->GetZaxis()->SetTitle("#sigmaDCA_xy (cm)");
457 gr2[0]->SetName("DCAXYResolSPTTanTPC");
458 gr2[0]->SetTitle("#sigma DCA_xy (TPC only)");
459 gr2[0]->GetHistogram()->Draw("colz");
461 gr2[0]->GetHistogram()->SetName("DCAXYResolSPTTanTPC");
462 aFolderObj->Add(gr2[0]->GetHistogram());
465 gr2[1] = AliMathBase::MakeStat2D(comp->fD1TanSPtTPC,4,2,5);
466 gr2[1]->GetXaxis()->SetTitle("Tan(#theta)");
467 gr2[1]->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");
468 gr2[1]->GetZaxis()->SetTitle("#sigmaDCA_z (cm)");
469 gr2[1]->SetName("DCAZResolSPTTanTPC");
470 gr2[1]->SetTitle("#sigma DCA_z (TPC only)");
471 gr2[1]->GetHistogram()->Draw("colz");
473 gr2[1]->GetHistogram()->SetName("DCAZResolSPTTanTPC");
474 aFolderObj->Add(gr2[1]->GetHistogram());
478 gr2[2] = AliMathBase::MakeStat2D(comp->fD0TanSPtTPCITS,4,2,5);
479 gr2[2]->GetXaxis()->SetTitle("Tan(#theta)");
480 gr2[2]->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");
481 gr2[2]->GetZaxis()->SetTitle("#sigmaDCA_xy (cm)");
482 gr2[2]->SetName("DCAXYResolSPTTanTPCITS");
483 gr2[2]->SetTitle("#sigma DCA_xy (TPC+ITS)");
484 gr2[2]->GetHistogram()->Draw("colz");
486 gr2[2]->GetHistogram()->SetName("DCAXYResolSPTTanTPCITS");
487 aFolderObj->Add(gr2[2]->GetHistogram());
490 gr2[3] = AliMathBase::MakeStat2D(comp->fD1TanSPtTPCITS,4,2,5);
491 gr2[3]->GetXaxis()->SetTitle("Tan(#theta)");
492 gr2[3]->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");
493 gr2[3]->GetZaxis()->SetTitle("#sigmaDCA_z (cm)");
494 gr2[3]->SetName("DCAZResolSPTTanTPCITS");
495 gr2[3]->SetTitle("#sigma DCA_z (TPC+ITS)");
496 gr2[3]->GetHistogram()->Draw("colz");
498 gr2[3]->GetHistogram()->SetName("DCAZResolSPTTanTPCITS");
499 aFolderObj->Add(gr2[3]->GetHistogram());
503 gr2[4] = AliMathBase::MakeStat2D(comp->fD0TanSPtTPC,4,2,4);
504 gr2[4]->GetXaxis()->SetTitle("Tan(#theta)");
505 gr2[4]->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");
506 gr2[4]->GetZaxis()->SetTitle("mean DCA_xy (cm)");
507 gr2[4]->SetName("DCAXYMeanSPTTanTPC");
508 gr2[4]->SetTitle("mean DCA_xy (TPC only)");
509 gr2[4]->GetHistogram()->Draw("colz");
511 gr2[4]->GetHistogram()->SetName("DCAXYMeanSPTTanTPC");
512 aFolderObj->Add(gr2[4]->GetHistogram());
515 gr2[5] = AliMathBase::MakeStat2D(comp->fD1TanSPtTPC,4,2,4);
516 gr2[5]->GetXaxis()->SetTitle("Tan(#theta)");
517 gr2[5]->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");
518 gr2[5]->GetZaxis()->SetTitle("mean DCA_z (cm)");
519 gr2[5]->SetName("DCAZMeanSPTTanTPC");
520 gr2[5]->SetTitle("mean DCA_z (TPC only)");
521 gr2[5]->GetHistogram()->Draw("colz");
523 gr2[5]->GetHistogram()->SetName("DCAZMeanSPTTanTPC");
524 aFolderObj->Add(gr2[5]->GetHistogram());
527 gr2[6] = AliMathBase::MakeStat2D(comp->fD0TanSPtTPCITS,4,2,4);
528 gr2[6]->GetXaxis()->SetTitle("Tan(#theta)");
529 gr2[6]->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");
530 gr2[6]->GetZaxis()->SetTitle("mean DCA_xy (cm)");
531 gr2[6]->SetName("DCAXYMeanSPTTanTPCITS");
532 gr2[6]->SetTitle("mean DCA_xy (TPC+ITS)");
533 gr2[6]->GetHistogram()->Draw("colz");
535 gr2[6]->GetHistogram()->SetName("DCAXYMeanSPTTanTPCITS");
536 aFolderObj->Add(gr2[6]->GetHistogram());
539 gr2[7] = AliMathBase::MakeStat2D(comp->fD1TanSPtTPCITS,4,2,4);
540 gr2[7]->GetXaxis()->SetTitle("Tan(#theta)");
541 gr2[7]->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");
542 gr2[7]->GetZaxis()->SetTitle("mean DCA_z (cm)");
543 gr2[7]->SetName("DCAZMeanSPTTanTPCITS");
544 gr2[7]->SetTitle("mean DCA_z (TPC+ITS)");
545 gr2[7]->GetHistogram()->Draw("colz");
547 gr2[7]->GetHistogram()->SetName("DCAZMeanSPTTanTPCITS");
548 aFolderObj->Add(gr2[7]->GetHistogram());
551 // export objects to analysis folder
552 fAnalysisFolder = ExportToFolder(aFolderObj);
554 // delete only TObjArray
555 if(aFolderObj) delete aFolderObj;
558 //_____________________________________________________________________________
559 TFolder* AliComparisonDCA::ExportToFolder(TObjArray * array)
561 // recreate folder avery time and export objects to new one
563 AliComparisonDCA * comp=this;
564 TFolder *folder = comp->GetAnalysisFolder();
567 TFolder *newFolder = 0;
569 Int_t size = array->GetSize();
572 // get name and title from old folder
573 name = folder->GetName();
574 title = folder->GetTitle();
580 newFolder = CreateFolder(name.Data(),title.Data());
581 newFolder->SetOwner();
583 // add objects to folder
585 newFolder->Add(array->At(i));
594 //_____________________________________________________________________________
595 TFolder* AliComparisonDCA::CreateFolder(TString name,TString title) {
596 // create folder for analysed histograms
598 folder = new TFolder(name.Data(),title.Data());