- V0s tree 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 #include <TDatabasePDG.h>\r
20 \r
21 #include "TChain.h"\r
22 #include "TTreeStream.h"\r
23 #include "TTree.h"\r
24 #include "TH1F.h"\r
25 #include "TCanvas.h"\r
26 #include "TList.h"\r
27 #include "TFile.h"\r
28 #include "TMatrixD.h"\r
29 #include "TRandom3.h"\r
30 \r
31 #include "AliHeader.h"  \r
32 #include "AliGenEventHeader.h"  \r
33 #include "AliInputEventHandler.h"  \r
34 #include "AliStack.h"  \r
35 #include "AliTrackReference.h"  \r
36 \r
37 #include "AliPhysicsSelection.h"\r
38 #include "AliAnalysisTask.h"\r
39 #include "AliAnalysisManager.h"\r
40 #include "AliESDEvent.h"\r
41 #include "AliESDfriend.h"\r
42 #include "AliMCEvent.h"\r
43 #include "AliESDInputHandler.h"\r
44 #include "AliESDVertex.h"\r
45 #include "AliTracker.h"\r
46 #include "AliGeomManager.h"\r
47 \r
48 #include "AliCentrality.h"\r
49 #include "AliESDVZERO.h"\r
50 #include "AliMultiplicity.h"\r
51 \r
52 #include "AliESDtrackCuts.h"\r
53 #include "AliMCEventHandler.h"\r
54 #include "AlidNdPt.h"\r
55 #include "AlidNdPtEventCuts.h"\r
56 #include "AlidNdPtAcceptanceCuts.h"\r
57 \r
58 #include "AlidNdPtTrackDumpTask.h"\r
59 #include "AliKFParticle.h"\r
60 #include "AliESDv0.h"\r
61 \r
62 using namespace std;\r
63 \r
64 ClassImp(AlidNdPtTrackDumpTask)\r
65 \r
66 //_____________________________________________________________________________\r
67 AlidNdPtTrackDumpTask::AlidNdPtTrackDumpTask(const char *name) \r
68   : AliAnalysisTaskSE(name)\r
69   , fESD(0)\r
70   , fMC(0)\r
71   , fESDfriend(0)\r
72   , fOutput(0)\r
73   , fPitList(0)\r
74   , fUseMCInfo(kFALSE)\r
75   , fdNdPtEventCuts(0)\r
76   , fdNdPtAcceptanceCuts(0)\r
77   , fdNdPtRecAcceptanceCuts(0)\r
78   , fEsdTrackCuts(0)\r
79   , fTrigger(AliTriggerAnalysis::kMB1) \r
80   , fAnalysisMode(AlidNdPtHelper::kTPC) \r
81   , fOutputSummary(0)\r
82   , fTreeSRedirector(0)\r
83   , fCentralityEstimator(0)\r
84   , fLowPtTrackDownscaligF(0)\r
85   , fLowPtV0DownscaligF(0)\r
86 {\r
87   // Constructor\r
88 \r
89   // Define input and output slots here\r
90   DefineOutput(0, TTree::Class());\r
91   //DefineOutput(1, TList::Class());\r
92 \r
93 }\r
94 \r
95 //_____________________________________________________________________________\r
96 AlidNdPtTrackDumpTask::~AlidNdPtTrackDumpTask()\r
97 {\r
98   if(fOutput) delete fOutput;  fOutput =0; \r
99   //if(fOutputSummary) delete fOutputSummary;  fOutputSummary =0; \r
100   if(fTreeSRedirector) delete fTreeSRedirector;  fTreeSRedirector =0; \r
101 \r
102   if(fdNdPtEventCuts) delete fdNdPtEventCuts; fdNdPtEventCuts=NULL; \r
103   if(fdNdPtAcceptanceCuts) delete fdNdPtAcceptanceCuts; fdNdPtAcceptanceCuts=NULL;\r
104   if(fdNdPtRecAcceptanceCuts) delete fdNdPtRecAcceptanceCuts; fdNdPtRecAcceptanceCuts=NULL;  \r
105   if(fEsdTrackCuts) delete fEsdTrackCuts; fEsdTrackCuts=NULL;\r
106 }\r
107 \r
108 //____________________________________________________________________________\r
109 Bool_t AlidNdPtTrackDumpTask::Notify()\r
110 {\r
111   static Int_t count = 0;\r
112   count++;\r
113   TTree *chain = (TChain*)GetInputData(0);\r
114   if(chain)\r
115   {\r
116     Printf("Processing %d. file: %s", count, chain->GetCurrentFile()->GetName());\r
117   }\r
118 \r
119 return kTRUE;\r
120 }\r
121 \r
122 //_____________________________________________________________________________\r
123 void AlidNdPtTrackDumpTask::UserCreateOutputObjects()\r
124 {\r
125   // Create histograms\r
126   // Called once\r
127   fOutput = new TList;\r
128   fOutput->SetOwner();\r
129 \r
130   //\r
131   // create output tree\r
132   //\r
133   fTreeSRedirector = new TTreeSRedirector("jotwinow_HighPt_TrackAndV0_Trees.root");\r
134 \r
135   PostData(0, fOutputSummary);\r
136   //PostData(1, fOutput);\r
137 }\r
138 \r
139 //_____________________________________________________________________________\r
140 void AlidNdPtTrackDumpTask::UserExec(Option_t *) \r
141 {\r
142   //\r
143   // Called for each event\r
144   //\r
145 \r
146   // ESD event\r
147   fESD = (AliESDEvent*) (InputEvent());\r
148   if (!fESD) {\r
149     Printf("ERROR: ESD event not available");\r
150     return;\r
151   }\r
152 \r
153   // MC event\r
154   if(fUseMCInfo) {\r
155     fMC = MCEvent();\r
156     if (!fMC) {\r
157       Printf("ERROR: MC event not available");\r
158       return;\r
159     }\r
160   }\r
161 \r
162   fESDfriend = static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));\r
163   if(!fESDfriend) {\r
164     Printf("ERROR: ESD friends not available");\r
165   }\r
166 \r
167   //\r
168   Process(fESD,fMC,fESDfriend);\r
169   ProcessV0(fESD,fMC,fESDfriend);\r
170 \r
171   // Post output data.\r
172   PostData(0, fOutputSummary);\r
173   //PostData(1, fOutput);\r
174 }\r
175 \r
176 //_____________________________________________________________________________\r
177 void AlidNdPtTrackDumpTask::Process(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)\r
178 {\r
179   //\r
180   // Process real and/or simulated events\r
181   //\r
182   if(!esdEvent) {\r
183     AliDebug(AliLog::kError, "esdEvent not available");\r
184     return;\r
185   }\r
186 \r
187   // get selection cuts\r
188   AlidNdPtEventCuts *evtCuts = GetEventCuts(); \r
189   AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
190   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
191 \r
192   if(!evtCuts || !accCuts  || !esdTrackCuts) {\r
193     AliDebug(AliLog::kError, "cuts not available");\r
194     return;\r
195   }\r
196 \r
197   // trigger selection\r
198   Bool_t isEventTriggered = kTRUE;\r
199   AliPhysicsSelection *physicsSelection = NULL;\r
200   AliTriggerAnalysis* triggerAnalysis = NULL;\r
201 \r
202   // \r
203   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
204   if (!inputHandler)\r
205   {\r
206     Printf("ERROR: Could not receive input handler");\r
207     return;\r
208   }\r
209    \r
210   // get file name\r
211   TTree *chain = (TChain*)GetInputData(0);\r
212   if(!chain) { \r
213     Printf("ERROR: Could not receive input chain");\r
214     return;\r
215   }\r
216   TObjString fileName(chain->GetCurrentFile()->GetName());\r
217 \r
218   // trigger\r
219   if(evtCuts->IsTriggerRequired())  \r
220   {\r
221     // always MB\r
222     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
223 \r
224     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
225     if(!physicsSelection) return;\r
226     //SetPhysicsTriggerSelection(physicsSelection);\r
227 \r
228     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
229       // set trigger (V0AND)\r
230       triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
231       if(!triggerAnalysis) return;\r
232       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
233     }\r
234   }\r
235 \r
236   // centrality determination\r
237   Float_t centralityF = -1;\r
238   AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
239   centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
240 \r
241   // use MC information\r
242   AliHeader* header = 0;\r
243   AliGenEventHeader* genHeader = 0;\r
244   AliStack* stack = 0;\r
245   TArrayF vtxMC(3);\r
246 \r
247   Int_t multMCTrueTracks = 0;\r
248   if(IsUseMCInfo())\r
249   {\r
250     //\r
251     if(!mcEvent) {\r
252       AliDebug(AliLog::kError, "mcEvent not available");\r
253       return;\r
254     }\r
255     // get MC event header\r
256     header = mcEvent->Header();\r
257     if (!header) {\r
258       AliDebug(AliLog::kError, "Header not available");\r
259       return;\r
260     }\r
261     // MC particle stack\r
262     stack = mcEvent->Stack();\r
263     if (!stack) {\r
264       AliDebug(AliLog::kError, "Stack not available");\r
265       return;\r
266     }\r
267 \r
268     // get MC vertex\r
269     genHeader = header->GenEventHeader();\r
270     if (!genHeader) {\r
271       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");\r
272       return;\r
273     }\r
274     genHeader->PrimaryVertex(vtxMC);\r
275 \r
276     // multipliticy of all MC primary tracks\r
277     // in Zv, pt and eta ranges)\r
278     multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);\r
279 \r
280   } // end bUseMC\r
281 \r
282   // get reconstructed vertex  \r
283   //const AliESDVertex* vtxESD = 0; \r
284   const AliESDVertex* vtxESD = 0; \r
285   if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {\r
286         vtxESD = esdEvent->GetPrimaryVertexTPC();\r
287   }\r
288   else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {\r
289      vtxESD = esdEvent->GetPrimaryVertexTracks();\r
290   }\r
291   else {\r
292         return;\r
293   }\r
294 \r
295   if(!vtxESD) return;\r
296 \r
297   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
298   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
299   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
300 \r
301   // check event cuts\r
302   if(isEventOK && isEventTriggered)\r
303   {\r
304     TRandom3 random;\r
305 \r
306     for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
307     {\r
308       AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
309       if(!track) continue;\r
310       if(track->Charge()==0) continue;\r
311       if(!esdTrackCuts->AcceptTrack(track)) continue;\r
312       if(!accCuts->AcceptTrack(track)) continue;\r
313 \r
314       // downscale low-pT tracks\r
315       Double_t scalempt= TMath::Min(track->Pt(),10.);\r
316       if(TMath::Exp(2*scalempt)<fLowPtTrackDownscaligF*random.Rndm()) continue;\r
317 \r
318       // Dump to the tree \r
319       // vertex\r
320       // TPC constrained tracks\r
321       // InnerParams constrained tracks\r
322       // TPC-ITS tracks\r
323       // ITSout-InnerParams tracks\r
324       // chi2 distance between TPC constrained and TPC-ITS tracks\r
325       // chi2 distance between TPC refitted constrained and TPC-ITS tracks\r
326       // chi2 distance between ITSout and InnerParams tracks\r
327       // MC information\r
328       \r
329       Double_t x[3]; track->GetXYZ(x);\r
330       Double_t b[3]; AliTracker::GetBxByBz(x,b);\r
331       Bool_t isOK = kFALSE;\r
332 \r
333       //\r
334       // Constrain TPC-only tracks (TPCinner) to vertex\r
335       //\r
336       AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());\r
337       if (!tpcInner) continue;\r
338       // transform to the track reference frame \r
339       isOK = tpcInner->Rotate(track->GetAlpha());\r
340       isOK = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
341       if(!isOK) continue;\r
342 \r
343       // clone TPCinner has to be deleted\r
344       AliExternalTrackParam * tpcInnerC = new AliExternalTrackParam(*(track->GetTPCInnerParam()));\r
345       if (!tpcInnerC) continue;\r
346  \r
347       // constrain TPCinner \r
348       //isOK = ConstrainTPCInner(tpcInnerC,vtxESD,esdEvent->GetMagneticField());\r
349       isOK = ConstrainTPCInner(tpcInnerC,vtxESD,b);\r
350 \r
351       // transform to the track reference frame \r
352       isOK = tpcInnerC->Rotate(track->GetAlpha());\r
353       isOK = tpcInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
354 \r
355       if(!isOK) {\r
356         if(tpcInnerC) delete tpcInnerC;\r
357         continue;\r
358       }\r
359 \r
360 \r
361       //\r
362       // Constrain TPC refitted tracks at inner TPC wall (InnerParams) \r
363       // to vertex\r
364       //\r
365       // clone track InnerParams has to be deleted\r
366       AliExternalTrackParam * trackInnerC =  new AliExternalTrackParam(*(track->GetInnerParam()));\r
367       if (!trackInnerC) continue;\r
368  \r
369       // constrain track InnerParams \r
370       isOK = ConstrainTrackInner(trackInnerC,vtxESD,track->GetMass(),b);\r
371 \r
372       // transform to the track reference frame \r
373       isOK = trackInnerC->Rotate(track->GetAlpha());\r
374       isOK = trackInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
375 \r
376       if(!isOK) {\r
377         if(trackInnerC) delete trackInnerC;\r
378         continue;\r
379       }\r
380       \r
381       //\r
382       // calculate chi2 between vi and vj vectors\r
383       // with covi and covj covariance matrices\r
384       // chi2ij = (vi-vj)^(T)*(covi+covj)^(-1)*(vi-vj)\r
385       //\r
386       TMatrixD deltaT(5,1), deltaTtrackC(5,1);\r
387       TMatrixD delta(1,5),  deltatrackC(1,5);\r
388       TMatrixD covarM(5,5), covarMtrackC(5,5);\r
389 \r
390       for (Int_t ipar=0; ipar<5; ipar++) {\r
391         deltaT(ipar,0)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
392         delta(0,ipar)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
393 \r
394         deltaTtrackC(ipar,0)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
395         deltatrackC(0,ipar)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
396 \r
397         for (Int_t jpar=0; jpar<5; jpar++) {\r
398           Int_t index=track->GetIndex(ipar,jpar);\r
399           covarM(ipar,jpar)=track->GetCovariance()[index]+tpcInnerC->GetCovariance()[index];\r
400           covarMtrackC(ipar,jpar)=track->GetCovariance()[index]+trackInnerC->GetCovariance()[index];\r
401         }\r
402       }\r
403       // chi2 distance TPC constrained and TPC+ITS\r
404       TMatrixD covarMInv = covarM.Invert();\r
405       TMatrixD mat2 = covarMInv*deltaT;\r
406       TMatrixD chi2 = delta*mat2; \r
407       //chi2.Print();\r
408 \r
409       // chi2 distance TPC refitted constrained and TPC+ITS\r
410       TMatrixD covarMInvtrackC = covarMtrackC.Invert();\r
411       TMatrixD mat2trackC = covarMInvtrackC*deltaTtrackC;\r
412       TMatrixD chi2trackC = deltatrackC*mat2trackC; \r
413       //chi2trackC.Print();\r
414 \r
415 \r
416       //\r
417       // Propagate ITSout to TPC inner wall \r
418       // and calculate chi2 distance to track (InnerParams)\r
419       //\r
420       const Double_t kTPCRadius=85; \r
421       const Double_t kStep=3; \r
422 \r
423       // clone track InnerParams has to be deleted\r
424       AliExternalTrackParam *trackInnerC2 = new AliExternalTrackParam(*(track->GetInnerParam()));\r
425       if (!trackInnerC2) continue;\r
426       if(!AliTracker::PropagateTrackToBxByBz(trackInnerC2,kTPCRadius,track->GetMass(),kStep,kFALSE))\r
427       {\r
428           if(trackInnerC2) delete trackInnerC2;\r
429           continue;\r
430       }\r
431 \r
432       AliExternalTrackParam *outerITSc = new AliExternalTrackParam();\r
433       if(!outerITSc) continue;\r
434 \r
435       TMatrixD chi2OuterITS(1,1);\r
436       chi2OuterITS(0,0) = 0;\r
437 \r
438       if(esdFriend && esdFriend->TestSkipBit()==kFALSE) \r
439       {\r
440         // propagate ITSout to TPC inner wall\r
441         AliESDfriendTrack *friendTrack = esdFriend->GetTrack(iTrack);\r
442 \r
443         if(friendTrack) \r
444         {\r
445           if( (outerITSc = new AliExternalTrackParam(*friendTrack->GetITSOut())) ) \r
446           {\r
447             if(AliTracker::PropagateTrackToBxByBz(outerITSc,kTPCRadius,track->GetMass(),kStep,kFALSE))\r
448             {\r
449               // transform ITSout to the track InnerParams reference frame \r
450               isOK = kFALSE;\r
451               isOK = outerITSc->Rotate(trackInnerC2->GetAlpha());\r
452               isOK = outerITSc->PropagateTo(trackInnerC2->GetX(),esdEvent->GetMagneticField());\r
453 \r
454               if(!isOK) {\r
455                 if(outerITSc) delete outerITSc;\r
456                 if(trackInnerC2) delete trackInnerC2;\r
457                 continue;\r
458               }\r
459               \r
460               //\r
461               // calculate chi2 between outerITS and innerParams\r
462               // cov without z-coordinate at the moment\r
463               //\r
464               TMatrixD deltaTouterITS(4,1);\r
465               TMatrixD deltaouterITS(1,4);\r
466               TMatrixD covarMouterITS(4,4);\r
467 \r
468               Int_t kipar = 0;\r
469               Int_t kjpar = 0;\r
470               for (Int_t ipar=0; ipar<5; ipar++) {\r
471                 if(ipar!=1) {\r
472                   deltaTouterITS(kipar,0)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];\r
473                   deltaouterITS(0,kipar)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];\r
474                 }\r
475 \r
476                 kjpar=0;\r
477                 for (Int_t jpar=0; jpar<5; jpar++) {\r
478                   Int_t index=outerITSc->GetIndex(ipar,jpar);\r
479                   if(ipar !=1 || jpar!=1) {\r
480                     covarMouterITS(kipar,kjpar)=outerITSc->GetCovariance()[index]+trackInnerC2->GetCovariance()[index];\r
481                   }\r
482                   if(jpar!=1)  kjpar++;\r
483                 }\r
484                 if(ipar!=1) kipar++;\r
485               }\r
486 \r
487               // chi2 distance ITSout and InnerParams\r
488               TMatrixD covarMInvouterITS = covarMouterITS.Invert();\r
489               TMatrixD mat2outerITS = covarMInvouterITS*deltaTouterITS;\r
490               chi2OuterITS = deltaouterITS*mat2outerITS; \r
491               //chi2OuterITS.Print();\r
492             } \r
493           }\r
494         }\r
495       }\r
496 \r
497       //\r
498       // MC info\r
499       //\r
500       TParticle *particle=NULL, *particleTPC=NULL, *particleITS=NULL;\r
501       TParticle *particleMother=NULL, *particleMotherTPC=NULL, *particleMotherITS=NULL;\r
502       Int_t mech=-1, mechTPC=-1, mechITS=-1;\r
503       Bool_t isPrim=kFALSE, isPrimTPC=kFALSE, isPrimITS=kFALSE;\r
504       Bool_t isFromStrangess=kFALSE, isFromStrangessTPC=kFALSE, isFromStrangessITS=kFALSE;\r
505       Bool_t isFromConversion=kFALSE, isFromConversionTPC=kFALSE, isFromConversionITS=kFALSE;\r
506       Bool_t isFromMaterial=kFALSE, isFromMaterialTPC=kFALSE, isFromMaterialITS=kFALSE;\r
507 \r
508       AliTrackReference *refTPCIn = NULL;\r
509       AliTrackReference *refITS = NULL;\r
510       AliExternalTrackParam *trackInnerC3 = new AliExternalTrackParam(*(track->GetInnerParam()));\r
511       if(!trackInnerC3) continue;\r
512 \r
513       if(IsUseMCInfo()) \r
514       {\r
515         if(!stack) return;\r
516 \r
517         //\r
518         // global track\r
519         //\r
520         Int_t label = TMath::Abs(track->GetLabel()); \r
521         particle = stack->Particle(label);\r
522         if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0.)\r
523         {\r
524           particleMother = GetMother(particle,stack);\r
525           mech = particle->GetUniqueID();\r
526           isPrim = stack->IsPhysicalPrimary(label);\r
527           isFromStrangess  = IsFromStrangeness(label,stack);\r
528           isFromConversion = IsFromConversion(label,stack);\r
529           isFromMaterial   = IsFromMaterial(label,stack);\r
530         }\r
531 \r
532         //\r
533         // TPC track\r
534         //\r
535         Int_t labelTPC = TMath::Abs(track->GetTPCLabel()); \r
536         particleTPC = stack->Particle(labelTPC);\r
537         if(particleTPC && particleTPC->GetPDG() && particleTPC->GetPDG()->Charge()!=0.)\r
538         {\r
539           particleMotherTPC = GetMother(particleTPC,stack);\r
540           mechTPC = particleTPC->GetUniqueID();\r
541           isPrimTPC = stack->IsPhysicalPrimary(labelTPC);\r
542           isFromStrangessTPC  = IsFromStrangeness(labelTPC,stack);\r
543           isFromConversionTPC = IsFromConversion(labelTPC,stack);\r
544           isFromMaterialTPC   = IsFromMaterial(labelTPC,stack);\r
545         }\r
546 \r
547         //\r
548         // store first track reference\r
549         // for TPC track\r
550         //\r
551         TParticle *part=0;\r
552         TClonesArray *trefs=0;\r
553         Int_t status = mcEvent->GetParticleAndTR(track->GetTPCLabel(), part, trefs);\r
554 \r
555         if(status>0 && part && trefs && part->GetPDG() && part->GetPDG()->Charge()!=0.) \r
556         {\r
557           Int_t nTrackRef = trefs->GetEntries();\r
558           //printf("nTrackRef %d \n",nTrackRef);\r
559 \r
560           Int_t countITS = 0;\r
561           for (Int_t iref = 0; iref < nTrackRef; iref++) \r
562           {\r
563             AliTrackReference *ref = (AliTrackReference *)trefs->At(iref);\r
564             //printf("ref %p \n",ref);\r
565             //if(ref) printf("ref->DetectorId() %d \n",ref->DetectorId());\r
566             //printf("AliTrackReference::kTPC  %d \n",AliTrackReference::kTPC);\r
567 \r
568 \r
569              // Ref. in the middle ITS \r
570             if(ref && ref->DetectorId()==AliTrackReference::kITS)\r
571             {\r
572               if(!refITS && countITS==2) {\r
573                 refITS = ref;\r
574                 //printf("refITS %p \n",refITS);\r
575               }\r
576               countITS++;\r
577             }\r
578 \r
579             // TPC\r
580             if(ref && ref->DetectorId()==AliTrackReference::kTPC)\r
581             {\r
582               if(!refTPCIn) {\r
583                 refTPCIn = ref;\r
584                 //printf("refTPCIn %p \n",refTPCIn);\r
585                 //break;\r
586               }\r
587             }\r
588           }\r
589 \r
590           // transform inner params to TrackRef\r
591           // reference frame\r
592           if(refTPCIn && trackInnerC3) \r
593           {\r
594             Double_t kRefPhi = TMath::ATan2(refTPCIn->Y(),refTPCIn->X());\r
595             isOK = kFALSE;\r
596             isOK = trackInnerC3->Rotate(kRefPhi);\r
597             isOK = AliTracker::PropagateTrackToBxByBz(trackInnerC3,refTPCIn->R(),track->GetMass(),kStep,kFALSE);\r
598 \r
599             if(!isOK){\r
600               if(trackInnerC3) delete trackInnerC3;\r
601               if(refTPCIn) delete refTPCIn;\r
602             }\r
603           }\r
604         }\r
605 \r
606         //\r
607         // ITS track\r
608         //\r
609         Int_t labelITS = TMath::Abs(track->GetITSLabel()); \r
610         particleITS = stack->Particle(labelITS);\r
611         if(particleITS && particleITS->GetPDG() && particleITS->GetPDG()->Charge()!=0.)\r
612         {\r
613           particleMotherITS = GetMother(particleITS,stack);\r
614           mechITS = particleITS->GetUniqueID();\r
615           isPrimITS = stack->IsPhysicalPrimary(labelITS);\r
616           isFromStrangessITS  = IsFromStrangeness(labelITS,stack);\r
617           isFromConversionITS = IsFromConversion(labelITS,stack);\r
618           isFromMaterialITS   = IsFromMaterial(labelITS,stack);\r
619         }\r
620       }\r
621 \r
622       //\r
623       Double_t vert[3] = {0}; \r
624       vert[0] = vtxESD->GetXv();\r
625       vert[1] = vtxESD->GetYv();\r
626       vert[2] = vtxESD->GetZv();\r
627       Int_t mult = vtxESD->GetNContributors();\r
628       Double_t bz = esdEvent->GetMagneticField();\r
629       Double_t runNumber = esdEvent->GetRunNumber();\r
630       Double_t evtTimeStamp = esdEvent->GetTimeStamp();\r
631       Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
632 \r
633 \r
634       //\r
635       if(!fTreeSRedirector) return;\r
636       (*fTreeSRedirector)<<"dNdPtTree"<<\r
637         "fileName.="<<&fileName<<\r
638         "runNumber="<<runNumber<<\r
639         "evtTimeStamp="<<evtTimeStamp<<\r
640         "evtNumberInFile="<<evtNumberInFile<<\r
641         "Bz="<<bz<<\r
642         "vertX="<<vert[0]<<\r
643         "vertY="<<vert[1]<<\r
644         "vertZ="<<vert[2]<<\r
645         "mult="<<mult<<\r
646         "esdTrack.="<<track<<\r
647         "extTPCInnerC.="<<tpcInnerC<<\r
648         "extInnerParamC.="<<trackInnerC<<\r
649         "extInnerParam.="<<trackInnerC2<<\r
650         "extOuterITS.="<<outerITSc<<\r
651         "extInnerParamRef.="<<trackInnerC3<<\r
652         "refTPCIn.="<<refTPCIn<<\r
653         "refITS.="<<refITS<<\r
654         "chi2TPCInnerC="<<chi2(0,0)<<\r
655         "chi2InnerC="<<chi2trackC(0,0)<<\r
656         "chi2OuterITS="<<chi2OuterITS(0,0)<<\r
657         "centralityF="<<centralityF<<\r
658         "particle.="<<particle<<\r
659         "particleMother.="<<particleMother<<\r
660         "mech="<<mech<<\r
661         "isPrim="<<isPrim<<\r
662         "isFromStrangess="<<isFromStrangess<<\r
663         "isFromConversion="<<isFromConversion<<\r
664         "isFromMaterial="<<isFromMaterial<<\r
665         "particleTPC.="<<particleTPC<<\r
666         "particleMotherTPC.="<<particleMotherTPC<<\r
667         "mechTPC="<<mechTPC<<\r
668         "isPrimTPC="<<isPrimTPC<<\r
669         "isFromStrangessTPC="<<isFromStrangessTPC<<\r
670         "isFromConversionTPC="<<isFromConversionTPC<<\r
671         "isFromMaterialTPC="<<isFromMaterialTPC<<\r
672         "particleITS.="<<particleITS<<\r
673         "particleMotherITS.="<<particleMotherITS<<\r
674         "mechITS="<<mechITS<<\r
675         "isPrimITS="<<isPrimITS<<\r
676         "isFromStrangessITS="<<isFromStrangessITS<<\r
677         "isFromConversionITS="<<isFromConversionITS<<\r
678         "isFromMaterialITS="<<isFromMaterialITS<<\r
679         "\n";\r
680 \r
681         if(tpcInnerC) delete tpcInnerC;\r
682         if(trackInnerC) delete trackInnerC;\r
683         if(trackInnerC2) delete trackInnerC2;\r
684         if(outerITSc) delete outerITSc;\r
685 \r
686         if(trackInnerC3) delete trackInnerC3;\r
687     }\r
688   }\r
689 \r
690   PostData(0, fOutputSummary);\r
691   //PostData(1, fOutput);\r
692 }\r
693 \r
694 //_____________________________________________________________________________\r
695 void AlidNdPtTrackDumpTask::ProcessV0(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)\r
696 {\r
697   //\r
698   // Process real and/or simulated events\r
699   //\r
700   if(!esdEvent) {\r
701     AliDebug(AliLog::kError, "esdEvent not available");\r
702     return;\r
703   }\r
704 \r
705   // get selection cuts\r
706   AlidNdPtEventCuts *evtCuts = GetEventCuts(); \r
707   AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
708   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
709 \r
710   if(!evtCuts || !accCuts  || !esdTrackCuts) {\r
711     AliDebug(AliLog::kError, "cuts not available");\r
712     return;\r
713   }\r
714 \r
715 \r
716 \r
717 \r
718   // trigger selection\r
719   Bool_t isEventTriggered = kTRUE;\r
720   AliPhysicsSelection *physicsSelection = NULL;\r
721   AliTriggerAnalysis* triggerAnalysis = NULL;\r
722 \r
723   // \r
724   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
725   if (!inputHandler)\r
726   {\r
727     Printf("ERROR: Could not receive input handler");\r
728     return;\r
729   }\r
730    \r
731   // get file name\r
732   TTree *chain = (TChain*)GetInputData(0);\r
733   if(!chain) { \r
734     Printf("ERROR: Could not receive input chain");\r
735     return;\r
736   }\r
737   TObjString fileName(chain->GetCurrentFile()->GetName());\r
738 \r
739   // trigger\r
740   if(evtCuts->IsTriggerRequired())  \r
741   {\r
742     // always MB\r
743     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
744 \r
745     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
746     if(!physicsSelection) return;\r
747     //SetPhysicsTriggerSelection(physicsSelection);\r
748 \r
749     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
750       // set trigger (V0AND)\r
751       triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
752       if(!triggerAnalysis) return;\r
753       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
754     }\r
755   }\r
756 \r
757   // centrality determination\r
758   Float_t centralityF = -1;\r
759   AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
760   centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
761 \r
762 \r
763   // get reconstructed vertex  \r
764   //const AliESDVertex* vtxESD = 0; \r
765   const AliESDVertex* vtxESD = 0; \r
766   if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {\r
767         vtxESD = esdEvent->GetPrimaryVertexTPC();\r
768   }\r
769   else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {\r
770      vtxESD = esdEvent->GetPrimaryVertexTracks();\r
771   }\r
772   else {\r
773         return;\r
774   }\r
775 \r
776   if(!vtxESD) return;\r
777 \r
778   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
779   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
780   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
781 \r
782   // check event cuts\r
783   if(isEventOK && isEventTriggered) {\r
784   //\r
785   // Dump the pt downscaled V0 into the tree\r
786   // \r
787   //\r
788   Int_t ntracks = esdEvent->GetNumberOfTracks();\r
789   Int_t nV0s = esdEvent->GetNumberOfV0s();\r
790   Int_t run = esdEvent->GetRunNumber();\r
791   Int_t time= esdEvent->GetTimeStamp();\r
792   Int_t evNr=esdEvent->GetEventNumberInFile();\r
793   \r
794   for (Int_t iv0=0; iv0<nV0s; iv0++){\r
795     AliESDv0 * v0 = esdEvent->GetV0(iv0);\r
796     if (!v0) continue;\r
797     AliESDtrack * track0 = esdEvent->GetTrack(v0->GetIndex(0));\r
798     AliESDtrack * track1 = esdEvent->GetTrack(v0->GetIndex(1));\r
799     if (!track0) continue;\r
800     if (!track1) continue;\r
801     if (track0->GetSign()<0) {\r
802       track1 = esdEvent->GetTrack(v0->GetIndex(0));\r
803       track0 = esdEvent->GetTrack(v0->GetIndex(1));\r
804     }\r
805     //\r
806     Bool_t isDownscaled = IsV0Downscaled(v0);\r
807     if (isDownscaled) continue;\r
808     AliKFParticle kfparticle; //\r
809     Int_t type=GetKFParticle(v0,esdEvent,kfparticle);\r
810     if (type==0) continue;   \r
811 \r
812     if(!fTreeSRedirector) return;\r
813     (*fTreeSRedirector)<<"V0s"<<\r
814       "isDownscaled="<<isDownscaled<<\r
815       "run="<<run<<\r
816       "fname="<<&fileName<<\r
817       "time="<<time<<\r
818       "evNr="<<evNr<<\r
819       "type="<<type<<\r
820       "ntracks="<<ntracks<<\r
821       "v0.="<<v0<<\r
822       "kf.="<<&kfparticle<<\r
823       "track0.="<<track0<<\r
824       "track1.="<<track1<<\r
825       "centralityF="<<centralityF<<\r
826       "\n";\r
827   }\r
828   }\r
829   PostData(0, fOutputSummary);\r
830 }\r
831 \r
832 //_____________________________________________________________________________\r
833 Int_t   AlidNdPtTrackDumpTask::GetKFParticle(AliESDv0 *const v0, AliESDEvent * const event, AliKFParticle & kfparticle)\r
834 {\r
835   //\r
836   // Create KF particle in case the V0 fullfill selection criteria\r
837   //\r
838   // Selection criteria\r
839   //  0. algorithm cut\r
840   //  1. track cut\r
841   //  3. chi2 cut\r
842   //  4. rough mass cut\r
843   //  5. Normalized pointing angle cut\r
844   //\r
845   const Double_t cutMass=0.2;\r
846   const Double_t kSigmaDCACut=3;\r
847   //\r
848   // 0.) algo cut - accept only on the fly\r
849   //\r
850   if (v0->GetOnFlyStatus() ==kFALSE) return 0;     \r
851   //\r
852   // 1.) track cut\r
853   // \r
854   AliESDtrack * track0 = event->GetTrack(v0->GetIndex(0));\r
855   AliESDtrack * track1 = event->GetTrack(v0->GetIndex(1));\r
856   /*\r
857     TCut cutD="abs(track0.fD/sqrt(track0.fCdd))>2&&abs(track1.fD/sqrt(track1.fCdd))>2";\r
858     TCut cutTheta="abs(track0.fP[3])<1&&abs(track1.fP[3])<1";\r
859     TCut cutNcl="track0.GetTPCClusterInfo(2,1)>100&&track1.GetTPCClusterInfo(2,1)>100";\r
860   */  \r
861   if (TMath::Abs(track0->GetTgl())>1) return 0;\r
862   if (TMath::Abs(track1->GetTgl())>1) return 0;\r
863   if ((track0->GetTPCClusterInfo(2,1))<100) return 0;\r
864   if ((track1->GetTPCClusterInfo(2,1))<100) return 0;\r
865   //if ((track0->GetITSclusters(0))<2) return 0;\r
866   //if ((track1->GetITSclusters(0))<2) return 0; \r
867   Float_t pos0[2]={0}, cov0[3]={0};\r
868   Float_t pos1[2]={0}, cov1[3]={0};\r
869   track0->GetImpactParameters(pos0,cov0);\r
870   track0->GetImpactParameters(pos1,cov1);\r
871   //\r
872   if (TMath::Abs(pos0[0])<kSigmaDCACut*TMath::Sqrt(cov0[0])) return 0;\r
873   if (TMath::Abs(pos1[0])<kSigmaDCACut*TMath::Sqrt(cov1[0])) return 0;\r
874   // \r
875   //\r
876   // 3.) Chi2 cut\r
877   //\r
878   Double_t chi2KF = v0->GetKFInfo(2,2,2);\r
879   if (chi2KF>25) return 0;\r
880   //\r
881   // 4.) Rough mass cut - 0.200 GeV\r
882   //\r
883   static Double_t masses[2]={-1};\r
884   if (masses[0]<0){\r
885     masses[0] = TDatabasePDG::Instance()->GetParticle("K_S0")->Mass();\r
886     masses[1] = TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass();\r
887   }\r
888   Double_t mass00=  v0->GetEffMass(0,0);\r
889   Double_t mass22=  v0->GetEffMass(2,2);\r
890   Double_t mass42=  v0->GetEffMass(4,2);\r
891   Double_t mass24=  v0->GetEffMass(2,4);\r
892   Bool_t massOK=kFALSE;\r
893   Int_t type=0;\r
894   Int_t ptype=0;\r
895   Double_t dmass=1;\r
896   Int_t p1=0, p2=0;\r
897   if (TMath::Abs(mass00-0)<cutMass) {\r
898     massOK=kTRUE; type+=1; \r
899     if (TMath::Abs(mass00-0)<dmass) {\r
900       ptype=1;\r
901       dmass=TMath::Abs(mass00-0);      \r
902       p1=0; p2=0;\r
903     } \r
904   }\r
905   if (TMath::Abs(mass24-masses[1])<cutMass) {\r
906     massOK=kTRUE; type+=2; \r
907     if (TMath::Abs(mass24-masses[1])<dmass){\r
908       dmass = TMath::Abs(mass24-masses[1]);\r
909       ptype=2;\r
910       p1=2; p2=4;\r
911     }\r
912   }\r
913   if (TMath::Abs(mass42-masses[1])<cutMass) {\r
914     massOK=kTRUE; type+=4;\r
915     if (TMath::Abs(mass42-masses[1])<dmass){\r
916       dmass = TMath::Abs(mass42-masses[1]);\r
917       ptype=4;\r
918       p1=4; p2=2;\r
919     }\r
920   }\r
921   if (TMath::Abs(mass22-masses[0])<cutMass) {\r
922     massOK=kTRUE; type+=8;\r
923     if (TMath::Abs(mass22-masses[0])<dmass){\r
924       dmass = TMath::Abs(mass22-masses[0]);\r
925       ptype=8;\r
926       p1=2; p2=2;\r
927     }\r
928   }\r
929   if (type==0) return 0;\r
930   //\r
931   const Int_t spdg[5]={kPositron,kMuonPlus,kPiPlus, kKPlus, kProton};\r
932   const AliExternalTrackParam *paramP = v0->GetParamP();\r
933   const AliExternalTrackParam *paramN = v0->GetParamN();\r
934   if (paramP->GetSign()<0){\r
935     paramP=v0->GetParamP();\r
936     paramN=v0->GetParamN();\r
937   }\r
938   //Double_t *pparam1 = (Double_t*)paramP->GetParameter();\r
939   //Double_t *pparam2 = (Double_t*)paramN->GetParameter();\r
940   //\r
941   AliKFParticle kfp1( *paramP, spdg[p1]  );\r
942   AliKFParticle kfp2( *paramN, -1 *spdg[p2]  );\r
943   AliKFParticle V0KF;\r
944   (V0KF)+=kfp1;\r
945   (V0KF)+=kfp2;\r
946   kfparticle=V0KF;\r
947   //\r
948   // Pointing angle\r
949   //\r
950   Double_t  errPhi    = V0KF.GetErrPhi();\r
951   Double_t  pointAngle= TMath::ACos(v0->GetV0CosineOfPointingAngle());\r
952   if (pointAngle/errPhi>10) return 0;  \r
953   //\r
954   return ptype;  \r
955 }\r
956 \r
957 //_____________________________________________________________________________\r
958 Bool_t AlidNdPtTrackDumpTask::IsV0Downscaled(AliESDv0 *const v0)\r
959 {\r
960   //\r
961   // Downscale randomly low pt V0\r
962   //\r
963   //return kFALSE;\r
964   Double_t maxPt= TMath::Max(v0->GetParamP()->Pt(), v0->GetParamN()->Pt());\r
965   Double_t scalempt= TMath::Min(maxPt,10.);\r
966   if (TMath::Exp(2*scalempt)<fLowPtV0DownscaligF*gRandom->Rndm()) return kTRUE;\r
967   return kFALSE;\r
968   /*\r
969   \r
970     TH1F his1("his1","his1",100,0,10);\r
971     TH1F his2("his2","his2",100,0,10);\r
972     {for (Int_t i=0; i<10000; i++){\r
973        Double_t rnd=gRandom->Exp(1);\r
974        Bool_t isDownscaled =TMath::Exp(rnd)<100*gRandom->Rndm();\r
975        his1->Fill(rnd); \r
976        if (!isDownscaled) his2->Fill(rnd); \r
977     }}\r
978 \r
979    */\r
980 \r
981 }\r
982 \r
983 \r
984 \r
985 \r
986 \r
987 \r
988 \r
989 \r
990 \r
991 //_____________________________________________________________________________\r
992 Bool_t AlidNdPtTrackDumpTask::ConstrainTPCInner(AliExternalTrackParam *const tpcInnerC, const AliESDVertex* vtx, Double_t b[3])\r
993 {\r
994  // Constrain TPC inner params constrained\r
995  //\r
996       if(!tpcInnerC) return kFALSE; \r
997       if(!vtx) return kFALSE;\r
998 \r
999       Double_t dz[2],cov[3];\r
1000       //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();\r
1001       //if(!tpcInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; \r
1002       //if(!tpcInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; \r
1003       if(!tpcInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; \r
1004 \r
1005 \r
1006       Double_t covar[6]; vtx->GetCovMatrix(covar);\r
1007       Double_t p[2]={tpcInnerC->GetParameter()[0]-dz[0],tpcInnerC->GetParameter()[1]-dz[1]};\r
1008       Double_t c[3]={covar[2],0.,covar[5]};\r
1009       Double_t chi2C=tpcInnerC->GetPredictedChi2(p,c);\r
1010       if (chi2C>kVeryBig) return kFALSE; \r
1011 \r
1012       if(!tpcInnerC->Update(p,c)) return kFALSE;\r
1013 \r
1014   return kTRUE;\r
1015 }\r
1016 \r
1017 //_____________________________________________________________________________\r
1018 Bool_t AlidNdPtTrackDumpTask::ConstrainTrackInner(AliExternalTrackParam *const trackInnerC, const AliESDVertex* vtx, Double_t mass, Double_t b[3])\r
1019 {\r
1020  // Constrain TPC inner params constrained\r
1021  //\r
1022       if(!trackInnerC) return kFALSE; \r
1023       if(!vtx) return kFALSE;\r
1024 \r
1025       const Double_t kRadius  = 2.8; \r
1026       const Double_t kMaxStep = 1.0; \r
1027 \r
1028       Double_t dz[2],cov[3];\r
1029 \r
1030       //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();\r
1031       //if(!trackInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; \r
1032       //if(!trackInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; \r
1033 \r
1034       if(!AliTracker::PropagateTrackToBxByBz(trackInnerC,kRadius,mass,kMaxStep,kFALSE)) return kFALSE;\r
1035       if(!trackInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; \r
1036 \r
1037       //\r
1038       Double_t covar[6]; vtx->GetCovMatrix(covar);\r
1039       Double_t p[2]={trackInnerC->GetParameter()[0]-dz[0],trackInnerC->GetParameter()[1]-dz[1]};\r
1040       Double_t c[3]={covar[2],0.,covar[5]};\r
1041       Double_t chi2C=trackInnerC->GetPredictedChi2(p,c);\r
1042       if (chi2C>kVeryBig) return kFALSE; \r
1043 \r
1044       if(!trackInnerC->Update(p,c)) return kFALSE;\r
1045 \r
1046   return kTRUE;\r
1047 }\r
1048 \r
1049 \r
1050 //_____________________________________________________________________________\r
1051 TParticle *AlidNdPtTrackDumpTask::GetMother(TParticle *const particle, AliStack *const stack) \r
1052 {\r
1053   if(!particle) return NULL;\r
1054   if(!stack) return NULL;\r
1055 \r
1056   Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
1057   TParticle* mother = NULL; \r
1058   mother = stack->Particle(motherLabel); \r
1059 \r
1060 return mother;\r
1061 }\r
1062 \r
1063 //_____________________________________________________________________________\r
1064 Bool_t AlidNdPtTrackDumpTask::IsFromConversion(const Int_t label, AliStack *const stack) \r
1065 {\r
1066   Bool_t isFromConversion = kFALSE;\r
1067 \r
1068   if(stack) {\r
1069     TParticle* particle = stack->Particle(label);\r
1070 \r
1071     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
1072     {\r
1073        Int_t mech = particle->GetUniqueID(); // production mechanism \r
1074        Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
1075 \r
1076        Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
1077        TParticle* mother = stack->Particle(motherLabel); \r
1078        if(mother) {\r
1079           Int_t motherPdg = mother->GetPdgCode();\r
1080 \r
1081           if(!isPrim && mech==5 && motherPdg==kGamma) { \r
1082             isFromConversion=kTRUE; \r
1083           }\r
1084        }\r
1085     } \r
1086   } \r
1087 \r
1088   return isFromConversion;\r
1089 }\r
1090 \r
1091 //_____________________________________________________________________________\r
1092 Bool_t AlidNdPtTrackDumpTask::IsFromMaterial(const Int_t label, AliStack *const stack) \r
1093 {\r
1094   Bool_t isFromMaterial = kFALSE;\r
1095 \r
1096   if(stack) {\r
1097     TParticle* particle = stack->Particle(label);\r
1098 \r
1099     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
1100     {\r
1101        Int_t mech = particle->GetUniqueID(); // production mechanism \r
1102        Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
1103 \r
1104        Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
1105        TParticle* mother = stack->Particle(motherLabel); \r
1106        if(mother) {\r
1107           if(!isPrim && mech==13) { \r
1108             isFromMaterial=kTRUE; \r
1109           }\r
1110        }\r
1111      } \r
1112   } \r
1113 \r
1114   return isFromMaterial;\r
1115 }\r
1116 \r
1117 //_____________________________________________________________________________\r
1118 Bool_t AlidNdPtTrackDumpTask::IsFromStrangeness(const Int_t label, AliStack *const stack) \r
1119 {\r
1120   Bool_t isFromStrangeness = kFALSE;\r
1121 \r
1122   if(stack) {\r
1123     TParticle* particle = stack->Particle(label);\r
1124 \r
1125     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
1126     {\r
1127        Int_t mech = particle->GetUniqueID(); // production mechanism \r
1128        Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
1129 \r
1130        Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
1131        TParticle* mother = stack->Particle(motherLabel); \r
1132        if(mother) {\r
1133           Int_t motherPdg = mother->GetPdgCode();\r
1134 \r
1135           // K+-, lambda, antilambda, K0s decays\r
1136           if(!isPrim && mech==4 && \r
1137               (TMath::Abs(motherPdg)==kKPlus || TMath::Abs(motherPdg)==kLambda0 || motherPdg==kK0Short))\r
1138           {\r
1139             isFromStrangeness = kTRUE;\r
1140           } \r
1141        }\r
1142     } \r
1143   } \r
1144 \r
1145   return isFromStrangeness;\r
1146 }\r
1147 \r
1148 \r
1149 //_____________________________________________________________________________\r
1150 void AlidNdPtTrackDumpTask::FinishTaskOutput() \r
1151 {\r
1152   //\r
1153   // Called one at the end \r
1154   // locally on working node\r
1155   //\r
1156 \r
1157   // Post output data.\r
1158   PostData(1, fOutput);\r
1159   //PostData(0, fOutputSummary);\r
1160 }\r
1161 \r
1162 //_____________________________________________________________________________\r
1163 void AlidNdPtTrackDumpTask::Terminate(Option_t *) \r
1164 {\r
1165   // Called one at the end \r
1166   if(fTreeSRedirector)  delete fTreeSRedirector; fTreeSRedirector=0;\r
1167   fOutputSummary = dynamic_cast<TTree*> (GetOutputData(0));\r
1168   if(fOutputSummary) delete fOutputSummary; fOutputSummary=0;\r
1169 \r
1170   TChain* chain = new TChain("dNdPtTree");\r
1171   if(!chain) return;\r
1172   chain->Add("jotwinow_HighPt_TrackAndV0_Trees.root");\r
1173   TTree *tree = chain->CopyTree("1");\r
1174   if (chain) { delete chain; chain=0; }\r
1175   if(!tree) return;\r
1176   tree->Print();\r
1177   fOutputSummary = tree;\r
1178 \r
1179   if (!fOutputSummary) {\r
1180     Printf("ERROR: AlidNdPtTrackDumpTask::Terminate(): Output data not avaiable GetOutputData(0)==0x0 ..." );\r
1181     return;\r
1182   }\r
1183 \r
1184   PostData(0, fOutputSummary);\r
1185   //PostData(1, fOutput);\r
1186 }\r