1 //------------------------------------------------------------------------------
2 // Implementation of AliComparisonDEdx 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) are stored in the output picture_dedx.root file.
9 // Author: J.Otwinowski 04/02/2008
10 //------------------------------------------------------------------------------
14 //after running analysis, read the file, and get component
15 gSystem->Load("libPWG1.so");
16 TFile f("Output.root");
17 AliComparisonDEdx * comp = (AliComparisonDEdx*)f.Get("AliComparisonDEdx");
19 // analyse comparison data (output stored in pictures_dedx.root)
22 // TPC track length parameterisation
23 TF1 fl("fl","((min(250./(abs(x+0.000001)),250)-90))",0,2); // length function
24 TF1 fl2("fl2","[0]/((min(250./(abs(x+0.000001)),250)-90))^[1]",0,2);
25 fl2.SetParameter(1,1);
26 fl2.SetParameter(0,1);
37 #include "TProfile2D.h"
42 #include "AliESDEvent.h"
44 #include "AliESDfriend.h"
45 #include "AliESDfriendTrack.h"
46 #include "AliRecInfoCuts.h"
47 #include "AliMCInfoCuts.h"
50 #include "AliMathBase.h"
51 #include "AliTreeDraw.h"
52 //#include "TStatToolkit.h"
54 #include "AliMCInfo.h"
55 #include "AliESDRecInfo.h"
56 #include "AliComparisonDEdx.h"
60 ClassImp(AliComparisonDEdx)
62 //_____________________________________________________________________________
63 AliComparisonDEdx::AliComparisonDEdx():
64 TNamed("AliComparisonDEdx","AliComparisonDEdx"),
68 fTPCSignalNormSPhi(0),
69 fTPCSignalNormTPhi(0),
71 fTPCSignalNormTanSPhi(0),
72 fTPCSignalNormTanTPhi(0),
73 fTPCSignalNormTanSPt(0),
86 //_____________________________________________________________________________
87 AliComparisonDEdx::~AliComparisonDEdx(){
89 if(fTPCSignalNormTan) delete fTPCSignalNormTan; fTPCSignalNormTan=0;
90 if(fTPCSignalNormSPhi) delete fTPCSignalNormSPhi; fTPCSignalNormSPhi=0;
91 if(fTPCSignalNormTPhi) delete fTPCSignalNormTPhi; fTPCSignalNormTPhi=0;
93 if(fTPCSignalNormTanSPhi) delete fTPCSignalNormTanSPhi; fTPCSignalNormTanSPhi=0;
94 if(fTPCSignalNormTanTPhi) delete fTPCSignalNormTanTPhi; fTPCSignalNormTanTPhi=0;
95 if(fTPCSignalNormTanSPt) delete fTPCSignalNormTanSPt; fTPCSignalNormTanSPt=0;
98 //_____________________________________________________________________________
99 void AliComparisonDEdx::InitHisto()
104 fTPCSignalNormTan = new TH2F("CdEdxTan","CdEdxTan",50, -2,2, 40,30,70);
105 fTPCSignalNormTan->SetXTitle("tan(#theta)");
106 fTPCSignalNormTan->SetYTitle("rec. dE/dx / calc. dE/dx");
108 fTPCSignalNormSPhi = new TH2F("CdEdxSPhi","CdEdxSPhi",10,0.0,1,40,30,70);
109 fTPCSignalNormSPhi->SetXTitle("sin(#phi)");
110 fTPCSignalNormSPhi->SetYTitle("rec. dE/dx / calc. dE/dx");
112 fTPCSignalNormTPhi = new TH2F("CdEdxTPhi","CdEdxTPhi",10,0.0,2,40,30,70);
113 fTPCSignalNormTPhi->SetXTitle("tan(#phi)");
114 fTPCSignalNormTPhi->SetYTitle("rec. dE/dx / calc. dE/dx");
116 fTPCSignalNormTanSPhi= new TH3F("CdEdxTanSPhi","CdEdxTanSPhi",20, -2,2, 10,0.0 ,1, 40,30,70);
117 fTPCSignalNormTanSPhi->SetXTitle("tan(#theta)");
118 fTPCSignalNormTanSPhi->SetYTitle("sin(#phi)");
119 fTPCSignalNormTanSPhi->SetZTitle("rec. dE/dx / calc. dE/dx");
121 fTPCSignalNormTanTPhi= new TH3F("CdEdxTanTPhi","CdEdxTanTPhi",20, -2,2, 10,0.0 ,1, 40,30,70);
122 fTPCSignalNormTanTPhi->SetXTitle("tan(#theta)");
123 fTPCSignalNormTanTPhi->SetYTitle("tan(#phi)");
124 fTPCSignalNormTanTPhi->SetZTitle("rec. dE/dx / calc. dE/dx");
126 fTPCSignalNormTanSPt= new TH3F("CdEdxTanSPt","CdEdxTanSPt",20, -2,2, 10,0.3 ,3, 40,30,70);
127 fTPCSignalNormTanSPt->SetXTitle("tan(#theta)");
128 fTPCSignalNormTanSPt->SetYTitle("#sqrt(p_{t})");
129 fTPCSignalNormTanSPt->SetZTitle("rec. dE/dx / calc. dE/dx");
132 void AliComparisonDEdx::InitCuts()
136 AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
138 AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
141 //_____________________________________________________________________________
142 void AliComparisonDEdx::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
144 // Fill dE/dx comparison information
146 Float_t mcpt = infoMC->GetParticle().Pt();
147 Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
148 Bool_t isPrim = infoMC->GetParticle().R()<fCutsMC->GetMaxR() && TMath::Abs(infoMC->GetParticle().Vz())<fCutsMC->GetMaxVz() ;
150 Float_t mprim = infoMC->GetPrim();
152 // Check selection cuts
153 if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return;
155 if (infoRC->GetStatus(1)!=3) return;
156 if (!infoRC->GetESDtrack()) return;
157 if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
158 if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;
159 //if (mprim>1.4) return;
160 //if (mprim<0.5) return;
161 if (mprim > fCutsMC->GetMaxTPCSignal()) return;
162 if (mprim < fCutsMC->GetMinTPCSignal()) return;
163 if (infoRC->GetESDtrack()->GetTPCsignalN()<fCutsRC->GetMinTPCsignalN()) return;
165 Float_t ratio = infoRC->GetESDtrack()->GetTPCsignal()/infoMC->GetPrim();
166 Float_t sphi = infoRC->GetESDtrack()->GetInnerParam()->GetSnp();
167 Float_t tphi = sphi/TMath::Sqrt(1-sphi*sphi);
169 if (TMath::Abs(infoMC->GetParticle().GetPdgCode()) != GetMCPdgCode()) return;
171 if (mcpt > GetMCPtMin()) {
172 fTPCSignalNormTan->Fill(tantheta,ratio); // only subset
175 //if (TMath::Abs(tantheta)<0.5){
176 if (TMath::Abs(tantheta) < GetMCAbsTanThetaMax()){
177 fTPCSignalNormSPhi->Fill(sphi,ratio); // only subset
178 fTPCSignalNormTPhi->Fill(tphi,ratio); // only subset
180 fTPCSignalNormTanSPhi->Fill(tantheta,sphi,ratio);
181 fTPCSignalNormTanTPhi->Fill(tantheta,tphi,ratio);
182 fTPCSignalNormTanSPt->Fill(tantheta,TMath::Sqrt(mcpt),ratio);
185 //_____________________________________________________________________________
186 Long64_t AliComparisonDEdx::Merge(TCollection* list)
188 // Merge list of objects (needed by PROOF)
196 TIterator* iter = list->MakeIterator();
199 // collection of generated histograms
201 while((obj = iter->Next()) != 0)
203 AliComparisonDEdx* entry = dynamic_cast<AliComparisonDEdx*>(obj);
205 Error("Add","Attempt to add object of class: %s to a %s",
206 obj->ClassName(),this->ClassName());
210 fTPCSignalNormTan->Add(entry->fTPCSignalNormTan);
211 fTPCSignalNormSPhi->Add(entry->fTPCSignalNormSPhi);
212 fTPCSignalNormTPhi->Add(entry->fTPCSignalNormTPhi);
214 fTPCSignalNormTanSPhi->Add(entry->fTPCSignalNormTanSPhi);
215 fTPCSignalNormTanTPhi->Add(entry->fTPCSignalNormTanTPhi);
216 fTPCSignalNormTanSPt->Add(entry->fTPCSignalNormTanSPt);
224 //_____________________________________________________________________________
225 void AliComparisonDEdx::Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
227 // Process comparison information
228 Process(infoMC,infoRC);
231 //_____________________________________________________________________________
232 TH1F* AliComparisonDEdx::MakeResol(TH2F * his, Int_t integ, Bool_t type)
234 // Make resolution histograms
236 if (!gPad) new TCanvas;
237 hisr = AliTreeDraw::CreateResHistoI(his,&hism,integ);
238 if (type) return hism;
243 //_____________________________________________________________________________
244 void AliComparisonDEdx::Analyse()
246 // Analyse output histograms
248 AliComparisonDEdx * comp=this;
252 TFile *fp = new TFile("pictures_dedx.root","recreate");
255 TCanvas * c = new TCanvas("TPCdedx","TPC dedx");
258 hiss = comp->MakeResol(comp->fTPCSignalNormTan,4,0);
259 hiss->SetXTitle("Tan(#theta)");
260 hiss->SetYTitle("#sigma_{dEdx}");
262 hiss->Write("TPCdEdxResolTan");
264 hiss = comp->MakeResol(comp->fTPCSignalNormTan,4,1);
265 hiss->SetXTitle("Tan(#theta)");
266 hiss->SetYTitle("<dEdx>");
268 hiss->Write("TPCdEdxMeanTan");
270 gr = AliMathBase::MakeStat2D(comp->fTPCSignalNormTanSPt,3,1,4);
271 gr->GetXaxis()->SetTitle("Tan(#theta)");
272 gr->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");
273 gr->GetZaxis()->SetTitle("<dEdx>");
275 gr->GetHistogram()->Write("TPCdEdxMeanTanPt_1");
277 gr = AliMathBase::MakeStat2D(comp->fTPCSignalNormTanSPt,3,1,5);
278 gr->GetXaxis()->SetTitle("Tan(#theta)");
279 gr->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");
280 gr->GetZaxis()->SetTitle("#sigma_{dEdx}");
282 gr->GetHistogram()->Write("TPCdEdxMeanTanPt_2");