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