1 //------------------------------------------------------------------------------
2 // Implementation of AliPerformanceDEdx 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 which is
8 // a data of AliPerformanceDEdx.
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 //AliPerformanceDEdx * compObj = (AliPerformanceDEdx*)f.Get("AliPerformanceDEdx");
20 AliPerformanceDEdx * compObj = (AliPerformanceDEdx*)coutput->FindObject("AliPerformanceDEdx");
22 // Analyse comparison data
25 // the output histograms/graphs will be stored in the folder "folderDEdx"
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_DEdx.root"."recreate");
31 compObj->Write(); // compObj->GetAnalysisFolder()->Write();
36 #include <TDirectory.h>
43 #include "AliPerformanceDEdx.h"
44 #include "AliESDEvent.h"
45 #include "AliMCEvent.h"
46 #include "AliESDtrack.h"
49 #include "AliMCInfoCuts.h"
50 #include "AliMathBase.h"
51 #include "AliRecInfoCuts.h"
52 #include "AliTreeDraw.h"
53 #include "AliHeader.h"
54 #include "AliGenEventHeader.h"
58 ClassImp(AliPerformanceDEdx)
60 //_____________________________________________________________________________
61 AliPerformanceDEdx::AliPerformanceDEdx():
62 AliPerformanceObject("AliPerformanceDEdx"),
74 // default constructor
78 //_____________________________________________________________________________
79 AliPerformanceDEdx::AliPerformanceDEdx(Char_t* name="AliPerformanceDEdx", Char_t* title="AliPerformanceDEdx",Int_t analysisMode=0, Bool_t hptGenerator=kFALSE):
80 AliPerformanceObject(name,title),
94 SetAnalysisMode(analysisMode);
95 SetHptGenerator(hptGenerator);
100 //_____________________________________________________________________________
101 AliPerformanceDEdx::~AliPerformanceDEdx()
104 if(fDeDxHisto) delete fDeDxHisto; fDeDxHisto=0;
105 if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
108 //_____________________________________________________________________________
109 void AliPerformanceDEdx::Init()
116 Double_t pMin = 1.e-2, pMax = 10.;
119 if (IsHptGenerator()) {
120 nPBins = 100; pMax = 100.;
121 binsP = CreateLogAxis(nPBins,pMin,pMax);
123 binsP = CreateLogAxis(nPBins,pMin,pMax);
129 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.};
130 Double_t pMin = 0., pMax = 10.;
132 if(IsHptGenerator() == kTRUE) {
134 pMin = 0.; pMax = 100.;
138 //dedx:alpha:y:z:snp:tgl:ncls:p
139 //dedx:phi:y:z:snp:tgl:ncls:p
140 //Int_t binsQA[8] = {300, 50, 50, 50, 50, 50, 80, nPBins};
141 //Double_t xminQA[8] = {0, -4,-20,-250, -1, -2, 0, pMin};
142 //Double_t xmaxQA[8] = {300, 4, 20, 250, 1, 2, 160, pMax};
143 Int_t binsQA[8] = {300, 144, 50, 50, 50, 50, 80, nPBins};
144 Double_t xminQA[8] = {0, -TMath::Pi(),-20,-250, -1, -2, 0, pMin};
145 Double_t xmaxQA[8] = {300, TMath::Pi(), 20, 250, 1, 2, 160, pMax};
147 //fDeDxHisto = new THnSparseF("fDeDxHisto","dedx:alpha:y:z:snp:tgl:ncls:momentum",8,binsQA,xminQA,xmaxQA);
148 fDeDxHisto = new THnSparseF("fDeDxHisto","dedx:phi:y:z:snp:tgl:ncls:momentum",8,binsQA,xminQA,xmaxQA);
149 fDeDxHisto->SetBinEdges(7,binsP);
151 fDeDxHisto->GetAxis(0)->SetTitle("dedx (a.u.)");
152 fDeDxHisto->GetAxis(1)->SetTitle("#phi (rad)");
153 fDeDxHisto->GetAxis(2)->SetTitle("y (cm)");
154 fDeDxHisto->GetAxis(3)->SetTitle("z (cm)");
155 fDeDxHisto->GetAxis(4)->SetTitle("sin#phi");
156 fDeDxHisto->GetAxis(5)->SetTitle("tan#lambda");
157 fDeDxHisto->GetAxis(6)->SetTitle("ncls");
158 fDeDxHisto->GetAxis(7)->SetTitle("p (GeV/c)");
163 AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
165 AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
168 fAnalysisFolder = CreateFolder("folderDEdx","Analysis de/dx Folder");
171 //_____________________________________________________________________________
172 void AliPerformanceDEdx::ProcessTPC(AliStack* const /*stack*/, AliESDtrack *const /*esdTrack*/)
174 // Fill dE/dx comparison information
175 AliDebug(AliLog::kWarning, "Warning: Not implemented");
178 //_____________________________________________________________________________
179 void AliPerformanceDEdx::ProcessInnerTPC(AliStack* const stack, AliESDtrack *const esdTrack)
181 if(!esdTrack) return;
183 const AliExternalTrackParam *innerParam = esdTrack->GetInnerParam();
184 if(!innerParam) return;
186 Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
187 esdTrack->GetImpactParametersTPC(dca,cov);
189 if((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
194 Double_t dcaToVertex = -1;
195 if( fCutsRC->GetDCAToVertex2D() )
197 dcaToVertex = TMath::Sqrt(dca[0]*dca[0]/fCutsRC->GetMaxDCAToVertexXY()/fCutsRC->GetMaxDCAToVertexXY() + dca[1]*dca[1]/fCutsRC->GetMaxDCAToVertexZ()/fCutsRC->GetMaxDCAToVertexZ());
199 if(fCutsRC->GetDCAToVertex2D() && dcaToVertex > 1) return;
200 if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[0]) > fCutsRC->GetMaxDCAToVertexXY()) return;
201 if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[1]) > fCutsRC->GetMaxDCAToVertexZ()) return;
203 Float_t dedx = esdTrack->GetTPCsignal();
204 Int_t ncls = esdTrack->GetTPCNcls();
206 Double_t pt = innerParam->Pt();
207 Double_t lam = TMath::ATan2(innerParam->Pz(),innerParam->Pt());
208 Double_t p = pt/TMath::Cos(lam);
209 //Double_t alpha = innerParam->GetAlpha();
210 Double_t phi = TMath::ATan2(innerParam->Py(),innerParam->Px());
211 //if(phi<0.) phi += 2.*TMath::Phi();
212 Double_t y = innerParam->GetY();
213 Double_t z = innerParam->GetZ();
214 Double_t snp = innerParam->GetSnp();
215 Double_t tgl = innerParam->GetTgl();
217 //Double_t vDeDxHisto[8] = {dedx,alpha,y,z,snp,tgl,ncls,p};
218 Double_t vDeDxHisto[8] = {dedx,phi,y,z,snp,tgl,ncls,p};
219 fDeDxHisto->Fill(vDeDxHisto);
224 //_____________________________________________________________________________
225 void AliPerformanceDEdx::ProcessTPCITS(AliStack* const /*stack*/, AliESDtrack *const /*esdTrack*/)
227 // Fill dE/dx comparison information
229 AliDebug(AliLog::kWarning, "Warning: Not implemented");
232 //_____________________________________________________________________________
233 void AliPerformanceDEdx::ProcessConstrained(AliStack* const /*stack*/, AliESDtrack *const /*esdTrack*/)
235 // Fill dE/dx comparison information
237 AliDebug(AliLog::kWarning, "Warning: Not implemented");
240 //_____________________________________________________________________________
241 Long64_t AliPerformanceDEdx::Merge(TCollection* const list)
243 // Merge list of objects (needed by PROOF)
251 TIterator* iter = list->MakeIterator();
254 // collection of generated histograms
256 while((obj = iter->Next()) != 0)
258 AliPerformanceDEdx* entry = dynamic_cast<AliPerformanceDEdx*>(obj);
259 if (entry == 0) continue;
261 fDeDxHisto->Add(entry->fDeDxHisto);
268 //_____________________________________________________________________________
269 void AliPerformanceDEdx::Exec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent, AliESDfriend *const esdFriend, const Bool_t bUseMC, const Bool_t bUseESDfriend)
271 // Process comparison information
275 AliDebug(AliLog::kError, "esdEvent not available");
278 AliHeader* header = 0;
279 AliGenEventHeader* genHeader = 0;
286 AliDebug(AliLog::kError, "mcEvent not available");
290 // get MC event header
291 header = mcEvent->Header();
293 AliDebug(AliLog::kError, "Header not available");
297 stack = mcEvent->Stack();
299 AliDebug(AliLog::kError, "Stack not available");
304 genHeader = header->GenEventHeader();
306 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
309 genHeader->PrimaryVertex(vtxMC);
316 AliDebug(AliLog::kError, "esdFriend not available");
322 Bool_t isEventTriggered = esdEvent->IsTriggerClassFired(GetTriggerClass());
323 if(!isEventTriggered) return;
325 // get TPC event vertex
326 const AliESDVertex *vtxESD = esdEvent->GetPrimaryVertexTPC();
327 if(vtxESD && (vtxESD->GetStatus()<=0)) return;
330 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
332 AliESDtrack *track = esdEvent->GetTrack(iTrack);
335 if(GetAnalysisMode() == 0) ProcessTPC(stack,track);
336 else if(GetAnalysisMode() == 1) ProcessTPCITS(stack,track);
337 else if(GetAnalysisMode() == 2) ProcessConstrained(stack,track);
338 else if(GetAnalysisMode() == 3) ProcessInnerTPC(stack,track);
340 printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
346 //_____________________________________________________________________________
347 void AliPerformanceDEdx::Analyse()
349 // Analyze comparison information and store output histograms
350 // in the folder "folderDEdx"
352 TH1::AddDirectory(kFALSE);
355 TObjArray *aFolderObj = new TObjArray;
360 for(Int_t i=1; i<8; i++) {
362 h2D = (TH2F*)fDeDxHisto->Projection(0,i);
364 sprintf(name,"h_dedx_%d_vs_%d",0,i);
366 sprintf(title,"%s vs %s",fDeDxHisto->GetAxis(0)->GetTitle(),fDeDxHisto->GetAxis(i)->GetTitle());
367 h2D->SetTitle(title);
368 h2D->GetXaxis()->SetTitle(fDeDxHisto->GetAxis(i)->GetTitle());
369 h2D->GetYaxis()->SetTitle(fDeDxHisto->GetAxis(0)->GetTitle());
371 if(i==7) h2D->SetBit(TH1::kLogX);
372 aFolderObj->Add(h2D);
375 // resolution histograms for mips
376 //dedx:phi:y:z:snp:tgl:ncls:p
377 fDeDxHisto->GetAxis(2)->SetRangeUser(-15.,14.999);
378 fDeDxHisto->GetAxis(3)->SetRangeUser(-120.,119.999);
379 fDeDxHisto->GetAxis(4)->SetRangeUser(-0.4, 0.399);
380 fDeDxHisto->GetAxis(5)->SetRangeUser(-0.9,0.89);
381 fDeDxHisto->GetAxis(6)->SetRangeUser(60.,160.);
382 fDeDxHisto->GetAxis(7)->SetRangeUser(0.4,0.499);
384 h1D=(TH1F*)fDeDxHisto->Projection(0);
385 h1D->SetName("dedx_mips");
387 h1D->SetTitle("dedx_mips");
388 h1D->GetXaxis()->SetTitle(fDeDxHisto->GetAxis(0)->GetTitle());
389 aFolderObj->Add(h1D);
392 TObjArray *arr[7] = {0};
395 for(Int_t i=1; i<8; i++)
397 arr[i] = new TObjArray;
398 f1[i] = new TF1("gaus","gaus");
399 //printf("i %d \n",i);
401 h2D = (TH2F*)fDeDxHisto->Projection(0,i);
403 f1[i]->SetRange(40,60); // should be pion peak
404 h2D->FitSlicesY(f1[i],0,-1,10,"QNR",arr[i]); // gaus fit of pion peak
406 h1D = (TH1F*)arr[i]->At(1);
407 sprintf(name,"mean_dedx_mips_vs_%d",i);
409 sprintf(title,"%s vs %s","mean_dedx_mips (a.u.)",fDeDxHisto->GetAxis(i)->GetTitle());
410 h1D->SetTitle(title);
411 h1D->GetXaxis()->SetTitle(fDeDxHisto->GetAxis(i)->GetTitle());
412 h1D->GetYaxis()->SetTitle("mean_dedx_mips (a.u.)");
413 //h1D->SetMinimum(40);
414 //h1D->SetMaximum(60);
416 aFolderObj->Add(h1D);
418 h1D = (TH1F*)arr[i]->At(2);
419 sprintf(name,"res_dedx_mips_vs_%d",i);
421 sprintf(title,"%s vs %s","res_dedx_mips (a.u)",fDeDxHisto->GetAxis(i)->GetTitle());
422 h1D->SetTitle(title);
423 h1D->GetXaxis()->SetTitle(fDeDxHisto->GetAxis(i)->GetTitle());
424 h1D->GetYaxis()->SetTitle("res_dedx_mips (a.u.)");
425 //h1D->SetMinimum(0);
426 //h1D->SetMaximum(6);
428 aFolderObj->Add(h1D);
431 // export objects to analysis folder
432 fAnalysisFolder = ExportToFolder(aFolderObj);
434 // delete only TObjArray
435 if(aFolderObj) delete aFolderObj;
438 //_____________________________________________________________________________
439 TFolder* AliPerformanceDEdx::ExportToFolder(TObjArray * array)
441 // recreate folder avery time and export objects to new one
443 AliPerformanceDEdx * comp=this;
444 TFolder *folder = comp->GetAnalysisFolder();
447 TFolder *newFolder = 0;
449 Int_t size = array->GetSize();
452 // get name and title from old folder
453 name = folder->GetName();
454 title = folder->GetTitle();
460 newFolder = CreateFolder(name.Data(),title.Data());
461 newFolder->SetOwner();
463 // add objects to folder
465 newFolder->Add(array->At(i));
474 //_____________________________________________________________________________
475 TFolder* AliPerformanceDEdx::CreateFolder(TString name,TString title) {
476 // create folder for analysed histograms
478 folder = new TFolder(name.Data(),title.Data());