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