Fix in AliTrackerBase::PropagateTo... routines: length integral was not correctly...
[u/mrichter/AliRoot.git] / PWGUD / 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 "dNdPt/AlidNdPt.h"\r
54 #include "dNdPt/AlidNdPtEventCuts.h"\r
55 #include "dNdPt/AlidNdPtAcceptanceCuts.h"\r
56 \r
57 #include "dNdPt/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   if(evtCuts->IsTriggerRequired())  \r
205   {\r
206     // always MB\r
207     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
208 \r
209     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
210     if(!physicsSelection) return;\r
211     //SetPhysicsTriggerSelection(physicsSelection);\r
212 \r
213     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
214       // set trigger (V0AND)\r
215       triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
216       if(!triggerAnalysis) return;\r
217       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
218     }\r
219   }\r
220 \r
221   // centrality determination\r
222   Float_t centralityF = -1;\r
223   AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
224   centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
225 \r
226   // use MC information\r
227   AliHeader* header = 0;\r
228   AliGenEventHeader* genHeader = 0;\r
229   AliStack* stack = 0;\r
230   TArrayF vtxMC(3);\r
231 \r
232   Int_t multMCTrueTracks = 0;\r
233   if(IsUseMCInfo())\r
234   {\r
235     //\r
236     if(!mcEvent) {\r
237       AliDebug(AliLog::kError, "mcEvent not available");\r
238       return;\r
239     }\r
240     // get MC event header\r
241     header = mcEvent->Header();\r
242     if (!header) {\r
243       AliDebug(AliLog::kError, "Header not available");\r
244       return;\r
245     }\r
246     // MC particle stack\r
247     stack = mcEvent->Stack();\r
248     if (!stack) {\r
249       AliDebug(AliLog::kError, "Stack not available");\r
250       return;\r
251     }\r
252 \r
253     // get MC vertex\r
254     genHeader = header->GenEventHeader();\r
255     if (!genHeader) {\r
256       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");\r
257       return;\r
258     }\r
259     genHeader->PrimaryVertex(vtxMC);\r
260 \r
261     // multipliticy of all MC primary tracks\r
262     // in Zv, pt and eta ranges)\r
263     multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);\r
264 \r
265   } // end bUseMC\r
266 \r
267   // get reconstructed vertex  \r
268   //const AliESDVertex* vtxESD = 0; \r
269   const AliESDVertex* vtxESD = 0; \r
270   if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {\r
271         vtxESD = esdEvent->GetPrimaryVertexTPC();\r
272   }\r
273   else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {\r
274      vtxESD = esdEvent->GetPrimaryVertexTracks();\r
275   }\r
276   else {\r
277         return;\r
278   }\r
279 \r
280   if(!vtxESD) return;\r
281 \r
282   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
283   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
284   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
285 \r
286   // check event cuts\r
287   if(isEventOK && isEventTriggered)\r
288   {\r
289     TRandom random;\r
290 \r
291     for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
292     {\r
293       AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
294       if(!track) continue;\r
295       if(track->Charge()==0) continue;\r
296       if(!esdTrackCuts->AcceptTrack(track)) continue;\r
297       if(!accCuts->AcceptTrack(track)) continue;\r
298 \r
299       // downscale low-pT tracks\r
300       if(TMath::Exp(2*track->Pt())<1000*random.Rndm()) continue;\r
301 \r
302       // Dump to the tree \r
303       // vertex\r
304       // TPC constrained tracks\r
305       // InnerParams constrained tracks\r
306       // TPC-ITS tracks\r
307       // ITSout-InnerParams tracks\r
308       // chi2 distance between TPC constrained and TPC-ITS tracks\r
309       // chi2 distance between TPC refitted constrained and TPC-ITS tracks\r
310       // chi2 distance between ITSout and InnerParams tracks\r
311       // MC information\r
312       \r
313       Double_t x[3]; track->GetXYZ(x);\r
314       Double_t b[3]; AliTracker::GetBxByBz(x,b);\r
315       Bool_t isOK = kFALSE;\r
316 \r
317       //\r
318       // Constrain TPC-only tracks (TPCinner) to vertex\r
319       //\r
320       AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());\r
321       if (!tpcInner) continue;\r
322       // transform to the track reference frame \r
323       isOK = tpcInner->Rotate(track->GetAlpha());\r
324       isOK = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
325       if(!isOK) continue;\r
326 \r
327       // clone TPCinner has to be deleted\r
328       AliExternalTrackParam * tpcInnerC = new AliExternalTrackParam(*(track->GetTPCInnerParam()));\r
329       if (!tpcInnerC) continue;\r
330  \r
331       // constrain TPCinner \r
332       //isOK = ConstrainTPCInner(tpcInnerC,vtxESD,esdEvent->GetMagneticField());\r
333       isOK = ConstrainTPCInner(tpcInnerC,vtxESD,b);\r
334 \r
335       // transform to the track reference frame \r
336       isOK = tpcInnerC->Rotate(track->GetAlpha());\r
337       isOK = tpcInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
338 \r
339       if(!isOK) {\r
340         if(tpcInnerC) delete tpcInnerC;\r
341         continue;\r
342       }\r
343 \r
344 \r
345       //\r
346       // Constrain TPC refitted tracks at inner TPC wall (InnerParams) \r
347       // to vertex\r
348       //\r
349       // clone track InnerParams has to be deleted\r
350       AliExternalTrackParam * trackInnerC =  new AliExternalTrackParam(*(track->GetInnerParam()));\r
351       if (!trackInnerC) continue;\r
352  \r
353       // constrain track InnerParams \r
354       isOK = ConstrainTrackInner(trackInnerC,vtxESD,track->GetMass(),b);\r
355 \r
356       // transform to the track reference frame \r
357       isOK = trackInnerC->Rotate(track->GetAlpha());\r
358       isOK = trackInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
359 \r
360       if(!isOK) {\r
361         if(trackInnerC) delete trackInnerC;\r
362         continue;\r
363       }\r
364       \r
365       //\r
366       // calculate chi2 between vi and vj vectors\r
367       // with covi and covj covariance matrices\r
368       // chi2ij = (vi-vj)^(T)*(covi+covj)^(-1)*(vi-vj)\r
369       //\r
370       TMatrixD deltaT(5,1), deltaTtrackC(5,1);\r
371       TMatrixD delta(1,5),  deltatrackC(1,5);\r
372       TMatrixD covarM(5,5), covarMtrackC(5,5);\r
373 \r
374       for (Int_t ipar=0; ipar<5; ipar++) {\r
375         deltaT(ipar,0)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
376         delta(0,ipar)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
377 \r
378         deltaTtrackC(ipar,0)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
379         deltatrackC(0,ipar)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
380 \r
381         for (Int_t jpar=0; jpar<5; jpar++) {\r
382           Int_t index=track->GetIndex(ipar,jpar);\r
383           covarM(ipar,jpar)=track->GetCovariance()[index]+tpcInnerC->GetCovariance()[index];\r
384           covarMtrackC(ipar,jpar)=track->GetCovariance()[index]+trackInnerC->GetCovariance()[index];\r
385         }\r
386       }\r
387       // chi2 distance TPC constrained and TPC+ITS\r
388       TMatrixD covarMInv = covarM.Invert();\r
389       TMatrixD mat2 = covarMInv*deltaT;\r
390       TMatrixD chi2 = delta*mat2; \r
391       //chi2.Print();\r
392 \r
393       // chi2 distance TPC refitted constrained and TPC+ITS\r
394       TMatrixD covarMInvtrackC = covarMtrackC.Invert();\r
395       TMatrixD mat2trackC = covarMInvtrackC*deltaTtrackC;\r
396       TMatrixD chi2trackC = deltatrackC*mat2trackC; \r
397       //chi2trackC.Print();\r
398 \r
399 \r
400       //\r
401       // Propagate ITSout to TPC inner wall \r
402       // and calculate chi2 distance to track (InnerParams)\r
403       //\r
404       const Double_t kTPCRadius=85; \r
405       const Double_t kStep=3; \r
406 \r
407       // clone track InnerParams has to be deleted\r
408       AliExternalTrackParam *trackInnerC2 = new AliExternalTrackParam(*(track->GetInnerParam()));\r
409       if (!trackInnerC2) continue;\r
410       if(!AliTracker::PropagateTrackToBxByBz(trackInnerC2,kTPCRadius,track->GetMass(),kStep,kFALSE))\r
411       {\r
412           if(trackInnerC2) delete trackInnerC2;\r
413           continue;\r
414       }\r
415 \r
416       AliExternalTrackParam *outerITSc = new AliExternalTrackParam();\r
417       if(!outerITSc) continue;\r
418 \r
419       TMatrixD chi2OuterITS(1,1);\r
420       chi2OuterITS(0,0) = 0;\r
421 \r
422       if(esdFriend && esdFriend->TestSkipBit()==kFALSE) \r
423       {\r
424         // propagate ITSout to TPC inner wall\r
425         AliESDfriendTrack *friendTrack = esdFriend->GetTrack(iTrack);\r
426 \r
427         if(friendTrack) \r
428         {\r
429           if( (outerITSc = new AliExternalTrackParam(*friendTrack->GetITSOut())) ) \r
430           {\r
431             if(AliTracker::PropagateTrackToBxByBz(outerITSc,kTPCRadius,track->GetMass(),kStep,kFALSE))\r
432             {\r
433               // transform ITSout to the track InnerParams reference frame \r
434               isOK = kFALSE;\r
435               isOK = outerITSc->Rotate(trackInnerC2->GetAlpha());\r
436               isOK = outerITSc->PropagateTo(trackInnerC2->GetX(),esdEvent->GetMagneticField());\r
437 \r
438               if(!isOK) {\r
439                 if(outerITSc) delete outerITSc;\r
440                 if(trackInnerC2) delete trackInnerC2;\r
441                 continue;\r
442               }\r
443               \r
444               //\r
445               // calculate chi2 between outerITS and innerParams\r
446               // cov without z-coordinate at the moment\r
447               //\r
448               TMatrixD deltaTouterITS(4,1);\r
449               TMatrixD deltaouterITS(1,4);\r
450               TMatrixD covarMouterITS(4,4);\r
451 \r
452               Int_t kipar = 0;\r
453               Int_t kjpar = 0;\r
454               for (Int_t ipar=0; ipar<5; ipar++) {\r
455                 if(ipar!=1) {\r
456                   deltaTouterITS(kipar,0)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];\r
457                   deltaouterITS(0,kipar)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];\r
458                 }\r
459 \r
460                 kjpar=0;\r
461                 for (Int_t jpar=0; jpar<5; jpar++) {\r
462                   Int_t index=outerITSc->GetIndex(ipar,jpar);\r
463                   if(ipar !=1 || jpar!=1) {\r
464                     covarMouterITS(kipar,kjpar)=outerITSc->GetCovariance()[index]+trackInnerC2->GetCovariance()[index];\r
465                   }\r
466                   if(jpar!=1)  kjpar++;\r
467                 }\r
468                 if(ipar!=1) kipar++;\r
469               }\r
470 \r
471               // chi2 distance ITSout and InnerParams\r
472               TMatrixD covarMInvouterITS = covarMouterITS.Invert();\r
473               TMatrixD mat2outerITS = covarMInvouterITS*deltaTouterITS;\r
474               chi2OuterITS = deltaouterITS*mat2outerITS; \r
475               //chi2OuterITS.Print();\r
476             } \r
477           }\r
478         }\r
479       }\r
480 \r
481       //\r
482       // MC info\r
483       //\r
484       TParticle *particle=NULL, *particleTPC=NULL, *particleITS=NULL;\r
485       TParticle *particleMother=NULL, *particleMotherTPC=NULL, *particleMotherITS=NULL;\r
486       Int_t mech=-1, mechTPC=-1, mechITS=-1;\r
487       Bool_t isPrim=kFALSE, isPrimTPC=kFALSE, isPrimITS=kFALSE;\r
488       Bool_t isFromStrangess=kFALSE, isFromStrangessTPC=kFALSE, isFromStrangessITS=kFALSE;\r
489       Bool_t isFromConversion=kFALSE, isFromConversionTPC=kFALSE, isFromConversionITS=kFALSE;\r
490       Bool_t isFromMaterial=kFALSE, isFromMaterialTPC=kFALSE, isFromMaterialITS=kFALSE;\r
491 \r
492       AliTrackReference *refTPCIn = NULL;\r
493       AliTrackReference *refITS = NULL;\r
494       AliExternalTrackParam *trackInnerC3 = new AliExternalTrackParam(*(track->GetInnerParam()));\r
495       if(!trackInnerC3) continue;\r
496 \r
497       if(IsUseMCInfo()) \r
498       {\r
499         if(!stack) return;\r
500 \r
501         //\r
502         // global track\r
503         //\r
504         Int_t label = TMath::Abs(track->GetLabel()); \r
505         particle = stack->Particle(label);\r
506         if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0.)\r
507         {\r
508           particleMother = GetMother(particle,stack);\r
509           mech = particle->GetUniqueID();\r
510           isPrim = stack->IsPhysicalPrimary(label);\r
511           isFromStrangess  = IsFromStrangeness(label,stack);\r
512           isFromConversion = IsFromConversion(label,stack);\r
513           isFromMaterial   = IsFromMaterial(label,stack);\r
514         }\r
515 \r
516         //\r
517         // TPC track\r
518         //\r
519         Int_t labelTPC = TMath::Abs(track->GetTPCLabel()); \r
520         particleTPC = stack->Particle(labelTPC);\r
521         if(particleTPC && particleTPC->GetPDG() && particleTPC->GetPDG()->Charge()!=0.)\r
522         {\r
523           particleMotherTPC = GetMother(particleTPC,stack);\r
524           mechTPC = particleTPC->GetUniqueID();\r
525           isPrimTPC = stack->IsPhysicalPrimary(labelTPC);\r
526           isFromStrangessTPC  = IsFromStrangeness(labelTPC,stack);\r
527           isFromConversionTPC = IsFromConversion(labelTPC,stack);\r
528           isFromMaterialTPC   = IsFromMaterial(labelTPC,stack);\r
529         }\r
530 \r
531         //\r
532         // store first track reference\r
533         // for TPC track\r
534         //\r
535         TParticle *part=0;\r
536         TClonesArray *trefs=0;\r
537         Int_t status = mcEvent->GetParticleAndTR(track->GetTPCLabel(), part, trefs);\r
538 \r
539         if(status>0 && part && trefs && part->GetPDG() && part->GetPDG()->Charge()!=0.) \r
540         {\r
541           Int_t nTrackRef = trefs->GetEntries();\r
542           //printf("nTrackRef %d \n",nTrackRef);\r
543 \r
544           Int_t countITS = 0;\r
545           for (Int_t iref = 0; iref < nTrackRef; iref++) \r
546           {\r
547             AliTrackReference *ref = (AliTrackReference *)trefs->At(iref);\r
548             //printf("ref %p \n",ref);\r
549             //if(ref) printf("ref->DetectorId() %d \n",ref->DetectorId());\r
550             //printf("AliTrackReference::kTPC  %d \n",AliTrackReference::kTPC);\r
551 \r
552 \r
553              // Ref. in the middle ITS \r
554             if(ref && ref->DetectorId()==AliTrackReference::kITS)\r
555             {\r
556               if(!refITS && countITS==2) {\r
557                 refITS = ref;\r
558                 //printf("refITS %p \n",refITS);\r
559               }\r
560               countITS++;\r
561             }\r
562 \r
563             // TPC\r
564             if(ref && ref->DetectorId()==AliTrackReference::kTPC)\r
565             {\r
566               if(!refTPCIn) {\r
567                 refTPCIn = ref;\r
568                 //printf("refTPCIn %p \n",refTPCIn);\r
569                 //break;\r
570               }\r
571             }\r
572           }\r
573 \r
574           // transform inner params to TrackRef\r
575           // reference frame\r
576           if(refTPCIn && trackInnerC3) \r
577           {\r
578             Double_t kRefPhi = TMath::ATan2(refTPCIn->Y(),refTPCIn->X());\r
579             isOK = kFALSE;\r
580             isOK = trackInnerC3->Rotate(kRefPhi);\r
581             isOK = AliTracker::PropagateTrackToBxByBz(trackInnerC3,refTPCIn->R(),track->GetMass(),kStep,kFALSE);\r
582 \r
583             if(!isOK){\r
584               if(trackInnerC3) delete trackInnerC3;\r
585               if(refTPCIn) delete refTPCIn;\r
586             }\r
587           }\r
588         }\r
589 \r
590         //\r
591         // ITS track\r
592         //\r
593         Int_t labelITS = TMath::Abs(track->GetITSLabel()); \r
594         particleITS = stack->Particle(labelITS);\r
595         if(particleITS && particleITS->GetPDG() && particleITS->GetPDG()->Charge()!=0.)\r
596         {\r
597           particleMotherITS = GetMother(particleITS,stack);\r
598           mechITS = particleITS->GetUniqueID();\r
599           isPrimITS = stack->IsPhysicalPrimary(labelITS);\r
600           isFromStrangessITS  = IsFromStrangeness(labelITS,stack);\r
601           isFromConversionITS = IsFromConversion(labelITS,stack);\r
602           isFromMaterialITS   = IsFromMaterial(labelITS,stack);\r
603         }\r
604       }\r
605 \r
606       //\r
607       Double_t vert[3] = {0}; \r
608       vert[0] = vtxESD->GetXv();\r
609       vert[1] = vtxESD->GetYv();\r
610       vert[2] = vtxESD->GetZv();\r
611       Int_t mult = vtxESD->GetNContributors();\r
612       Double_t bz = esdEvent->GetMagneticField();\r
613       Double_t runNumber = esdEvent->GetRunNumber();\r
614       Double_t evtTimeStamp = esdEvent->GetTimeStamp();\r
615 \r
616       //\r
617       if(!fTreeSRedirector) return;\r
618       (*fTreeSRedirector)<<"dNdPtTree"<<\r
619         "runNumber="<<runNumber<<\r
620         "evtTimeStamp="<<evtTimeStamp<<\r
621         "Bz="<<bz<<\r
622         "vertX="<<vert[0]<<\r
623         "vertY="<<vert[1]<<\r
624         "vertZ="<<vert[2]<<\r
625         "mult="<<mult<<\r
626         "esdTrack.="<<track<<\r
627         "extTPCInnerC.="<<tpcInnerC<<\r
628         "extInnerParamC.="<<trackInnerC<<\r
629         "extInnerParam.="<<trackInnerC2<<\r
630         "extOuterITS.="<<outerITSc<<\r
631         "extInnerParamRef.="<<trackInnerC3<<\r
632         "refTPCIn.="<<refTPCIn<<\r
633         "refITS.="<<refITS<<\r
634         "chi2TPCInnerC="<<chi2(0,0)<<\r
635         "chi2InnerC="<<chi2trackC(0,0)<<\r
636         "chi2OuterITS="<<chi2OuterITS(0,0)<<\r
637         "centralityF="<<centralityF<<\r
638         "particle.="<<particle<<\r
639         "particleMother.="<<particleMother<<\r
640         "mech="<<mech<<\r
641         "isPrim="<<isPrim<<\r
642         "isFromStrangess="<<isFromStrangess<<\r
643         "isFromConversion="<<isFromConversion<<\r
644         "isFromMaterial="<<isFromMaterial<<\r
645         "particleTPC.="<<particleTPC<<\r
646         "particleMotherTPC.="<<particleMotherTPC<<\r
647         "mechTPC="<<mechTPC<<\r
648         "isPrimTPC="<<isPrimTPC<<\r
649         "isFromStrangessTPC="<<isFromStrangessTPC<<\r
650         "isFromConversionTPC="<<isFromConversionTPC<<\r
651         "isFromMaterialTPC="<<isFromMaterialTPC<<\r
652         "particleITS.="<<particleITS<<\r
653         "particleMotherITS.="<<particleMotherITS<<\r
654         "mechITS="<<mechITS<<\r
655         "isPrimITS="<<isPrimITS<<\r
656         "isFromStrangessITS="<<isFromStrangessITS<<\r
657         "isFromConversionITS="<<isFromConversionITS<<\r
658         "isFromMaterialITS="<<isFromMaterialITS<<\r
659         "\n";\r
660 \r
661         if(tpcInnerC) delete tpcInnerC;\r
662         if(trackInnerC) delete trackInnerC;\r
663         if(trackInnerC2) delete trackInnerC2;\r
664         if(outerITSc) delete outerITSc;\r
665 \r
666         if(trackInnerC3) delete trackInnerC3;\r
667     }\r
668   }\r
669 \r
670   PostData(0, fOutputSummary);\r
671   //PostData(1, fOutput);\r
672 }\r
673 \r
674 \r
675 //_____________________________________________________________________________\r
676 Bool_t AlidNdPtTrackDumpTask::ConstrainTPCInner(AliExternalTrackParam *const tpcInnerC, const AliESDVertex* vtx, Double_t b[3])\r
677 {\r
678  // Constrain TPC inner params constrained\r
679  //\r
680       if(!tpcInnerC) return kFALSE; \r
681       if(!vtx) return kFALSE;\r
682 \r
683       Double_t dz[2],cov[3];\r
684       //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();\r
685       //if(!tpcInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; \r
686       //if(!tpcInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; \r
687       if(!tpcInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; \r
688 \r
689 \r
690       Double_t covar[6]; vtx->GetCovMatrix(covar);\r
691       Double_t p[2]={tpcInnerC->GetParameter()[0]-dz[0],tpcInnerC->GetParameter()[1]-dz[1]};\r
692       Double_t c[3]={covar[2],0.,covar[5]};\r
693       Double_t chi2C=tpcInnerC->GetPredictedChi2(p,c);\r
694       if (chi2C>kVeryBig) return kFALSE; \r
695 \r
696       if(!tpcInnerC->Update(p,c)) return kFALSE;\r
697 \r
698   return kTRUE;\r
699 }\r
700 \r
701 //_____________________________________________________________________________\r
702 Bool_t AlidNdPtTrackDumpTask::ConstrainTrackInner(AliExternalTrackParam *const trackInnerC, const AliESDVertex* vtx, Double_t mass, Double_t b[3])\r
703 {\r
704  // Constrain TPC inner params constrained\r
705  //\r
706       if(!trackInnerC) return kFALSE; \r
707       if(!vtx) return kFALSE;\r
708 \r
709       const Double_t kRadius  = 2.8; \r
710       const Double_t kMaxStep = 1.0; \r
711 \r
712       Double_t dz[2],cov[3];\r
713 \r
714       //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();\r
715       //if(!trackInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; \r
716       //if(!trackInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; \r
717 \r
718       if(!AliTracker::PropagateTrackToBxByBz(trackInnerC,kRadius,mass,kMaxStep,kFALSE)) return kFALSE;\r
719       if(!trackInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; \r
720 \r
721       //\r
722       Double_t covar[6]; vtx->GetCovMatrix(covar);\r
723       Double_t p[2]={trackInnerC->GetParameter()[0]-dz[0],trackInnerC->GetParameter()[1]-dz[1]};\r
724       Double_t c[3]={covar[2],0.,covar[5]};\r
725       Double_t chi2C=trackInnerC->GetPredictedChi2(p,c);\r
726       if (chi2C>kVeryBig) return kFALSE; \r
727 \r
728       if(!trackInnerC->Update(p,c)) return kFALSE;\r
729 \r
730   return kTRUE;\r
731 }\r
732 \r
733 \r
734 //_____________________________________________________________________________\r
735 TParticle *AlidNdPtTrackDumpTask::GetMother(TParticle *const particle, AliStack *const stack) \r
736 {\r
737   if(!particle) return NULL;\r
738   if(!stack) return NULL;\r
739 \r
740   Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
741   TParticle* mother = NULL; \r
742   mother = stack->Particle(motherLabel); \r
743 \r
744 return mother;\r
745 }\r
746 \r
747 //_____________________________________________________________________________\r
748 Bool_t AlidNdPtTrackDumpTask::IsFromConversion(const Int_t label, AliStack *const stack) \r
749 {\r
750   Bool_t isFromConversion = kFALSE;\r
751 \r
752   if(stack) {\r
753     TParticle* particle = stack->Particle(label);\r
754 \r
755     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
756     {\r
757        Int_t mech = particle->GetUniqueID(); // production mechanism \r
758        Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
759 \r
760        Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
761        TParticle* mother = stack->Particle(motherLabel); \r
762        if(mother) {\r
763           Int_t motherPdg = mother->GetPdgCode();\r
764 \r
765           if(!isPrim && mech==5 && motherPdg==kGamma) { \r
766             isFromConversion=kTRUE; \r
767           }\r
768        }\r
769     } \r
770   } \r
771 \r
772   return isFromConversion;\r
773 }\r
774 \r
775 //_____________________________________________________________________________\r
776 Bool_t AlidNdPtTrackDumpTask::IsFromMaterial(const Int_t label, AliStack *const stack) \r
777 {\r
778   Bool_t isFromMaterial = kFALSE;\r
779 \r
780   if(stack) {\r
781     TParticle* particle = stack->Particle(label);\r
782 \r
783     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
784     {\r
785        Int_t mech = particle->GetUniqueID(); // production mechanism \r
786        Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
787 \r
788        Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
789        TParticle* mother = stack->Particle(motherLabel); \r
790        if(mother) {\r
791           if(!isPrim && mech==13) { \r
792             isFromMaterial=kTRUE; \r
793           }\r
794        }\r
795      } \r
796   } \r
797 \r
798   return isFromMaterial;\r
799 }\r
800 \r
801 //_____________________________________________________________________________\r
802 Bool_t AlidNdPtTrackDumpTask::IsFromStrangeness(const Int_t label, AliStack *const stack) \r
803 {\r
804   Bool_t isFromStrangeness = kFALSE;\r
805 \r
806   if(stack) {\r
807     TParticle* particle = stack->Particle(label);\r
808 \r
809     if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
810     {\r
811        Int_t mech = particle->GetUniqueID(); // production mechanism \r
812        Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
813 \r
814        Int_t motherLabel = TMath::Abs(particle->GetMother(0));  \r
815        TParticle* mother = stack->Particle(motherLabel); \r
816        if(mother) {\r
817           Int_t motherPdg = mother->GetPdgCode();\r
818 \r
819           // K+-, lambda, antilambda, K0s decays\r
820           if(!isPrim && mech==4 && \r
821               (TMath::Abs(motherPdg)==kKPlus || TMath::Abs(motherPdg)==kLambda0 || motherPdg==kK0Short))\r
822           {\r
823             isFromStrangeness = kTRUE;\r
824           } \r
825        }\r
826     } \r
827   } \r
828 \r
829   return isFromStrangeness;\r
830 }\r
831 \r
832 \r
833 //_____________________________________________________________________________\r
834 void AlidNdPtTrackDumpTask::FinishTaskOutput() \r
835 {\r
836   //\r
837   // Called one at the end \r
838   // locally on working node\r
839   //\r
840 \r
841   // Post output data.\r
842   PostData(1, fOutput);\r
843   //PostData(0, fOutputSummary);\r
844 }\r
845 \r
846 //_____________________________________________________________________________\r
847 void AlidNdPtTrackDumpTask::Terminate(Option_t *) \r
848 {\r
849   // Called one at the end \r
850   \r
851   // check output data\r
852   fOutputSummary = dynamic_cast<TTree*> (GetOutputData(0));\r
853   if(fOutputSummary) delete fOutputSummary; fOutputSummary=0;\r
854   if(fTreeSRedirector)  delete fTreeSRedirector; fTreeSRedirector=0;\r
855 \r
856   TChain* chain = new TChain("dNdPtTree");\r
857   if(!chain) return;\r
858   chain->Add("dNdPtOutliersAnalysisPbPb.root");\r
859   TTree *tree = chain->CopyTree("1");\r
860   if (chain) { delete chain; chain=0; }\r
861   if(!tree) return;\r
862   tree->Print();\r
863   fOutputSummary = tree;\r
864 \r
865   if (!fOutputSummary) {\r
866     Printf("ERROR: AlidNdPtTrackDumpTask::Terminate(): Output data not avaiable GetOutputData(0)==0x0 ..." );\r
867     return;\r
868   }\r
869   \r
870 \r
871 \r
872   PostData(0, fOutputSummary);\r
873   //PostData(1, fOutput);\r
874 }\r