]>
Commit | Line | Data |
---|---|---|
7cc34f08 | 1 | //------------------------------------------------------------------------------ |
2 | // Implementation of AliPerformanceMC 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 member of AliPerformanceMC. | |
9 | // | |
10 | // Author: J.Otwinowski 04/02/2008 | |
11 | //------------------------------------------------------------------------------ | |
12 | ||
13 | /* | |
14 | ||
15 | // after running comparison task, read the file, and get component | |
2bfe5463 | 16 | gROOT->LoadMacro("$ALICE_ROOT/PWGPP/Macros/LoadMyLibs.C"); |
7cc34f08 | 17 | LoadMyLibs(); |
18 | ||
19 | TFile f("Output.root"); | |
20 | AliPerformanceMC * compObj = (AliPerformanceMC*)coutput->FindObject("AliPerformanceMC"); | |
21 | ||
22 | // analyse comparison data | |
23 | compObj->Analyse(); | |
24 | ||
25 | // the output histograms/graphs will be stored in the folder "folderRes" | |
26 | compObj->GetAnalysisFolder()->ls("*"); | |
27 | ||
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_MC.root","recreate"); | |
31 | compObj->Write(); // compObj->GetAnalysisFolder()->Write(); | |
32 | fout.Close(); | |
33 | ||
34 | */ | |
35 | ||
36 | #include "TCanvas.h" | |
37 | #include "TH1.h" | |
38 | #include "TH2.h" | |
39 | #include "TAxis.h" | |
40 | #include "TF1.h" | |
41 | ||
42 | #include "AliPerformanceMC.h" | |
43 | #include "AliESDEvent.h" | |
44 | #include "AliESDVertex.h" | |
45 | #include "AliESDtrack.h" | |
46 | #include "AliESDfriendTrack.h" | |
47 | #include "AliESDfriend.h" | |
48 | #include "AliLog.h" | |
49 | #include "AliMCEvent.h" | |
50 | #include "AliMCParticle.h" | |
51 | #include "AliHeader.h" | |
52 | #include "AliGenEventHeader.h" | |
53 | #include "AliStack.h" | |
54 | #include "AliMCInfoCuts.h" | |
55 | #include "AliRecInfoCuts.h" | |
56 | #include "AliTracker.h" | |
57 | #include "AliTreeDraw.h" | |
58 | #include "AliTPCseed.h" | |
59 | ||
60 | using namespace std; | |
61 | ||
62 | ClassImp(AliPerformanceMC) | |
63 | ||
64 | //_____________________________________________________________________________ | |
4823e4d4 | 65 | AliPerformanceMC::AliPerformanceMC(const Char_t* name, const Char_t* title, Int_t analysisMode, Bool_t hptGenerator): |
7cc34f08 | 66 | AliPerformanceObject(name,title), |
67 | fResolHisto(0), | |
68 | fPullHisto(0), | |
69 | ||
70 | // Cuts | |
71 | fCutsRC(0), | |
72 | fCutsMC(0), | |
73 | ||
74 | // histogram folder | |
75 | fAnalysisFolder(0) | |
76 | { | |
77 | // named constructor | |
78 | // | |
79 | SetAnalysisMode(analysisMode); | |
80 | SetHptGenerator(hptGenerator); | |
81 | ||
82 | Init(); | |
83 | } | |
84 | ||
85 | //_____________________________________________________________________________ | |
86 | AliPerformanceMC::~AliPerformanceMC() | |
87 | { | |
88 | // destructor | |
89 | ||
90 | if(fResolHisto) delete fResolHisto; fResolHisto=0; | |
91 | if(fPullHisto) delete fPullHisto; fPullHisto=0; | |
92 | ||
93 | if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0; | |
94 | } | |
95 | ||
96 | //_____________________________________________________________________________ | |
97 | void AliPerformanceMC::Init(){ | |
98 | ||
99 | // | |
100 | // histogram bining | |
101 | // | |
102 | ||
103 | // set pt bins | |
104 | Int_t nPtBins = 50; | |
ef1ed531 | 105 | Double_t ptMin = 1.e-2, ptMax = 20.; |
7cc34f08 | 106 | |
107 | Double_t *binsPt = 0; | |
ef1ed531 | 108 | |
7cc34f08 | 109 | if (IsHptGenerator()) { |
ef1ed531 | 110 | ptMax = 100.; |
111 | } | |
112 | binsPt = CreateLogAxis(nPtBins,ptMin,ptMax); | |
113 | ||
7cc34f08 | 114 | |
115 | Double_t yMin = -0.02, yMax = 0.02; | |
116 | Double_t zMin = -12.0, zMax = 12.0; | |
117 | if(GetAnalysisMode() == 3) { // TrackRef coordinate system | |
118 | yMin = -100.; yMax = 100.; | |
119 | zMin = -100.; zMax = 100.; | |
120 | } | |
121 | ||
122 | // res_y:res_z:res_phi,res_lambda:res_pt:y:z:eta:phi:pt | |
123 | Int_t binsResolHisto[10]={100,100,100,100,100,25,50,90,30,nPtBins}; | |
124 | Double_t minResolHisto[10]={-0.2,-0.2,-0.002,-0.002,-0.002, yMin, zMin, 0., -1.5, ptMin}; | |
125 | Double_t maxResolHisto[10]={ 0.2, 0.2, 0.002, 0.002, 0.002, yMax, zMax, 2.*TMath::Pi(), 1.5, ptMax}; | |
126 | ||
127 | fResolHisto = new THnSparseF("fResolHisto","res_y:res_z:res_phi:res_lambda:res_pt:y:z:phi:eta:pt",10,binsResolHisto,minResolHisto,maxResolHisto); | |
128 | fResolHisto->SetBinEdges(9,binsPt); | |
129 | ||
130 | fResolHisto->GetAxis(0)->SetTitle("y-y_{mc} (cm)"); | |
131 | fResolHisto->GetAxis(1)->SetTitle("z-z_{mc} (cm)"); | |
132 | fResolHisto->GetAxis(2)->SetTitle("#phi-#phi_{mc} (rad)"); | |
133 | fResolHisto->GetAxis(3)->SetTitle("#lambda-#lambda_{mc} (rad)"); | |
134 | fResolHisto->GetAxis(4)->SetTitle("(p_{T}/p_{Tmc}-1)"); | |
135 | fResolHisto->GetAxis(5)->SetTitle("y_{mc} (cm)"); | |
136 | fResolHisto->GetAxis(6)->SetTitle("z_{mc} (cm)"); | |
137 | fResolHisto->GetAxis(7)->SetTitle("#phi_{mc} (rad)"); | |
138 | fResolHisto->GetAxis(8)->SetTitle("#eta_{mc}"); | |
139 | fResolHisto->GetAxis(9)->SetTitle("p_{Tmc} (GeV/c)"); | |
140 | fResolHisto->Sumw2(); | |
141 | ||
142 | //pull_y:pull_z:pull_snp:pull_tgl:pull_1pt:y:z:snp:tgl:1pt | |
143 | Int_t binsPullHisto[10]={100,100,100,100,100,50,50,50,50,nPtBins}; | |
144 | Double_t minPullHisto[10]={-5.,-5.,-5.,-5.,-5., yMin, zMin,-1., -2.0, ptMin}; | |
145 | Double_t maxPullHisto[10]={ 5., 5., 5., 5., 5., yMax, zMax, 1., 2.0, ptMax}; | |
146 | fPullHisto = new THnSparseF("fPullHisto","pull_y:pull_z:pull_y:pull_z:pull_snp:pull_tgl:pull_1pt:y:z:snp:tgl:1pt",10,binsPullHisto,minPullHisto,maxPullHisto); | |
147 | ||
148 | fPullHisto->GetAxis(0)->SetTitle("(y-y_{mc})/#sigma"); | |
149 | fPullHisto->GetAxis(1)->SetTitle("(z-z_{mc})/#sigma"); | |
150 | fPullHisto->GetAxis(2)->SetTitle("(sin#phi-sin#phi_{mc})/#sigma"); | |
151 | fPullHisto->GetAxis(3)->SetTitle("(tan#lambda-tan#lambda_{mc})/#sigma"); | |
152 | fPullHisto->GetAxis(4)->SetTitle("(p_{Tmc}/p_{T}-1)/#sigma"); | |
153 | fPullHisto->GetAxis(5)->SetTitle("y_{mc} (cm)"); | |
154 | fPullHisto->GetAxis(6)->SetTitle("z_{mc} (cm)"); | |
155 | fPullHisto->GetAxis(7)->SetTitle("sin#phi_{mc}"); | |
156 | fPullHisto->GetAxis(8)->SetTitle("tan#lambda_{mc}"); | |
157 | fPullHisto->GetAxis(9)->SetTitle("1/p_{Tmc} (GeV/c)^{-1}"); | |
158 | fPullHisto->Sumw2(); | |
159 | ||
160 | // Init cuts | |
161 | if(!fCutsMC) | |
162 | AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object"); | |
163 | if(!fCutsRC) | |
164 | AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object"); | |
165 | ||
166 | // init folder | |
167 | fAnalysisFolder = CreateFolder("folderMC","Analysis MC Folder"); | |
168 | } | |
169 | ||
170 | //_____________________________________________________________________________ | |
171 | void AliPerformanceMC::Exec(AliMCEvent* const mcEvent, AliESDEvent *const /*esdEvent*/, AliESDfriend *const /*esdFriend*/, const Bool_t bUseMC, const Bool_t /*bUseESDfriend*/) | |
172 | { | |
173 | // Process pure MC information | |
174 | // | |
175 | AliHeader* header = 0; | |
176 | AliGenEventHeader* genHeader = 0; | |
177 | AliStack* stack = 0; | |
178 | TArrayF vtxMC(3); | |
179 | ||
180 | if(bUseMC) | |
181 | { | |
182 | if(!mcEvent) { | |
183 | Error("Exec","mcEvent not available"); | |
184 | return; | |
185 | } | |
186 | // get MC event header | |
187 | header = mcEvent->Header(); | |
188 | if (!header) { | |
189 | Error("Exec","Header not available"); | |
190 | return; | |
191 | } | |
192 | // MC particle stack | |
193 | stack = mcEvent->Stack(); | |
194 | if (!stack) { | |
195 | Error("Exec","Stack not available"); | |
196 | return; | |
197 | } | |
198 | // get MC vertex | |
199 | genHeader = header->GenEventHeader(); | |
200 | if (!genHeader) { | |
201 | Error("Exec","Could not retrieve genHeader from Header"); | |
202 | return; | |
203 | } | |
204 | genHeader->PrimaryVertex(vtxMC); | |
205 | } | |
206 | else { | |
207 | Error("Exec","MC information required!"); | |
208 | return; | |
209 | } | |
210 | ||
211 | Int_t nPart = mcEvent->GetNumberOfTracks(); | |
212 | if (nPart==0) return; | |
213 | ||
214 | //TParticle * part = new TParticle; | |
215 | //TClonesArray * trefs = new TClonesArray("AliTrackReference"); | |
216 | TParticle *part=0; | |
217 | TClonesArray *trefs=0; | |
218 | ||
219 | // Process events | |
220 | for (Int_t iPart = 0; iPart < nPart; iPart++) | |
221 | { | |
222 | Int_t status = mcEvent->GetParticleAndTR(iPart, part, trefs); | |
223 | if (status<0) continue; | |
224 | if(!part) continue; | |
225 | if(!trefs) continue; | |
226 | ||
227 | AliMCParticle *mcParticle = (AliMCParticle*) mcEvent->GetTrack(iPart); | |
228 | if(!mcParticle) continue; | |
229 | ||
230 | TParticle *particle = mcParticle->Particle(); | |
231 | if(!particle) continue; | |
232 | ||
233 | Bool_t prim = stack->IsPhysicalPrimary(iPart); | |
234 | ||
235 | // Only 5 charged particle species (e,mu,pi,K,p) | |
236 | if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) break; | |
237 | ||
238 | // exclude electrons | |
239 | if (fCutsMC->GetEP()==TMath::Abs(particle->GetPdgCode())) return; | |
240 | ||
241 | AliTrackReference *refTPCIn = 0; | |
242 | AliTrackReference *refTPCOut = 0; | |
243 | AliTrackReference *refITSIn = 0; | |
244 | AliTrackReference *refITSOut = 0; | |
245 | AliTrackReference *refTRDIn = 0; | |
246 | AliTrackReference *refTRDOut = 0; | |
247 | ||
248 | Float_t rmax = 0.; | |
249 | Int_t nTrackRef = mcParticle->GetNumberOfTrackReferences(); | |
250 | if(nTrackRef < 5) continue; | |
251 | ||
252 | for (Int_t iref = 0; iref < nTrackRef; iref++) | |
253 | { | |
254 | AliTrackReference *ref = mcParticle->GetTrackReference(iref); | |
255 | if(!ref) continue; | |
256 | ||
257 | Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py(); | |
258 | if(dir < 0.) break; // looping tracks (return back) | |
259 | if (ref->R()<rmax) break; | |
260 | ||
261 | // pt threshold | |
262 | if(ref->Pt() < fCutsRC->GetPtMin()) break; | |
263 | ||
264 | // TPC | |
265 | if(ref->DetectorId()==AliTrackReference::kTPC) | |
266 | { | |
267 | if(!refTPCIn) { | |
268 | refTPCIn = ref; | |
269 | } | |
270 | else { | |
271 | if(ref->R() > refTPCIn->R()) | |
272 | refTPCOut = ref; | |
273 | } | |
274 | } | |
275 | // ITS | |
276 | if(ref->DetectorId()==AliTrackReference::kITS) | |
277 | { | |
278 | if(!refITSIn) { | |
279 | refITSIn = ref; | |
280 | } | |
281 | else { | |
282 | if(ref->R() > refITSIn->R()) | |
283 | refITSOut = ref; | |
284 | } | |
285 | } | |
286 | // TRD | |
287 | if(ref->DetectorId()==AliTrackReference::kTRD) | |
288 | { | |
289 | if(!refTRDIn) { | |
290 | refTRDIn = ref; | |
291 | } | |
292 | else { | |
293 | if(ref->R() > refTRDIn->R()) | |
294 | refTRDOut = ref; | |
295 | } | |
296 | } | |
297 | if (ref->R()>rmax) rmax=ref->R(); | |
298 | } | |
299 | ||
300 | if(GetAnalysisMode() == 0) { | |
301 | // Propagate refTPCIn to DCA to prim. vertex and make comparison | |
302 | // only primary particles | |
303 | if(refTPCIn) { | |
304 | if(prim) ProcessTPC(refTPCIn,particle); | |
305 | } | |
306 | } | |
307 | ||
308 | if(GetAnalysisMode() == 3) { | |
309 | // Propagate refTPCOut to refTPCIn using AliTrack::PropagateTrackToBxByBz() | |
310 | // | |
311 | if(refTPCIn && refTPCOut) | |
312 | ProcessInnerTPC(refTPCIn,refTPCOut,particle); | |
313 | } | |
314 | ||
315 | if(GetAnalysisMode() == 4) { | |
316 | // propagate refTPCIn to refTPCOut using AliExternalTrackParam::PropagateToBxByBz() | |
317 | ProcessOuterTPCExt(part,trefs); | |
318 | } | |
319 | } | |
320 | ||
321 | //if(part) delete part; | |
322 | //if(trefs) delete trefs; | |
323 | } | |
324 | ||
325 | //_____________________________________________________________________________ | |
326 | void AliPerformanceMC::ProcessTPC(AliTrackReference* const refIn, TParticle *const particle) { | |
327 | // | |
328 | // Test propagation from refIn to DCA | |
329 | // | |
330 | if(!refIn) return; | |
331 | if(!particle) return; | |
332 | ||
333 | AliExternalTrackParam *track = 0; | |
334 | Double_t radius = 2.8; // cm | |
335 | Double_t mass = particle->GetMass(); | |
336 | Double_t step=1; | |
337 | // | |
338 | track=MakeTrack(refIn,particle); | |
339 | if (!track) return; | |
340 | // | |
341 | //AliTracker::PropagateTrackTo(track, radius+step, mass, step, kTRUE,0.99); | |
342 | //AliTracker::PropagateTrackTo(track, radius+0.5, mass, step*0.1, kTRUE,0.99); | |
343 | AliTracker::PropagateTrackToBxByBz(track, radius+step, mass, step, kTRUE,0.99); | |
344 | AliTracker::PropagateTrackToBxByBz(track, radius+0.5, mass, step*0.1, kTRUE,0.99); | |
345 | Double_t xyz[3] = {particle->Vx(),particle->Vy(),particle->Vz()}; | |
346 | Double_t sxyz[3]={0.0,0.0,0.0}; | |
347 | AliESDVertex vertex(xyz,sxyz); | |
348 | Double_t dca[2], cov[3]; | |
349 | //Bool_t isOK = track->PropagateToDCA(&vertex,AliTracker::GetBz(),10,dca,cov); | |
350 | Double_t field[3]; track->GetBxByBz(field); | |
351 | Bool_t isOK = track->PropagateToDCABxByBz(&vertex,field,10,dca,cov); | |
352 | if(!isOK) return; | |
353 | ||
354 | // Fill histogram | |
355 | Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC; | |
356 | Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC; | |
357 | ||
358 | Float_t mceta = particle->Eta(); | |
359 | Float_t mcphi = particle->Phi(); | |
360 | if(mcphi<0) mcphi += 2.*TMath::Pi(); | |
361 | Float_t mcpt = particle->Pt(); | |
362 | Float_t mcsnp = TMath::Sin(TMath::ATan2(particle->Py(),particle->Px())); | |
363 | Float_t mctgl = TMath::Tan(TMath::ATan2(particle->Pz(),particle->Pt())); | |
364 | ||
365 | //if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ()) | |
366 | //{ | |
367 | if(mcpt == 0) return; | |
368 | ||
369 | //deltaYTPC= track->GetY()-particle->Vy(); | |
370 | deltaYTPC= track->GetY(); // it is already rotated | |
371 | deltaZTPC = track->GetZ()-particle->Vz(); | |
372 | deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(particle->Pz(),particle->Pt()); | |
373 | deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(particle->Py(),particle->Px()); | |
374 | deltaPtTPC = (track->Pt()-mcpt) / mcpt; | |
375 | ||
376 | /* | |
377 | printf("------------------------------ \n"); | |
378 | printf("trY %f, trZ %f, trTheta %f ,trPhi %f, trPt %f \n", track->GetY(), track->GetZ(), TMath::ATan2(track->Pz(),track->Pt()), TMath::ATan2(track->Py(),track->Px()), track->Pt() ); | |
379 | printf("partY %f, partZ %f, partTheta %f ,partPhi %f, partPt %f \n", particle->Vy(), particle->Vz(), TMath::ATan2(particle->Pz(),particle->Pt()), TMath::ATan2(particle->Py(),particle->Px()), mcpt ); | |
380 | */ | |
381 | ||
382 | ||
383 | pullYTPC= (track->GetY()-particle->Vy()) / TMath::Sqrt(track->GetSigmaY2()); | |
384 | pullZTPC = (track->GetZ()-particle->Vz()) / TMath::Sqrt(track->GetSigmaZ2()); | |
385 | ||
386 | pullPhiTPC = (track->GetSnp() - mcsnp) / TMath::Sqrt(track->GetSigmaSnp2()); | |
387 | pullLambdaTPC = (track->GetTgl() - mctgl) / TMath::Sqrt(track->GetSigmaTgl2()); | |
388 | ||
389 | if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2()); | |
390 | else pull1PtTPC = 0.; | |
391 | ||
392 | Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,particle->Vy(),particle->Vz(),mcphi,mceta,mcpt}; | |
393 | fResolHisto->Fill(vResolHisto); | |
394 | ||
395 | Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mcsnp,mctgl,1./mcpt}; | |
396 | fPullHisto->Fill(vPullHisto); | |
397 | //} | |
398 | if(track) delete track; | |
399 | } | |
400 | ||
401 | //_____________________________________________________________________________ | |
402 | void AliPerformanceMC::ProcessInnerTPC(AliTrackReference* const refIn, AliTrackReference *const refOut, TParticle *const particle) { | |
403 | // | |
404 | // Test propagation from Out to In | |
405 | // | |
406 | if(!refIn) return; | |
407 | if(!refOut) return; | |
408 | if(!particle) return; | |
409 | ||
410 | AliExternalTrackParam *track=MakeTrack(refOut,particle); | |
411 | if (!track) return; | |
412 | ||
413 | Double_t mass = particle->GetMass(); | |
414 | Double_t step=1; | |
415 | ||
416 | Float_t alphaIn= TMath::ATan2(refIn->Y(),refIn->X()); | |
417 | //Float_t alphaOut= TMath::ATan2(refOut->Y(),refOut->X()); | |
418 | //track->Rotate(alphaIn-alphaOut); | |
419 | track->Rotate(alphaIn); | |
420 | //printf("alphaIn %f, alphaOut %f \n",alphaIn,alphaOut); | |
421 | ||
422 | // | |
423 | //Bool_t isOK = AliTracker::PropagateTrackTo(track, refIn->R(), mass, step, kTRUE,0.99); | |
424 | //Bool_t isOK = AliTracker::PropagateTrackTo(track, refIn->R(), mass, step, kFALSE,0.99); | |
425 | Bool_t isOK = AliTracker::PropagateTrackToBxByBz(track, refIn->R(), mass, step, kFALSE,0.99); | |
426 | if(!isOK) return; | |
427 | ||
428 | // calculate alpha angle | |
429 | //Double_t xyz[3] = {refIn->X(),refIn->Y(),refIn0->Z()}; | |
430 | //Double_t alpha = TMath::ATan2(xyz[1],xyz[0]); | |
431 | //if (alpha<0) alpha += TMath::TwoPi(); | |
432 | ||
433 | // rotate outer track to inner | |
434 | // and propagate to the radius of the first track referenco of TPC | |
435 | //Double_t trRadius = TMath::Sqrt(xyz[1] * xyz[1] + xyz[0] * xyz[0]); | |
436 | //Bool_t isOK = track->Propagate(alpha,trRadius,AliTracker::GetBz()); | |
437 | //if(!isOK) return; | |
438 | ||
439 | Float_t mceta = -TMath::Log(TMath::Tan(0.5 * refIn->Theta())); | |
440 | Float_t mcphi = refIn->Phi(); | |
441 | if(mcphi<0) mcphi += 2.*TMath::Pi(); | |
442 | Float_t mcpt = refIn->Pt(); | |
443 | Float_t mcsnp = TMath::Sin(TMath::ATan2(refIn->Py(),refIn->Px())); | |
444 | Float_t mctgl = TMath::Tan(TMath::ATan2(refIn->Pz(),refIn->Pt())); | |
445 | ||
446 | Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC; | |
447 | Float_t pull1PtTPC=0., pullYTPC=0., pullZTPC=0., pullPhiTPC=0., pullLambdaTPC=0.; | |
448 | ||
449 | if(mcpt == 0) return; | |
450 | ||
451 | //deltaYTPC = track->GetY()-refIn->Y(); | |
452 | deltaYTPC = track->GetY(); | |
453 | deltaZTPC = track->GetZ()-refIn->Z(); | |
454 | deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(refIn->Pz(),refIn->Pt()); | |
455 | deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(refIn->Py(),refIn->Px()); | |
456 | deltaPtTPC = (track->Pt()-mcpt) / mcpt; | |
457 | ||
458 | /* | |
459 | printf("------------------------------ \n"); | |
460 | printf("trX %f, trY %f, trZ %f, trTheta %f ,trPhi %f, trPt %f \n", track->GetX(), track->GetY(), track->GetZ(), TMath::ATan2(track->Pz(),track->Pt()), TMath::ATan2(track->Py(),track->Px()), track->Pt() ); | |
461 | printf("partX %f, partY %f, partZ %f, partTheta %f ,partPhi %f, partPt %f \n",refIn->X(), refIn->Y(), refIn->Z(), TMath::ATan2(refIn->Pz(),refIn->Pt()), TMath::ATan2(refIn->Py(),refIn->Px()), mcpt ); | |
462 | */ | |
463 | ||
464 | if(TMath::Sqrt(track->GetSigmaY2())) pullYTPC = track->GetY() / TMath::Sqrt(track->GetSigmaY2()); | |
465 | if(TMath::Sqrt(track->GetSigmaZ2())) pullZTPC = (track->GetZ()-refIn->Z()) / TMath::Sqrt(track->GetSigmaZ2()); | |
466 | if(TMath::Sqrt(track->GetSigmaSnp2())) pullPhiTPC = (track->GetSnp() - mcsnp) / TMath::Sqrt(track->GetSigmaSnp2()); | |
467 | if(TMath::Sqrt(track->GetSigmaTgl2())) pullLambdaTPC = (track->GetTgl() - mctgl) / TMath::Sqrt(track->GetSigmaTgl2()); | |
468 | if(TMath::Sqrt(track->GetSigma1Pt2())) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2()); | |
469 | ||
470 | //printf("pullYTPC %f,pullZTPC %f,pullPhiTPC %f,pullLambdaTPC %f,pull1PtTPC %f,refIn->Y() %f,refIn->Z() %f,mcsnp %f,mctgl %f,1./mcpt %f \n",pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,refIn->Y(),refIn->Z(),mcsnp,mctgl,1./mcpt); | |
471 | ||
472 | Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,refIn->Y(),refIn->Z(),mcphi,mceta,mcpt}; | |
473 | fResolHisto->Fill(vResolHisto); | |
474 | ||
475 | Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,refIn->Y(),refIn->Z(),mcsnp,mctgl,1./mcpt}; | |
476 | fPullHisto->Fill(vPullHisto); | |
477 | ||
478 | if(track) delete track; | |
479 | } | |
480 | ||
481 | //_____________________________________________________________________________ | |
482 | void AliPerformanceMC::ProcessOuterTPCExt(TParticle *const part, TClonesArray * const trefs) | |
483 | { | |
484 | const Int_t kMinRefs=5; | |
485 | Int_t nrefs = trefs->GetEntries(); | |
486 | if (nrefs<kMinRefs) return; // we should have enough references | |
487 | Int_t iref0 =-1; | |
488 | Int_t iref1 =-1; | |
489 | ||
490 | for (Int_t iref=0; iref<nrefs; iref++){ | |
491 | AliTrackReference * ref = (AliTrackReference*)trefs->At(iref); | |
492 | if (!ref) continue; | |
493 | Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py(); | |
494 | if (dir<0) break; | |
495 | if (ref->DetectorId()!=AliTrackReference::kTPC) continue; | |
496 | if (iref0<0) iref0 = iref; | |
497 | iref1 = iref; | |
498 | } | |
499 | if (iref1-iref0<kMinRefs) return; | |
c6d45264 | 500 | Double_t covar[15]; |
501 | for (Int_t icov=0; icov<15; icov++) covar[icov]=0; | |
7cc34f08 | 502 | covar[0]=1; |
503 | covar[2]=1; | |
504 | covar[5]=1; | |
505 | covar[9]=1; | |
506 | covar[14]=1; | |
507 | ||
508 | AliTrackReference * refIn = (AliTrackReference*)trefs->At(iref0); | |
509 | AliTrackReference * refOut = (AliTrackReference*)trefs->At(iref1); | |
510 | AliExternalTrackParam *paramPropagate= MakeTrack(refIn,part); | |
511 | AliExternalTrackParam *paramUpdate = MakeTrack(refIn,part); | |
512 | paramUpdate->AddCovariance(covar); | |
513 | ||
514 | Bool_t isOKP=kTRUE; | |
515 | Bool_t isOKU=kTRUE; | |
516 | AliMagF * field = (AliMagF*) TGeoGlobalMagField::Instance()->GetField(); | |
517 | if(!field) return; | |
518 | for (Int_t iref = iref0; iref<=iref1; iref++){ | |
519 | AliTrackReference * ref = (AliTrackReference*)trefs->At(iref); | |
520 | if(!ref) continue; | |
521 | Float_t alphaC= TMath::ATan2(ref->Y(),ref->X()); | |
522 | Double_t pos[3] = {ref->X(), ref->Y(), ref->Z()}; | |
523 | Double_t mag[3]; | |
524 | field->Field(pos,mag); | |
525 | isOKP&=paramPropagate->Rotate(alphaC); | |
526 | isOKU&=paramUpdate->Rotate(alphaC); | |
527 | for (Float_t xref= paramPropagate->GetX(); xref<ref->R(); xref++){ | |
528 | //isOKP&=paramPropagate->PropagateTo(xref, mag[2]); | |
529 | //isOKU&=paramUpdate->PropagateTo(xref, mag[2]); | |
530 | isOKP&=paramPropagate->PropagateToBxByBz(xref, mag); | |
531 | isOKU&=paramUpdate->PropagateToBxByBz(xref, mag); | |
532 | } | |
533 | //isOKP&=paramPropagate->PropagateTo(ref->R(), mag[2]); | |
534 | //isOKU&=paramUpdate->PropagateTo(ref->R(), mag[2]); | |
535 | isOKP&=paramPropagate->PropagateToBxByBz(ref->R(), mag); | |
536 | isOKU&=paramUpdate->PropagateToBxByBz(ref->R(), mag); | |
537 | Double_t clpos[2] = {0, ref->Z()}; | |
538 | Double_t clcov[3] = { 0.005,0,0.005}; | |
539 | isOKU&= paramUpdate->Update(clpos, clcov); | |
540 | } | |
541 | ||
542 | Float_t mceta = -TMath::Log(TMath::Tan(0.5 * refOut->Theta())); | |
543 | Float_t mcphi = refOut->Phi(); | |
544 | if(mcphi<0) mcphi += 2.*TMath::Pi(); | |
545 | Float_t mcpt = refOut->Pt(); | |
546 | Float_t mcsnp = TMath::Sin(TMath::ATan2(refOut->Py(),refOut->Px())); | |
547 | Float_t mctgl = TMath::Tan(TMath::ATan2(refOut->Pz(),refOut->Pt())); | |
548 | ||
549 | Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC; | |
550 | Float_t pull1PtTPC=0., pullYTPC=0., pullZTPC=0., pullPhiTPC=0., pullLambdaTPC=0.; | |
551 | ||
552 | if(mcpt == 0) return; | |
553 | ||
554 | //deltaYTPC = track->GetY()-refIn->Y(); | |
555 | deltaYTPC = paramPropagate->GetY(); | |
556 | deltaZTPC = paramPropagate->GetZ()-refOut->Z(); | |
557 | deltaLambdaTPC = TMath::ATan2(paramPropagate->Pz(),paramPropagate->Pt())-TMath::ATan2(refOut->Pz(),refOut->Pt()); | |
558 | deltaPhiTPC = TMath::ATan2(paramPropagate->Py(),paramPropagate->Px())-TMath::ATan2(refOut->Py(),refOut->Px()); | |
559 | deltaPtTPC = (paramPropagate->Pt()-mcpt) / mcpt; | |
560 | ||
561 | if(TMath::Sqrt(paramPropagate->GetSigmaY2())) pullYTPC = paramPropagate->GetY() / TMath::Sqrt(paramPropagate->GetSigmaY2()); | |
562 | if(TMath::Sqrt(paramPropagate->GetSigmaZ2())) pullZTPC = (paramPropagate->GetZ()-refOut->Z()) / TMath::Sqrt(paramPropagate->GetSigmaZ2()); | |
563 | if(TMath::Sqrt(paramPropagate->GetSigmaSnp2())) pullPhiTPC = (paramPropagate->GetSnp() - mcsnp) / TMath::Sqrt(paramPropagate->GetSigmaSnp2()); | |
564 | if(TMath::Sqrt(paramPropagate->GetSigmaTgl2())) pullLambdaTPC = (paramPropagate->GetTgl() - mctgl) / TMath::Sqrt(paramPropagate->GetSigmaTgl2()); | |
565 | if(TMath::Sqrt(paramPropagate->GetSigma1Pt2())) pull1PtTPC = (paramPropagate->OneOverPt()-1./mcpt) / TMath::Sqrt(paramPropagate->GetSigma1Pt2()); | |
566 | ||
567 | //printf("pullYTPC %f,pullZTPC %f,pullPhiTPC %f,pullLambdaTPC %f,pull1PtTPC %f,refIn->Y() %f,refIn->Z() %f,mcsnp %f,mctgl %f,1./mcpt %f \n",pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,refIn->Y(),refIn->Z(),mcsnp,mctgl,1./mcpt); | |
568 | ||
569 | Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,refIn->Y(),refIn->Z(),mcphi,mceta,mcpt}; | |
570 | fResolHisto->Fill(vResolHisto); | |
571 | ||
572 | Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,refIn->Y(),refIn->Z(),mcsnp,mctgl,1./mcpt}; | |
573 | fPullHisto->Fill(vPullHisto); | |
574 | } | |
575 | ||
576 | //_____________________________________________________________________________ | |
577 | AliExternalTrackParam * AliPerformanceMC::MakeTrack(const AliTrackReference* const ref, TParticle* const part) | |
578 | { | |
579 | // | |
580 | // Make track out of the track ref | |
581 | // part - TParticle used to determine chargr | |
582 | // the covariance matrix - equal 0 - starting from ideal MC position | |
583 | Double_t xyz[3]={ref->X(),ref->Y(),ref->Z()}; | |
584 | Double_t pxyz[3]={ref->Px(),ref->Py(),ref->Pz()}; | |
585 | Double_t cv[21]; | |
586 | for (Int_t i=0; i<21;i++) cv[i]=0; | |
587 | if (!part->GetPDG()) return 0; | |
588 | AliExternalTrackParam * param = new AliExternalTrackParam(xyz,pxyz,cv,TMath::Nint(part->GetPDG()->Charge()/3.)); | |
589 | return param; | |
590 | } | |
591 | ||
592 | //_____________________________________________________________________________ | |
593 | TH1F* AliPerformanceMC::MakeResol(TH2F * his, Int_t integ, Bool_t type, Int_t cut){ | |
594 | // Create resolution histograms | |
595 | ||
596 | TH1F *hisr, *hism; | |
597 | if (!gPad) new TCanvas; | |
598 | hisr = AliTreeDraw::CreateResHistoII(his,&hism,integ,kTRUE,cut); | |
599 | if (type) return hism; | |
600 | else | |
601 | return hisr; | |
602 | } | |
603 | ||
604 | //_____________________________________________________________________________ | |
605 | void AliPerformanceMC::Analyse() { | |
606 | // Analyse comparison information and store output histograms | |
607 | // in the folder "folderRes" | |
608 | // | |
609 | TH1::AddDirectory(kFALSE); | |
610 | TH1F *h=0; | |
611 | TH2F *h2D=0; | |
612 | TObjArray *aFolderObj = new TObjArray; | |
18e588e9 | 613 | if(!aFolderObj) return; |
7cc34f08 | 614 | |
615 | // write results in the folder | |
616 | TCanvas * c = new TCanvas("Phi resol Tan","Phi resol Tan"); | |
617 | c->cd(); | |
618 | ||
619 | char name[256]; | |
620 | char title[256]; | |
621 | ||
622 | for(Int_t i=0; i<5; i++) | |
623 | { | |
624 | for(Int_t j=5; j<10; j++) | |
625 | { | |
626 | //if(j!=8) fResolHisto->GetAxis(8)->SetRangeUser(-0.9,0.89); // eta window | |
627 | fResolHisto->GetAxis(9)->SetRangeUser(0.16,10.); // pt threshold | |
628 | if(j!=8) fResolHisto->GetAxis(8)->SetRangeUser(0.,0.9); // eta window | |
629 | //fResolHisto->GetAxis(9)->SetRangeUser(0.16,3.); // pt window | |
630 | if(GetAnalysisMode() == 3) fResolHisto->GetAxis(5)->SetRangeUser(-80.,80.); // y range | |
631 | ||
632 | h2D = (TH2F*)fResolHisto->Projection(i,j); | |
633 | ||
634 | h = AliPerformanceMC::MakeResol(h2D,1,0,20); | |
18e588e9 | 635 | snprintf(name,256,"h_res_%d_vs_%d",i,j); |
7cc34f08 | 636 | h->SetName(name); |
637 | ||
638 | h->GetXaxis()->SetTitle(fResolHisto->GetAxis(j)->GetTitle()); | |
18e588e9 | 639 | snprintf(title,256,"%s %s",fResolHisto->GetAxis(i)->GetTitle(),"(resolution)"); |
7cc34f08 | 640 | h->GetYaxis()->SetTitle(title); |
18e588e9 | 641 | snprintf(title,256,"%s vs %s",title,fResolHisto->GetAxis(j)->GetTitle()); |
7cc34f08 | 642 | h->SetTitle(title); |
643 | ||
644 | if(j==9) h->SetBit(TH1::kLogX); | |
645 | aFolderObj->Add(h); | |
646 | ||
647 | h = AliPerformanceMC::MakeResol(h2D,1,1,20); | |
648 | //h = (TH1F*)arr->At(1); | |
18e588e9 | 649 | snprintf(name,256,"h_mean_res_%d_vs_%d",i,j); |
7cc34f08 | 650 | h->SetName(name); |
651 | ||
652 | h->GetXaxis()->SetTitle(fResolHisto->GetAxis(j)->GetTitle()); | |
18e588e9 | 653 | snprintf(title,256,"%s %s",fResolHisto->GetAxis(i)->GetTitle(),"(mean)"); |
7cc34f08 | 654 | h->GetYaxis()->SetTitle(title); |
655 | ||
18e588e9 | 656 | snprintf(title,256,"%s vs %s",title,fResolHisto->GetAxis(j)->GetTitle()); |
7cc34f08 | 657 | h->SetTitle(title); |
658 | ||
659 | if(j==9) h->SetBit(TH1::kLogX); | |
660 | aFolderObj->Add(h); | |
661 | ||
662 | fResolHisto->GetAxis(8)->SetRangeUser(-1.5,1.5); | |
663 | fResolHisto->GetAxis(9)->SetRangeUser(0.0,10.); | |
664 | ||
665 | // | |
666 | if(j!=8) fPullHisto->GetAxis(8)->SetRangeUser(-1.0,0.99); // theta window | |
667 | else fPullHisto->GetAxis(8)->SetRangeUser(-1.5,1.5); // theta window | |
668 | fPullHisto->GetAxis(9)->SetRangeUser(0.0,10.); // pt threshold | |
669 | if(GetAnalysisMode() == 3) fPullHisto->GetAxis(5)->SetRangeUser(-80.,80.); // y range | |
670 | ||
671 | h2D = (TH2F*)fPullHisto->Projection(i,j); | |
672 | ||
673 | h = AliPerformanceMC::MakeResol(h2D,1,0,20); | |
18e588e9 | 674 | snprintf(name,256,"h_pull_%d_vs_%d",i,j); |
7cc34f08 | 675 | h->SetName(name); |
676 | ||
677 | h->GetXaxis()->SetTitle(fPullHisto->GetAxis(j)->GetTitle()); | |
18e588e9 | 678 | snprintf(title,256,"%s %s",fPullHisto->GetAxis(i)->GetTitle(),"(resolution)"); |
7cc34f08 | 679 | h->GetYaxis()->SetTitle(title); |
18e588e9 | 680 | snprintf(title,256,"%s vs %s",title,fPullHisto->GetAxis(j)->GetTitle()); |
7cc34f08 | 681 | h->SetTitle(title); |
682 | ||
683 | //if(j==9) h->SetBit(TH1::kLogX); | |
684 | aFolderObj->Add(h); | |
685 | ||
686 | h = AliPerformanceMC::MakeResol(h2D,1,1,20); | |
18e588e9 | 687 | snprintf(name,256,"h_mean_pull_%d_vs_%d",i,j); |
7cc34f08 | 688 | h->SetName(name); |
689 | ||
690 | h->GetXaxis()->SetTitle(fPullHisto->GetAxis(j)->GetTitle()); | |
18e588e9 | 691 | snprintf(title,256,"%s %s",fPullHisto->GetAxis(i)->GetTitle(),"(mean)"); |
7cc34f08 | 692 | h->GetYaxis()->SetTitle(title); |
18e588e9 | 693 | snprintf(title,256,"%s vs %s",title,fPullHisto->GetAxis(j)->GetTitle()); |
7cc34f08 | 694 | h->SetTitle(title); |
695 | ||
696 | //if(j==9) h->SetBit(TH1::kLogX); | |
697 | aFolderObj->Add(h); | |
698 | } | |
699 | } | |
700 | ||
701 | // export objects to analysis folder | |
702 | fAnalysisFolder = ExportToFolder(aFolderObj); | |
703 | ||
704 | // delete only TObjArray | |
705 | if(aFolderObj) delete aFolderObj; | |
706 | } | |
707 | ||
708 | //_____________________________________________________________________________ | |
709 | TFolder* AliPerformanceMC::ExportToFolder(TObjArray * array) | |
710 | { | |
711 | // recreate folder avery time and export objects to new one | |
712 | // | |
713 | AliPerformanceMC * comp=this; | |
714 | TFolder *folder = comp->GetAnalysisFolder(); | |
715 | ||
716 | TString name, title; | |
717 | TFolder *newFolder = 0; | |
718 | Int_t i = 0; | |
719 | Int_t size = array->GetSize(); | |
720 | ||
721 | if(folder) { | |
722 | // get name and title from old folder | |
723 | name = folder->GetName(); | |
724 | title = folder->GetTitle(); | |
725 | ||
726 | // delete old one | |
727 | delete folder; | |
728 | ||
729 | // create new one | |
730 | newFolder = CreateFolder(name.Data(),title.Data()); | |
731 | newFolder->SetOwner(); | |
732 | ||
733 | // add objects to folder | |
734 | while(i < size) { | |
735 | newFolder->Add(array->At(i)); | |
736 | i++; | |
737 | } | |
738 | } | |
739 | ||
740 | return newFolder; | |
741 | } | |
742 | ||
743 | //_____________________________________________________________________________ | |
744 | Long64_t AliPerformanceMC::Merge(TCollection* const list) | |
745 | { | |
746 | // Merge list of objects (needed by PROOF) | |
747 | ||
748 | if (!list) | |
749 | return 0; | |
750 | ||
751 | if (list->IsEmpty()) | |
752 | return 1; | |
753 | ||
754 | TIterator* iter = list->MakeIterator(); | |
755 | TObject* obj = 0; | |
756 | ||
757 | // collection of generated histograms | |
758 | Int_t count=0; | |
759 | while((obj = iter->Next()) != 0) | |
760 | { | |
761 | AliPerformanceMC* entry = dynamic_cast<AliPerformanceMC*>(obj); | |
762 | if (entry == 0) continue; | |
763 | ||
764 | fResolHisto->Add(entry->fResolHisto); | |
765 | fPullHisto->Add(entry->fPullHisto); | |
766 | ||
767 | count++; | |
768 | } | |
769 | ||
770 | return count; | |
771 | } | |
772 | ||
773 | //_____________________________________________________________________________ | |
774 | TFolder* AliPerformanceMC::CreateFolder(TString name,TString title) { | |
775 | // create folder for analysed histograms | |
776 | // | |
777 | TFolder *folder = 0; | |
778 | folder = new TFolder(name.Data(),title.Data()); | |
779 | ||
780 | return folder; | |
781 | } |