]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/AliComparisonDEdx.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Trigger Lut generation analysis task and how to (Bogdan)
[u/mrichter/AliRoot.git] / PWG1 / AliComparisonDEdx.cxx
CommitLineData
6d1c79ca 1//------------------------------------------------------------------------------\r
2// Implementation of AliComparisonDEdx class. It keeps information from \r
3// comparison of reconstructed and MC particle tracks. In addtion, \r
4// it keeps selection cuts used during comparison. The comparison \r
5// information is stored in the ROOT histograms. Analysis of these \r
6// histograms can be done by using Analyse() class function. The result of \r
7// the analysis (histograms) are stored in the output picture_dedx.root file.\r
8// \r
9// Author: J.Otwinowski 04/02/2008 \r
10//------------------------------------------------------------------------------\r
11\r
12\r
13/*\r
14 //after running analysis, read the file, and get component\r
15 gSystem->Load("libPWG1.so");\r
16 TFile f("Output.root");\r
17 AliComparisonDEdx * comp = (AliComparisonDEdx*)f.Get("AliComparisonDEdx");\r
18\r
19 // analyse comparison data (output stored in pictures_dedx.root)\r
20 comp->Analyse();\r
21\r
22 // TPC track length parameterisation\r
23 TF1 fl("fl","((min(250./(abs(x+0.000001)),250)-90))",0,2); // length function\r
24 TF1 fl2("fl2","[0]/((min(250./(abs(x+0.000001)),250)-90))^[1]",0,2);\r
25 fl2.SetParameter(1,1);\r
26 fl2.SetParameter(0,1);\r
27*/\r
28\r
29#include <iostream>\r
30\r
31#include "TFile.h"\r
32#include "TCint.h"\r
33#include "TH3F.h"\r
34#include "TH2F.h"\r
35#include "TF1.h"\r
36#include "TProfile.h"\r
37#include "TProfile2D.h"\r
38#include "TGraph2D.h"\r
39#include "TCanvas.h"\r
40#include "TGraph.h"\r
41//\r
42#include "AliESDEvent.h"\r
43#include "AliESD.h"\r
44#include "AliESDfriend.h"\r
45#include "AliESDfriendTrack.h"\r
46#include "AliRecInfoCuts.h" \r
47#include "AliMCInfoCuts.h" \r
48#include "AliLog.h" \r
49//\r
50#include "AliMathBase.h"\r
51#include "AliTreeDraw.h"\r
52//#include "TStatToolkit.h"\r
53\r
54#include "AliMCInfo.h" \r
55#include "AliESDRecInfo.h" \r
56#include "AliComparisonDEdx.h" \r
57\r
58using namespace std;\r
59\r
60ClassImp(AliComparisonDEdx)\r
61\r
62//_____________________________________________________________________________\r
63AliComparisonDEdx::AliComparisonDEdx():\r
64 TNamed("AliComparisonDEdx","AliComparisonDEdx"),\r
65\r
66 // dEdx \r
67 fTPCSignalNormTan(0), \r
68 fTPCSignalNormSPhi(0),\r
69 fTPCSignalNormTPhi(0), \r
70 //\r
71 fTPCSignalNormTanSPhi(0),\r
72 fTPCSignalNormTanTPhi(0),\r
73 fTPCSignalNormTanSPt(0), \r
74 \r
75 // Cuts \r
76 fCutsRC(0), \r
77 fCutsMC(0),\r
78 fMCPtMin(0),\r
79 fMCAbsTanThetaMax(0),\r
80 fMCPdgCode(0)\r
81{\r
82 InitHisto();\r
83 InitCuts();\r
84}\r
85\r
86//_____________________________________________________________________________\r
87AliComparisonDEdx::~AliComparisonDEdx(){\r
88 \r
89 if(fTPCSignalNormTan) delete fTPCSignalNormTan; fTPCSignalNormTan=0; \r
90 if(fTPCSignalNormSPhi) delete fTPCSignalNormSPhi; fTPCSignalNormSPhi=0;\r
91 if(fTPCSignalNormTPhi) delete fTPCSignalNormTPhi; fTPCSignalNormTPhi=0;\r
92 //\r
93 if(fTPCSignalNormTanSPhi) delete fTPCSignalNormTanSPhi; fTPCSignalNormTanSPhi=0;\r
94 if(fTPCSignalNormTanTPhi) delete fTPCSignalNormTanTPhi; fTPCSignalNormTanTPhi=0;\r
95 if(fTPCSignalNormTanSPt) delete fTPCSignalNormTanSPt; fTPCSignalNormTanSPt=0;\r
96}\r
97\r
98//_____________________________________________________________________________\r
99void AliComparisonDEdx::InitHisto()\r
100{\r
101 // Init histograms\r
102 \r
103 // TPC dEdx\r
104 fTPCSignalNormTan = new TH2F("CdEdxTan","CdEdxTan",50, -2,2, 40,30,70); \r
105 fTPCSignalNormTan->SetXTitle("tan(#theta)");\r
106 fTPCSignalNormTan->SetYTitle("rec. dE/dx / calc. dE/dx");\r
107\r
108 fTPCSignalNormSPhi = new TH2F("CdEdxSPhi","CdEdxSPhi",10,0.0,1,40,30,70);\r
109 fTPCSignalNormSPhi->SetXTitle("sin(#phi)");\r
110 fTPCSignalNormSPhi->SetYTitle("rec. dE/dx / calc. dE/dx");\r
111\r
112 fTPCSignalNormTPhi = new TH2F("CdEdxTPhi","CdEdxTPhi",10,0.0,2,40,30,70); \r
113 fTPCSignalNormTPhi->SetXTitle("tan(#phi)");\r
114 fTPCSignalNormTPhi->SetYTitle("rec. dE/dx / calc. dE/dx");\r
115\r
116 fTPCSignalNormTanSPhi= new TH3F("CdEdxTanSPhi","CdEdxTanSPhi",20, -2,2, 10,0.0 ,1, 40,30,70);\r
117 fTPCSignalNormTanSPhi->SetXTitle("tan(#theta)");\r
118 fTPCSignalNormTanSPhi->SetYTitle("sin(#phi)");\r
119 fTPCSignalNormTanSPhi->SetZTitle("rec. dE/dx / calc. dE/dx");\r
120\r
121 fTPCSignalNormTanTPhi= new TH3F("CdEdxTanTPhi","CdEdxTanTPhi",20, -2,2, 10,0.0 ,1, 40,30,70);\r
122 fTPCSignalNormTanTPhi->SetXTitle("tan(#theta)");\r
123 fTPCSignalNormTanTPhi->SetYTitle("tan(#phi)");\r
124 fTPCSignalNormTanTPhi->SetZTitle("rec. dE/dx / calc. dE/dx");\r
125\r
126 fTPCSignalNormTanSPt= new TH3F("CdEdxTanSPt","CdEdxTanSPt",20, -2,2, 10,0.3 ,3, 40,30,70); \r
127 fTPCSignalNormTanSPt->SetXTitle("tan(#theta)");\r
128 fTPCSignalNormTanSPt->SetYTitle("#sqrt{p_{t}}");\r
129 fTPCSignalNormTanSPt->SetZTitle("rec. dE/dx / calc. dE/dx");\r
130}\r
131\r
132void AliComparisonDEdx::InitCuts()\r
133{\r
134 // Init cuts\r
135 if(!fCutsMC) \r
136 AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");\r
137 if(!fCutsRC) \r
138 AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");\r
139}\r
140\r
141//_____________________________________________________________________________\r
142void AliComparisonDEdx::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC){\r
143\r
144 // Fill dE/dx comparison information\r
145 \r
146 Float_t mcpt = infoMC->GetParticle().Pt();\r
147 Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);\r
148 Float_t mprim = infoMC->GetPrim();\r
149\r
150 // distance to Prim. vertex \r
151 const Double_t* dv = infoMC->GetVDist(); \r
152\r
153 Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();\r
154 \r
155 // Check selection cuts \r
156 if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return; \r
157 if (!isPrim) return;\r
158 if (infoRC->GetStatus(1)!=3) return;\r
159 if (!infoRC->GetESDtrack()) return; \r
160 if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;\r
161 if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;\r
162 //if (mprim>1.4) return;\r
163 //if (mprim<0.5) return;\r
164 if (mprim > fCutsMC->GetMaxTPCSignal()) return;\r
165 if (mprim < fCutsMC->GetMinTPCSignal()) return;\r
166 if (infoRC->GetESDtrack()->GetTPCsignalN()<fCutsRC->GetMinTPCsignalN()) return;\r
167 //\r
168 Float_t ratio = infoRC->GetESDtrack()->GetTPCsignal()/infoMC->GetPrim();\r
169 Float_t sphi = infoRC->GetESDtrack()->GetInnerParam()->GetSnp();\r
170 Float_t tphi = sphi/TMath::Sqrt(1-sphi*sphi);\r
171\r
172 if (TMath::Abs(infoMC->GetParticle().GetPdgCode()) != GetMCPdgCode()) return;\r
173 //if (mcpt>0.5) {\r
174 if (mcpt > GetMCPtMin()) {\r
175 fTPCSignalNormTan->Fill(tantheta,ratio); // only subset\r
176 }\r
177\r
178 //if (TMath::Abs(tantheta)<0.5){\r
179 if (TMath::Abs(tantheta) < GetMCAbsTanThetaMax()){\r
180 fTPCSignalNormSPhi->Fill(sphi,ratio); // only subset\r
181 fTPCSignalNormTPhi->Fill(tphi,ratio); // only subset\r
182 }\r
183 fTPCSignalNormTanSPhi->Fill(tantheta,sphi,ratio); \r
184 fTPCSignalNormTanTPhi->Fill(tantheta,tphi,ratio); \r
185 fTPCSignalNormTanSPt->Fill(tantheta,TMath::Sqrt(mcpt),ratio);\r
186}\r
187\r
188//_____________________________________________________________________________\r
189Long64_t AliComparisonDEdx::Merge(TCollection* list) \r
190{\r
191 // Merge list of objects (needed by PROOF)\r
192\r
193 if (!list)\r
194 return 0;\r
195\r
196 if (list->IsEmpty())\r
197 return 1;\r
198\r
199 TIterator* iter = list->MakeIterator();\r
200 TObject* obj = 0;\r
201\r
202 // collection of generated histograms\r
203 Int_t count=0;\r
204 while((obj = iter->Next()) != 0) \r
205 {\r
206 AliComparisonDEdx* entry = dynamic_cast<AliComparisonDEdx*>(obj);\r
207 if (entry == 0) { \r
208 Error("Add","Attempt to add object of class: %s to a %s",\r
209 obj->ClassName(),this->ClassName());\r
210 return -1;\r
211 }\r
212\r
213 fTPCSignalNormTan->Add(entry->fTPCSignalNormTan);\r
214 fTPCSignalNormSPhi->Add(entry->fTPCSignalNormSPhi);\r
215 fTPCSignalNormTPhi->Add(entry->fTPCSignalNormTPhi);\r
216 //\r
217 fTPCSignalNormTanSPhi->Add(entry->fTPCSignalNormTanSPhi);\r
218 fTPCSignalNormTanTPhi->Add(entry->fTPCSignalNormTanTPhi);\r
219 fTPCSignalNormTanSPt->Add(entry->fTPCSignalNormTanSPt);\r
220\r
221 count++;\r
222 }\r
223\r
224return count;\r
225}\r
226\r
227//_____________________________________________________________________________\r
228void AliComparisonDEdx::Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC)\r
229{\r
230 // Process comparison information\r
231 Process(infoMC,infoRC);\r
232}\r
233\r
234//_____________________________________________________________________________\r
235TH1F* AliComparisonDEdx::MakeResol(TH2F * his, Int_t integ, Bool_t type)\r
236{\r
237 // Make resolution histograms\r
238 TH1F *hisr, *hism;\r
239 if (!gPad) new TCanvas;\r
240 hisr = AliTreeDraw::CreateResHistoI(his,&hism,integ);\r
241 if (type) return hism;\r
242 else \r
243 return hisr;\r
244}\r
245\r
246//_____________________________________________________________________________\r
247void AliComparisonDEdx::Analyse()\r
248{\r
249 // Analyse output histograms\r
250 \r
251 AliComparisonDEdx * comp=this;\r
252 TH1F *hiss=0;\r
253 TGraph2D * gr=0;\r
254\r
255 TFile *fp = new TFile("pictures_dedx.root","recreate");\r
256 fp->cd();\r
257\r
258 TCanvas * c = new TCanvas("TPCdedx","TPC dedx");\r
259 c->cd();\r
260\r
261 hiss = comp->MakeResol(comp->fTPCSignalNormTan,4,0);\r
262 hiss->SetXTitle("Tan(#theta)");\r
263 hiss->SetYTitle("#sigma_{dEdx}");\r
264 hiss->Draw();\r
265 hiss->Write("TPCdEdxResolTan");\r
266 //\r
267 hiss = comp->MakeResol(comp->fTPCSignalNormTan,4,1); \r
268 hiss->SetXTitle("Tan(#theta)");\r
269 hiss->SetYTitle("<dEdx>");\r
270 hiss->Draw(); \r
271 hiss->Write("TPCdEdxMeanTan");\r
272 //\r
273 gr = AliMathBase::MakeStat2D(comp->fTPCSignalNormTanSPt,3,1,4);\r
274 gr->GetXaxis()->SetTitle("Tan(#theta)");\r
275 gr->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");\r
276 gr->GetZaxis()->SetTitle("<dEdx>");\r
277 gr->Draw("colz"); \r
278 gr->GetHistogram()->Write("TPCdEdxMeanTanPt_1");\r
279 //\r
280 gr = AliMathBase::MakeStat2D(comp->fTPCSignalNormTanSPt,3,1,5);\r
281 gr->GetXaxis()->SetTitle("Tan(#theta)");\r
282 gr->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");\r
283 gr->GetZaxis()->SetTitle("#sigma_{dEdx}");\r
284 gr->Draw("colz"); \r
285 gr->GetHistogram()->Write("TPCdEdxMeanTanPt_2");\r
286\r
287 fp->Close();\r
288}\r