- classes for pT spectra charged hadrons analysis added
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / ChargedHadrons / dNdPt / AlidNdPtTrackDumpTask.cxx
1 /**************************************************************************\r
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
3 *                                                                        *\r
4 * Author: The ALICE Off-line Project.                                    *\r
5 * Contributors are mentioned in the code where appropriate.              *\r
6 *                                                                        *\r
7 * Permission to use, copy, modify and distribute this software and its   *\r
8 * documentation strictly for non-commercial purposes is hereby granted   *\r
9 * without fee, provided that the above copyright notice appears in all   *\r
10 * copies and that both the copyright notice and this permission notice   *\r
11 * appear in the supporting documentation. The authors make no claims     *\r
12 * about the suitability of this software for any purpose. It is          *\r
13 * provided "as is" without express or implied warranty.                  *\r
14 **************************************************************************/\r
15 \r
16 #include "iostream"\r
17 \r
18 #include <TPDGCode.h>\r
19 \r
20 #include "TChain.h"\r
21 #include "TTreeStream.h"\r
22 #include "TTree.h"\r
23 #include "TH1F.h"\r
24 #include "TCanvas.h"\r
25 #include "TList.h"\r
26 #include "TFile.h"\r
27 #include "TMatrixD.h"\r
28 #include "TRandom.h"\r
29 \r
30 #include "AliHeader.h"  \r
31 #include "AliGenEventHeader.h"  \r
32 #include "AliInputEventHandler.h"  \r
33 #include "AliStack.h"  \r
34 #include "AliTrackReference.h"  \r
35 \r
36 #include "AliPhysicsSelection.h"\r
37 #include "AliAnalysisTask.h"\r
38 #include "AliAnalysisManager.h"\r
39 #include "AliESDEvent.h"\r
40 #include "AliESDfriend.h"\r
41 #include "AliMCEvent.h"\r
42 #include "AliESDInputHandler.h"\r
43 #include "AliESDVertex.h"\r
44 #include "AliTracker.h"\r
45 #include "AliGeomManager.h"\r
46 \r
47 #include "AliCentrality.h"\r
48 #include "AliESDVZERO.h"\r
49 #include "AliMultiplicity.h"\r
50 \r
51 #include "AliESDtrackCuts.h"\r
52 #include "AliMCEventHandler.h"\r
53 #include "AlidNdPt.h"\r
54 #include "AlidNdPtEventCuts.h"\r
55 #include "AlidNdPtAcceptanceCuts.h"\r
56 \r
57 #include "AlidNdPtTrackDumpTask.h"\r
58 \r
59 using namespace std;\r
60 \r
61 ClassImp(AlidNdPtTrackDumpTask)\r
62 \r
63 //_____________________________________________________________________________\r
64 AlidNdPtTrackDumpTask::AlidNdPtTrackDumpTask(const char *name) \r
65   : AliAnalysisTaskSE(name)\r
66   , fESD(0)\r
67   , fMC(0)\r
68   , fESDfriend(0)\r
69   , fOutput(0)\r
70   , fPitList(0)\r
71   , fUseMCInfo(kFALSE)\r
72   , fdNdPtEventCuts(0)\r
73   , fdNdPtAcceptanceCuts(0)\r
74   , fdNdPtRecAcceptanceCuts(0)\r
75   , fEsdTrackCuts(0)\r
76   , fTrigger(AliTriggerAnalysis::kMB1) \r
77   , fAnalysisMode(AlidNdPtHelper::kTPC) \r
78   , fOutputSummary(0)\r
79   , fTreeSRedirector(0)\r
80   , fCentralityEstimator(0)\r
81 {\r
82   // Constructor\r
83 \r
84   // Define input and output slots here\r
85   DefineOutput(0, TTree::Class());\r
86   //DefineOutput(1, TList::Class());\r
87 \r
88 }\r
89 \r
90 //_____________________________________________________________________________\r
91 AlidNdPtTrackDumpTask::~AlidNdPtTrackDumpTask()\r
92 {\r
93   if(fOutput) delete fOutput;  fOutput =0; \r
94   //if(fOutputSummary) delete fOutputSummary;  fOutputSummary =0; \r
95   if(fTreeSRedirector) delete fTreeSRedirector;  fTreeSRedirector =0; \r
96 \r
97   if(fdNdPtEventCuts) delete fdNdPtEventCuts; fdNdPtEventCuts=NULL; \r
98   if(fdNdPtAcceptanceCuts) delete fdNdPtAcceptanceCuts; fdNdPtAcceptanceCuts=NULL;\r
99   if(fdNdPtRecAcceptanceCuts) delete fdNdPtRecAcceptanceCuts; fdNdPtRecAcceptanceCuts=NULL;  \r
100   if(fEsdTrackCuts) delete fEsdTrackCuts; fEsdTrackCuts=NULL;\r
101 }\r
102 \r
103 //____________________________________________________________________________\r
104 Bool_t AlidNdPtTrackDumpTask::Notify()\r
105 {\r
106   static Int_t count = 0;\r
107   count++;\r
108   TTree *chain = (TChain*)GetInputData(0);\r
109   if(chain)\r
110   {\r
111     Printf("Processing %d. file: %s", count, chain->GetCurrentFile()->GetName());\r
112   }\r
113 \r
114 return kTRUE;\r
115 }\r
116 \r
117 //_____________________________________________________________________________\r
118 void AlidNdPtTrackDumpTask::UserCreateOutputObjects()\r
119 {\r
120   // Create histograms\r
121   // Called once\r
122   fOutput = new TList;\r
123   fOutput->SetOwner();\r
124 \r
125   //\r
126   // create output tree\r
127   //\r
128   fTreeSRedirector = new TTreeSRedirector("dNdPtOutliersAnalysisPbPb.root");\r
129 \r
130   PostData(0, fOutputSummary);\r
131   //PostData(1, fOutput);\r
132 }\r
133 \r
134 //_____________________________________________________________________________\r
135 void AlidNdPtTrackDumpTask::UserExec(Option_t *) \r
136 {\r
137   //\r
138   // Called for each event\r
139   //\r
140 \r
141   // ESD event\r
142   fESD = (AliESDEvent*) (InputEvent());\r
143   if (!fESD) {\r
144     Printf("ERROR: ESD event not available");\r
145     return;\r
146   }\r
147 \r
148   // MC event\r
149   if(fUseMCInfo) {\r
150     fMC = MCEvent();\r
151     if (!fMC) {\r
152       Printf("ERROR: MC event not available");\r
153       return;\r
154     }\r
155   }\r
156 \r
157   fESDfriend = static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));\r
158   if(!fESDfriend) {\r
159     Printf("ERROR: ESD friends not available");\r
160   }\r
161 \r
162   //\r
163   Process(fESD,fMC,fESDfriend);\r
164 \r
165   // Post output data.\r
166   PostData(0, fOutputSummary);\r
167   //PostData(1, fOutput);\r
168 }\r
169 \r
170 //_____________________________________________________________________________\r
171 void AlidNdPtTrackDumpTask::Process(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)\r
172 {\r
173   //\r
174   // Process real and/or simulated events\r
175   //\r
176   if(!esdEvent) {\r
177     AliDebug(AliLog::kError, "esdEvent not available");\r
178     return;\r
179   }\r
180 \r
181   // get selection cuts\r
182   AlidNdPtEventCuts *evtCuts = GetEventCuts(); \r
183   AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
184   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
185 \r
186   if(!evtCuts || !accCuts  || !esdTrackCuts) {\r
187     AliDebug(AliLog::kError, "cuts not available");\r
188     return;\r
189   }\r
190 \r
191   // trigger selection\r
192   Bool_t isEventTriggered = kTRUE;\r
193   AliPhysicsSelection *physicsSelection = NULL;\r
194   AliTriggerAnalysis* triggerAnalysis = NULL;\r
195 \r
196   // \r
197   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
198   if (!inputHandler)\r
199   {\r
200     Printf("ERROR: Could not receive input handler");\r
201     return;\r
202   }\r
203    \r
204   // get file name\r
205   TTree *chain = (TChain*)GetInputData(0);\r
206   if(!chain) { \r
207     Printf("ERROR: Could not receive input chain");\r
208     return;\r
209   }\r
210   TObjString fileName(chain->GetCurrentFile()->GetName());\r
211 \r
212   // trigger\r
213   if(evtCuts->IsTriggerRequired())  \r
214   {\r
215     // always MB\r
216     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
217 \r
218     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
219     if(!physicsSelection) return;\r
220     //SetPhysicsTriggerSelection(physicsSelection);\r
221 \r
222     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
223       // set trigger (V0AND)\r
224       triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
225       if(!triggerAnalysis) return;\r
226       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
227     }\r
228   }\r
229 \r
230   // centrality determination\r
231   Float_t centralityF = -1;\r
232   AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
233   centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
234 \r
235   // use MC information\r
236   AliHeader* header = 0;\r
237   AliGenEventHeader* genHeader = 0;\r
238   AliStack* stack = 0;\r
239   TArrayF vtxMC(3);\r
240 \r
241   Int_t multMCTrueTracks = 0;\r
242   if(IsUseMCInfo())\r
243   {\r
244     //\r
245     if(!mcEvent) {\r
246       AliDebug(AliLog::kError, "mcEvent not available");\r
247       return;\r
248     }\r
249     // get MC event header\r
250     header = mcEvent->Header();\r
251     if (!header) {\r
252       AliDebug(AliLog::kError, "Header not available");\r
253       return;\r
254     }\r
255     // MC particle stack\r
256     stack = mcEvent->Stack();\r
257     if (!stack) {\r
258       AliDebug(AliLog::kError, "Stack not available");\r
259       return;\r
260     }\r
261 \r
262     // get MC vertex\r
263     genHeader = header->GenEventHeader();\r
264     if (!genHeader) {\r
265       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");\r
266       return;\r
267     }\r
268     genHeader->PrimaryVertex(vtxMC);\r
269 \r
270     // multipliticy of all MC primary tracks\r
271     // in Zv, pt and eta ranges)\r
272     multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);\r
273 \r
274   } // end bUseMC\r
275 \r
276   // get reconstructed vertex  \r
277   //const AliESDVertex* vtxESD = 0; \r
278   const AliESDVertex* vtxESD = 0; \r
279   if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {\r
280         vtxESD = esdEvent->GetPrimaryVertexTPC();\r
281   }\r
282   else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {\r
283      vtxESD = esdEvent->GetPrimaryVertexTracks();\r
284   }\r
285   else {\r
286         return;\r
287   }\r
288 \r
289   if(!vtxESD) return;\r
290 \r
291   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
292   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
293   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
294 \r
295   // check event cuts\r
296   if(isEventOK && isEventTriggered)\r
297   {\r
298     TRandom random;\r
299 \r
300     for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
301     {\r
302       AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
303       if(!track) continue;\r
304       if(track->Charge()==0) continue;\r
305       if(!esdTrackCuts->AcceptTrack(track)) continue;\r
306       if(!accCuts->AcceptTrack(track)) continue;\r
307 \r
308       // downscale low-pT tracks\r
309       if(TMath::Exp(2*track->Pt())<1000*random.Rndm()) continue;\r
310 \r
311       // Dump to the tree \r
312       // vertex\r
313       // TPC constrained tracks\r
314       // InnerParams constrained tracks\r
315       // TPC-ITS tracks\r
316       // ITSout-InnerParams tracks\r
317       // chi2 distance between TPC constrained and TPC-ITS tracks\r
318       // chi2 distance between TPC refitted constrained and TPC-ITS tracks\r
319       // chi2 distance between ITSout and InnerParams tracks\r
320       // MC information\r
321       \r
322       Double_t x[3]; track->GetXYZ(x);\r
323       Double_t b[3]; AliTracker::GetBxByBz(x,b);\r
324       Bool_t isOK = kFALSE;\r
325 \r
326       //\r
327       // Constrain TPC-only tracks (TPCinner) to vertex\r
328       //\r
329       AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());\r
330       if (!tpcInner) continue;\r
331       // transform to the track reference frame \r
332       isOK = tpcInner->Rotate(track->GetAlpha());\r
333       isOK = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
334       if(!isOK) continue;\r
335 \r
336       // clone TPCinner has to be deleted\r
337       AliExternalTrackParam * tpcInnerC = new AliExternalTrackParam(*(track->GetTPCInnerParam()));\r
338       if (!tpcInnerC) continue;\r
339  \r
340       // constrain TPCinner \r
341       //isOK = ConstrainTPCInner(tpcInnerC,vtxESD,esdEvent->GetMagneticField());\r
342       isOK = ConstrainTPCInner(tpcInnerC,vtxESD,b);\r
343 \r
344       // transform to the track reference frame \r
345       isOK = tpcInnerC->Rotate(track->GetAlpha());\r
346       isOK = tpcInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
347 \r
348       if(!isOK) {\r
349         if(tpcInnerC) delete tpcInnerC;\r
350         continue;\r
351       }\r
352 \r
353 \r
354       //\r
355       // Constrain TPC refitted tracks at inner TPC wall (InnerParams) \r
356       // to vertex\r
357       //\r
358       // clone track InnerParams has to be deleted\r
359       AliExternalTrackParam * trackInnerC =  new AliExternalTrackParam(*(track->GetInnerParam()));\r
360       if (!trackInnerC) continue;\r
361  \r
362       // constrain track InnerParams \r
363       isOK = ConstrainTrackInner(trackInnerC,vtxESD,track->GetMass(),b);\r
364 \r
365       // transform to the track reference frame \r
366       isOK = trackInnerC->Rotate(track->GetAlpha());\r
367       isOK = trackInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
368 \r
369       if(!isOK) {\r
370         if(trackInnerC) delete trackInnerC;\r
371         continue;\r
372       }\r
373       \r
374       //\r
375       // calculate chi2 between vi and vj vectors\r
376       // with covi and covj covariance matrices\r
377       // chi2ij = (vi-vj)^(T)*(covi+covj)^(-1)*(vi-vj)\r
378       //\r
379       TMatrixD deltaT(5,1), deltaTtrackC(5,1);\r
380       TMatrixD delta(1,5),  deltatrackC(1,5);\r
381       TMatrixD covarM(5,5), covarMtrackC(5,5);\r
382 \r
383       for (Int_t ipar=0; ipar<5; ipar++) {\r
384         deltaT(ipar,0)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
385         delta(0,ipar)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
386 \r
387         deltaTtrackC(ipar,0)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
388         deltatrackC(0,ipar)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
389 \r
390         for (Int_t jpar=0; jpar<5; jpar++) {\r
391           Int_t index=track->GetIndex(ipar,jpar);\r
392           covarM(ipar,jpar)=track->GetCovariance()[index]+tpcInnerC->GetCovariance()[index];\r
393           covarMtrackC(ipar,jpar)=track->GetCovariance()[index]+trackInnerC->GetCovariance()[index];\r
394         }\r
395       }\r
396       // chi2 distance TPC constrained and TPC+ITS\r
397       TMatrixD covarMInv = covarM.Invert();\r
398       TMatrixD mat2 = covarMInv*deltaT;\r
399       TMatrixD chi2 = delta*mat2; \r
400       //chi2.Print();\r
401 \r
402       // chi2 distance TPC refitted constrained and TPC+ITS\r
403       TMatrixD covarMInvtrackC = covarMtrackC.Invert();\r
404       TMatrixD mat2trackC = covarMInvtrackC*deltaTtrackC;\r
405       TMatrixD chi2trackC = deltatrackC*mat2trackC; \r
406       //chi2trackC.Print();\r
407 \r
408 \r
409       //\r
410       // Propagate ITSout to TPC inner wall \r
411       // and calculate chi2 distance to track (InnerParams)\r
412       //\r
413       const Double_t kTPCRadius=85; \r
414       const Double_t kStep=3; \r
415 \r
416       // clone track InnerParams has to be deleted\r
417       AliExternalTrackParam *trackInnerC2 = new AliExternalTrackParam(*(track->GetInnerParam()));\r
418       if (!trackInnerC2) continue;\r
419       if(!AliTracker::PropagateTrackToBxByBz(trackInnerC2,kTPCRadius,track->GetMass(),kStep,kFALSE))\r
420       {\r
421           if(trackInnerC2) delete trackInnerC2;\r
422           continue;\r
423       }\r
424 \r
425       AliExternalTrackParam *outerITSc = new AliExternalTrackParam();\r
426       if(!outerITSc) continue;\r
427 \r
428       TMatrixD chi2OuterITS(1,1);\r
429       chi2OuterITS(0,0) = 0;\r
430 \r
431       if(esdFriend && esdFriend->TestSkipBit()==kFALSE) \r
432       {\r
433         // propagate ITSout to TPC inner wall\r
434         AliESDfriendTrack *friendTrack = esdFriend->GetTrack(iTrack);\r
435 \r
436         if(friendTrack) \r
437         {\r
438           if( (outerITSc = new AliExternalTrackParam(*friendTrack->GetITSOut())) ) \r
439           {\r
440             if(AliTracker::PropagateTrackToBxByBz(outerITSc,kTPCRadius,track->GetMass(),kStep,kFALSE))\r
441             {\r
442               // transform ITSout to the track InnerParams reference frame \r
443               isOK = kFALSE;\r
444               isOK = outerITSc->Rotate(trackInnerC2->GetAlpha());\r
445               isOK = outerITSc->PropagateTo(trackInnerC2->GetX(),esdEvent->GetMagneticField());\r
446 \r
447               if(!isOK) {\r
448                 if(outerITSc) delete outerITSc;\r
449                 if(trackInnerC2) delete trackInnerC2;\r
450                 continue;\r
451               }\r
452               \r
453               //\r
454               // calculate chi2 between outerITS and innerParams\r
455               // cov without z-coordinate at the moment\r
456               //\r
457               TMatrixD deltaTouterITS(4,1);\r
458               TMatrixD deltaouterITS(1,4);\r
459               TMatrixD covarMouterITS(4,4);\r
460 \r
461               Int_t kipar = 0;\r
462               Int_t kjpar = 0;\r
463               for (Int_t ipar=0; ipar<5; ipar++) {\r
464                 if(ipar!=1) {\r
465                   deltaTouterITS(kipar,0)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];\r
466                   deltaouterITS(0,kipar)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];\r
467                 }\r
468 \r
469                 kjpar=0;\r
470                 for (Int_t jpar=0; jpar<5; jpar++) {\r
471                   Int_t index=outerITSc->GetIndex(ipar,jpar);\r
472                   if(ipar !=1 || jpar!=1) {\r
473                     covarMouterITS(kipar,kjpar)=outerITSc->GetCovariance()[index]+trackInnerC2->GetCovariance()[index];\r
474                   }\r
475                   if(jpar!=1)  kjpar++;\r
476                 }\r
477                 if(ipar!=1) kipar++;\r
478               }\r
479 \r
480               // chi2 distance ITSout and InnerParams\r
481               TMatrixD covarMInvouterITS = covarMouterITS.Invert();\r
482               TMatrixD mat2outerITS = covarMInvouterITS*deltaTouterITS;\r
483               chi2OuterITS = deltaouterITS*mat2outerITS; \r
484               //chi2OuterITS.Print();\r
485             } \r
486           }\r
487         }\r
488       }\r
489 \r
490       //\r
491       // MC info\r
492       //\r
493       TParticle *particle=NULL, *particleTPC=NULL, *particleITS=NULL;\r
494       TParticle *particleMother=NULL, *particleMotherTPC=NULL, *particleMotherITS=NULL;\r
495       Int_t mech=-1, mechTPC=-1, mechITS=-1;\r
496       Bool_t isPrim=kFALSE, isPrimTPC=kFALSE, isPrimITS=kFALSE;\r
497       Bool_t isFromStrangess=kFALSE, isFromStrangessTPC=kFALSE, isFromStrangessITS=kFALSE;\r
498       Bool_t isFromConversion=kFALSE, isFromConversionTPC=kFALSE, isFromConversionITS=kFALSE;\r
499       Bool_t isFromMaterial=kFALSE, isFromMaterialTPC=kFALSE, isFromMaterialITS=kFALSE;\r
500 \r
501       AliTrackReference *refTPCIn = NULL;\r
502       AliTrackReference *refITS = NULL;\r
503       AliExternalTrackParam *trackInnerC3 = new AliExternalTrackParam(*(track->GetInnerParam()));\r
504       if(!trackInnerC3) continue;\r
505 \r
506       if(IsUseMCInfo()) \r
507       {\r
508         if(!stack) return;\r
509 \r
510         //\r
511         // global track\r
512         //\r
513         Int_t label = TMath::Abs(track->GetLabel()); \r
514         particle = stack->Particle(label);\r
515         if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0.)\r
516         {\r
517           particleMother = GetMother(particle,stack);\r
518           mech = particle->GetUniqueID();\r
519           isPrim = stack->IsPhysicalPrimary(label);\r
520           isFromStrangess  = IsFromStrangeness(label,stack);\r
521           isFromConversion = IsFromConversion(label,stack);\r
522           isFromMaterial   = IsFromMaterial(label,stack);\r
523         }\r
524 \r
525         //\r
526         // TPC track\r
527         //\r
528         Int_t labelTPC = TMath::Abs(track->GetTPCLabel()); \r
529         particleTPC = stack->Particle(labelTPC);\r
530         if(particleTPC && particleTPC->GetPDG() && particleTPC->GetPDG()->Charge()!=0.)\r
531         {\r
532           particleMotherTPC = GetMother(particleTPC,stack);\r
533           mechTPC = particleTPC->GetUniqueID();\r
534           isPrimTPC = stack->IsPhysicalPrimary(labelTPC);\r
535           isFromStrangessTPC  = IsFromStrangeness(labelTPC,stack);\r
536           isFromConversionTPC = IsFromConversion(labelTPC,stack);\r
537           isFromMaterialTPC   = IsFromMaterial(labelTPC,stack);\r
538         }\r
539 \r
540         //\r
541         // store first track reference\r
542         // for TPC track\r
543         //\r
544         TParticle *part=0;\r
545         TClonesArray *trefs=0;\r
546         Int_t status = mcEvent->GetParticleAndTR(track->GetTPCLabel(), part, trefs);\r
547 \r
548         if(status>0 && part && trefs && part->GetPDG() && part->GetPDG()->Charge()!=0.) \r
549         {\r
550           Int_t nTrackRef = trefs->GetEntries();\r
551           //printf("nTrackRef %d \n",nTrackRef);\r
552 \r
553           Int_t countITS = 0;\r
554           for (Int_t iref = 0; iref < nTrackRef; iref++) \r
555           {\r
556             AliTrackReference *ref = (AliTrackReference *)trefs->At(iref);\r
557             //printf("ref %p \n",ref);\r
558             //if(ref) printf("ref->DetectorId() %d \n",ref->DetectorId());\r
559             //printf("AliTrackReference::kTPC  %d \n",AliTrackReference::kTPC);\r
560 \r
561 \r
562              // Ref. in the middle ITS \r
563             if(ref && ref->DetectorId()==AliTrackReference::kITS)\r
564             {\r
565               if(!refITS && countITS==2) {\r
566                 refITS = ref;\r
567                 //printf("refITS %p \n",refITS);\r
568               }\r
569               countITS++;\r
570             }\r
571 \r
572             // TPC\r
573             if(ref && ref->DetectorId()==AliTrackReference::kTPC)\r
574             {\r
575               if(!refTPCIn) {\r
576                 refTPCIn = ref;\r
577                 //printf("refTPCIn %p \n",refTPCIn);\r
578                 //break;\r
579               }\r
580             }\r
581           }\r
582 \r
583           // transform inner params to TrackRef\r
584           // reference frame\r
585           if(refTPCIn && trackInnerC3) \r
586           {\r
587             Double_t kRefPhi = TMath::ATan2(refTPCIn->Y(),refTPCIn->X());\r
588             isOK = kFALSE;\r
589             isOK = trackInnerC3->Rotate(kRefPhi);\r
590             isOK = AliTracker::PropagateTrackToBxByBz(trackInnerC3,refTPCIn->R(),track->GetMass(),kStep,kFALSE);\r
591 \r
592             if(!isOK){\r
593               if(trackInnerC3) delete trackInnerC3;\r
594               if(refTPCIn) delete refTPCIn;\r
595             }\r
596           }\r
597         }\r
598 \r
599         //\r
600         // ITS track\r
601         //\r
602         Int_t labelITS = TMath::Abs(track->GetITSLabel()); \r
603         particleITS = stack->Particle(labelITS);\r
604         if(particleITS && particleITS->GetPDG() && particleITS->GetPDG()->Charge()!=0.)\r
605         {\r
606           particleMotherITS = GetMother(particleITS,stack);\r
607           mechITS = particleITS->GetUniqueID();\r
608           isPrimITS = stack->IsPhysicalPrimary(labelITS);\r
609           isFromStrangessITS  = IsFromStrangeness(labelITS,stack);\r
610           isFromConversionITS = IsFromConversion(labelITS,stack);\r
611           isFromMaterialITS   = IsFromMaterial(labelITS,stack);\r
612         }\r
613       }\r
614 \r
615       //\r
616       Double_t vert[3] = {0}; \r
617       vert[0] = vtxESD->GetXv();\r
618       vert[1] = vtxESD->GetYv();\r
619       vert[2] = vtxESD->GetZv();\r
620       Int_t mult = vtxESD->GetNContributors();\r
621       Double_t bz = esdEvent->GetMagneticField();\r
622       Double_t runNumber = esdEvent->GetRunNumber();\r
623       Double_t evtTimeStamp = esdEvent->GetTimeStamp();\r
624       Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
625 \r
626 \r
627       //\r
628       if(!fTreeSRedirector) return;\r
629       (*fTreeSRedirector)<<"dNdPtTree"<<\r
630         "fileName.="<<&fileName<<\r
631         "runNumber="<<runNumber<<\r
632         "evtTimeStamp="<<evtTimeStamp<<\r
633         "evtNumberInFile="<<evtNumberInFile<<\r
634         "Bz="<<bz<<\r
635         "vertX="<<vert[0]<<\r
636         "vertY="<<vert[1]<<\r
637         "vertZ="<<vert[2]<<\r
638         "mult="<<mult<<\r
639         "esdTrack.="<<track<<\r
640         "extTPCInnerC.="<<tpcInnerC<<\r
641         "extInnerParamC.="<<trackInnerC<<\r
642         "extInnerParam.="<<trackInnerC2<<\r
643         "extOuterITS.="<<outerITSc<<\r
644         "extInnerParamRef.="<<trackInnerC3<<\r
645         "refTPCIn.="<<refTPCIn<<\r
646         "refITS.="<<refITS<<\r
647         "chi2TPCInnerC="<<chi2(0,0)<<\r
648         "chi2InnerC="<<chi2trackC(0,0)<<\r
649         "chi2OuterITS="<<chi2OuterITS(0,0)<<\r
650         "centralityF="<<centralityF<<\r
651         "particle.="<<particle<<\r
652         "particleMother.="<<particleMother<<\r
653         "mech="<<mech<<\r
654         "isPrim="<<isPrim<<\r
655         "isFromStrangess="<<isFromStrangess<<\r
656         "isFromConversion="<<isFromConversion<<\r
657         "isFromMaterial="<<isFromMaterial<<\r
658         "particleTPC.="<<particleTPC<<\r
659         "particleMotherTPC.="<<particleMotherTPC<<\r
660         "mechTPC="<<mechTPC<<\r
661         "isPrimTPC="<<isPrimTPC<<\r
662         "isFromStrangessTPC="<<isFromStrangessTPC<<\r
663         "isFromConversionTPC="<<isFromConversionTPC<<\r
664         "isFromMaterialTPC="<<isFromMaterialTPC<<\r
665         "particleITS.="<<particleITS<<\r
666         "particleMotherITS.="<<particleMotherITS<<\r
667         "mechITS="<<mechITS<<\r
668         "isPrimITS="<<isPrimITS<<\r
669         "isFromStrangessITS="<<isFromStrangessITS<<\r
670         "isFromConversionITS="<<isFromConversionITS<<\r
671         "isFromMaterialITS="<<isFromMaterialITS<<\r
672         "\n";\r
673 \r
674         if(tpcInnerC) delete tpcInnerC;\r
675         if(trackInnerC) delete trackInnerC;\r
676         if(trackInnerC2) delete trackInnerC2;\r
677         if(outerITSc) delete outerITSc;\r
678 \r
679         if(trackInnerC3) delete trackInnerC3;\r
680     }\r
681   }\r
682 \r
683   PostData(0, fOutputSummary);\r
684   //PostData(1, fOutput);\r
685 }\r
686 \r
687 \r
688 //_____________________________________________________________________________\r
689 Bool_t AlidNdPtTrackDumpTask::ConstrainTPCInner(AliExternalTrackParam *const tpcInnerC, const AliESDVertex* vtx, Double_t b[3])\r
690 {\r
691  // Constrain TPC inner params constrained\r
692  //\r
693       if(!tpcInnerC) return kFALSE; \r
694       if(!vtx) return kFALSE;\r
695 \r
696       Double_t dz[2],cov[3];\r
697       //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();\r
698       //if(!tpcInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; \r
699       //if(!tpcInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; \r
700       if(!tpcInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; \r
701 \r
702 \r
703       Double_t covar[6]; vtx->GetCovMatrix(covar);\r
704       Double_t p[2]={tpcInnerC->GetParameter()[0]-dz[0],tpcInnerC->GetParameter()[1]-dz[1]};\r
705       Double_t c[3]={covar[2],0.,covar[5]};\r
706       Double_t chi2C=tpcInnerC->GetPredictedChi2(p,c);\r
707       if (chi2C>kVeryBig) return kFALSE; \r
708 \r
709       if(!tpcInnerC->Update(p,c)) return kFALSE;\r
710 \r
711   return kTRUE;\r
712 }\r
713 \r
714 //_____________________________________________________________________________\r
715 Bool_t AlidNdPtTrackDumpTask::ConstrainTrackInner(AliExternalTrackParam *const trackInnerC, const AliESDVertex* vtx, Double_t mass, Double_t b[3])\r
716 {\r
717  // Constrain TPC inner params constrained\r
718  //\r
719       if(!trackInnerC) return kFALSE; \r
720       if(!vtx) return kFALSE;\r
721 \r
722       const Double_t kRadius  = 2.8; \r
723       const Double_t kMaxStep = 1.0; \r
724 \r
725       Double_t dz[2],cov[3];\r
726 \r
727       //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();\r
728       //if(!trackInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; \r
729       //if(!trackInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; \r
730 \r
731       if(!AliTracker::PropagateTrackToBxByBz(trackInnerC,kRadius,mass,kMaxStep,kFALSE)) return kFALSE;\r
732       if(!trackInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; \r
733 \r
734       //\r
735       Double_t covar[6]; vtx->GetCovMatrix(covar);\r
736       Double_t p[2]={trackInnerC->GetParameter()[0]-dz[0],trackInnerC->GetParameter()[1]-dz[1]};\r
737       Double_t c[3]={covar[2],0.,covar[5]};\r
738       Double_t chi2C=trackInnerC->GetPredictedChi2(p,c);\r
739       if (chi2C>kVeryBig) return kFALSE; \r
740 \r
741       if(!trackInnerC->Update(p,c)) return kFALSE;\r
742 \r
743   return kTRUE;\r
744 }\r
745 \r
746 \r
747 //_____________________________________________________________________________\r
748 TParticle *AlidNdPtTrackDumpTask::GetMother(TParticle *const particle, AliStack *const stack) \r
749 {\r
750   if(!particle) return NULL;\r
751   if(!stack) return NULL;\r
752 \r
753   Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
754   TParticle* mother = NULL; \r
755   mother = stack->Particle(motherLabel); \r
756 \r
757 return mother;\r
758 }\r
759 \r
760 //_____________________________________________________________________________\r
761 Bool_t AlidNdPtTrackDumpTask::IsFromConversion(const Int_t label, AliStack *const stack) \r
762 {\r
763   Bool_t isFromConversion = kFALSE;\r
764 \r
765   if(stack) {\r
766     TParticle* particle = stack->Particle(label);\r
767 \r
768     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
769     {\r
770        Int_t mech = particle->GetUniqueID(); // production mechanism \r
771        Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
772 \r
773        Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
774        TParticle* mother = stack->Particle(motherLabel); \r
775        if(mother) {\r
776           Int_t motherPdg = mother->GetPdgCode();\r
777 \r
778           if(!isPrim && mech==5 && motherPdg==kGamma) { \r
779             isFromConversion=kTRUE; \r
780           }\r
781        }\r
782     } \r
783   } \r
784 \r
785   return isFromConversion;\r
786 }\r
787 \r
788 //_____________________________________________________________________________\r
789 Bool_t AlidNdPtTrackDumpTask::IsFromMaterial(const Int_t label, AliStack *const stack) \r
790 {\r
791   Bool_t isFromMaterial = kFALSE;\r
792 \r
793   if(stack) {\r
794     TParticle* particle = stack->Particle(label);\r
795 \r
796     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
797     {\r
798        Int_t mech = particle->GetUniqueID(); // production mechanism \r
799        Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
800 \r
801        Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
802        TParticle* mother = stack->Particle(motherLabel); \r
803        if(mother) {\r
804           if(!isPrim && mech==13) { \r
805             isFromMaterial=kTRUE; \r
806           }\r
807        }\r
808      } \r
809   } \r
810 \r
811   return isFromMaterial;\r
812 }\r
813 \r
814 //_____________________________________________________________________________\r
815 Bool_t AlidNdPtTrackDumpTask::IsFromStrangeness(const Int_t label, AliStack *const stack) \r
816 {\r
817   Bool_t isFromStrangeness = kFALSE;\r
818 \r
819   if(stack) {\r
820     TParticle* particle = stack->Particle(label);\r
821 \r
822     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
823     {\r
824        Int_t mech = particle->GetUniqueID(); // production mechanism \r
825        Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
826 \r
827        Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
828        TParticle* mother = stack->Particle(motherLabel); \r
829        if(mother) {\r
830           Int_t motherPdg = mother->GetPdgCode();\r
831 \r
832           // K+-, lambda, antilambda, K0s decays\r
833           if(!isPrim && mech==4 && \r
834               (TMath::Abs(motherPdg)==kKPlus || TMath::Abs(motherPdg)==kLambda0 || motherPdg==kK0Short))\r
835           {\r
836             isFromStrangeness = kTRUE;\r
837           } \r
838        }\r
839     } \r
840   } \r
841 \r
842   return isFromStrangeness;\r
843 }\r
844 \r
845 \r
846 //_____________________________________________________________________________\r
847 void AlidNdPtTrackDumpTask::FinishTaskOutput() \r
848 {\r
849   //\r
850   // Called one at the end \r
851   // locally on working node\r
852   //\r
853 \r
854   // Post output data.\r
855   PostData(1, fOutput);\r
856   //PostData(0, fOutputSummary);\r
857 }\r
858 \r
859 //_____________________________________________________________________________\r
860 void AlidNdPtTrackDumpTask::Terminate(Option_t *) \r
861 {\r
862   // Called one at the end \r
863   \r
864   // check output data\r
865   fOutputSummary = dynamic_cast<TTree*> (GetOutputData(0));\r
866   if(fOutputSummary) delete fOutputSummary; fOutputSummary=0;\r
867   if(fTreeSRedirector)  delete fTreeSRedirector; fTreeSRedirector=0;\r
868 \r
869   TChain* chain = new TChain("dNdPtTree");\r
870   if(!chain) return;\r
871   chain->Add("dNdPtOutliersAnalysisPbPb.root");\r
872   TTree *tree = chain->CopyTree("1");\r
873   if (chain) { delete chain; chain=0; }\r
874   if(!tree) return;\r
875   tree->Print();\r
876   fOutputSummary = tree;\r
877 \r
878   if (!fOutputSummary) {\r
879     Printf("ERROR: AlidNdPtTrackDumpTask::Terminate(): Output data not avaiable GetOutputData(0)==0x0 ..." );\r
880     return;\r
881   }\r
882   \r
883 \r
884 \r
885   PostData(0, fOutputSummary);\r
886   //PostData(1, fOutput);\r
887 }\r