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