]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGPP/TPC/AliPerformanceMC.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
#102886: Various fixes to the the code in EVE, STEER, PWGPP, cmake
[u/mrichter/AliRoot.git] / PWGPP / TPC / AliPerformanceMC.cxx
CommitLineData
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
60using namespace std;
61
62ClassImp(AliPerformanceMC)
63
64//_____________________________________________________________________________
4823e4d4 65AliPerformanceMC::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//_____________________________________________________________________________
86AliPerformanceMC::~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//_____________________________________________________________________________
97void 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//_____________________________________________________________________________
171void 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//_____________________________________________________________________________
326void 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//_____________________________________________________________________________
402void 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//_____________________________________________________________________________
482void 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//_____________________________________________________________________________
577AliExternalTrackParam * 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//_____________________________________________________________________________
593TH1F* 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//_____________________________________________________________________________
605void 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//_____________________________________________________________________________
709TFolder* 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
740return newFolder;
741}
742
743//_____________________________________________________________________________
744Long64_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
770return count;
771}
772
773//_____________________________________________________________________________
774TFolder* AliPerformanceMC::CreateFolder(TString name,TString title) {
775// create folder for analysed histograms
776//
777TFolder *folder = 0;
778 folder = new TFolder(name.Data(),title.Data());
779
780 return folder;
781}