ab201e1566888170741954ba51279c5d33643097
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / train / AliAnalysisTaskHighPtDeDx.cxx
1 /*
2   30/03/2012
3   Update: the bug on V0's was fixed, also an option to save VZERO cells was added
4   For the AOD analysis the filter 32 was replaced by the proper golden: 1024
5   // 1024 1<<10                        
6   trackFilter->AddCuts(esdTrackCutsH2); // add r_aa cuts
7
8   30/05/2012
9   Include neutral particles in MC truth
10   Add number of tpc shared clusters
11
12   10/01/13. 
13   Add a flag to see if the particle comes from material or WD
14   
15   20/02/13
16   Correct the bug on the multiplicity estimator. Before it was filled two times, one for the v0ana and other for the track ana. SPD vertex
17
18   10/04/13
19   Store all events with different vertex status.
20   Now the vertex is from TPC if not from SPD. Similar as in TOF
21
22   16/04/13
23   Add quality check to SPD vertex. Status -2
24
25   23/04/13
26   Store MC truth, |eta_lab|<2.4, also store y
27
28
29 */
30
31 #include "AliAnalysisTaskHighPtDeDx.h"
32 #include <bitset>
33 // ROOT includes
34 #include <TList.h>
35 #include <TTree.h>
36 #include <TMath.h>
37 #include <TH1.h>
38 #include <TParticle.h>
39 #include <TFile.h>
40
41 // AliRoot includes
42 #include <AliAnalysisManager.h>
43 #include <AliAnalysisFilter.h>
44 #include <AliESDInputHandler.h>
45 #include <AliESDEvent.h>
46 #include <AliESDVertex.h>
47 #include <AliLog.h>
48 #include <AliExternalTrackParam.h>
49 #include <AliESDtrackCuts.h>
50 #include <AliESDVZERO.h>
51 #include <AliAODVZERO.h>
52
53 #include <AliMCEventHandler.h>
54 #include <AliMCEvent.h>
55 #include <AliStack.h>
56
57 #include <TTreeStream.h>
58
59 #include <AliHeader.h>
60 #include <AliGenPythiaEventHeader.h>
61 #include <AliGenDPMjetEventHeader.h>
62
63 #include "AliCentrality.h" 
64 #include <AliESDv0.h>
65 #include <AliKFVertex.h>
66 #include <AliAODVertex.h>
67
68 #include <AliAODTrack.h> 
69 #include <AliAODPid.h> 
70 #include <AliAODMCHeader.h> 
71 #include <iostream>
72
73
74 // STL includes
75 #include <iostream>
76 using namespace std;
77
78
79 //
80 // Responsible:
81 // Alexandru Dobrin (Wayne State) 
82 // Peter Christiansen (Lund)
83 //
84
85 /* 
86 To do:
87
88 Separate the code into two
89
90 */
91
92
93
94
95 ClassImp(AliAnalysisTaskHighPtDeDx)
96
97 const Double_t AliAnalysisTaskHighPtDeDx::fgkClight = 2.99792458e-2;
98
99 //_____________________________________________________________________________
100 AliAnalysisTaskHighPtDeDx::AliAnalysisTaskHighPtDeDx():
101   AliAnalysisTaskSE(),
102   fESD(0x0),
103   fAOD(0x0),
104   fMC(0x0),
105   fMCStack(0x0),
106   fMCArray(0x0),
107   fTrackFilter(0x0),
108   fTrackFilterGolden(0x0),
109   fTrackFilterTPC(0x0),
110   fV0ArrayGlobalPar(0x0), //
111   fV0ArrayTPCPar(0x0),//
112   fAnalysisType("ESD"),
113   fAnalysisMC(kFALSE),
114   fAnalysisPbPb(kFALSE),
115   fVZEROBranch(kFALSE),
116   fTPCBranch(kFALSE),
117   ftrigBit1(0x0),
118   ftrigBit2(0x0),
119   fRandom(0x0),
120   fEvent(0x0),
121   fTrackArrayGlobalPar(0x0),
122   fTrackArrayTPCPar(0x0),
123   fTrackArrayMC(0x0),
124   fVZEROArray(0x0),
125   fVtxCut(10.0),  
126   fEtaCut(0.9),  
127   fEtaCutStack(1.2),  
128   fMinPt(0.1),
129   fMinPtV0(0.1),
130   fMinCent(0.0),
131   fMaxCent(100.0),
132   fMassCut(0.1),//
133   fLowPtFraction(0.01),
134   fTreeOption(0),
135   fRequireRecV0(kTRUE), //
136   fStoreMcIn(kFALSE),//
137   fMcProcessType(-999),
138   fTriggeredEventMB(-999),
139   fVtxStatus(-999),
140   fZvtx(-999),
141   fZvtxMC(-999),
142   fRun(-999),
143   fEventId(-999),
144   fListOfObjects(0), 
145   fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),
146   fTree(0x0),
147   fn1(0),
148   fn2(0),
149   fcent(0)
150
151
152 {
153   // Default constructor (should not be used)
154
155   //  fRandom = new TRandom(0); // 0 means random seed
156 }
157
158 //______________________________________________________________________________
159 AliAnalysisTaskHighPtDeDx::AliAnalysisTaskHighPtDeDx(const char *name):
160   AliAnalysisTaskSE(name),
161   fESD(0x0),
162   fAOD(0x0),
163   fMC(0x0),
164   fMCStack(0x0),
165   fMCArray(0x0),
166   fTrackFilter(0x0),
167   fTrackFilterGolden(0x0),
168   fTrackFilterTPC(0x0),
169   fV0ArrayGlobalPar(0x0), //
170   fV0ArrayTPCPar(0x0),//
171   fAnalysisType("ESD"),
172   fAnalysisMC(kFALSE),
173   fAnalysisPbPb(kFALSE),
174   fVZEROBranch(kFALSE),
175   fTPCBranch(kFALSE),
176   ftrigBit1(0x0),
177   ftrigBit2(0x0),
178   fRandom(0x0),
179   fEvent(0x0),
180   fTrackArrayGlobalPar(0x0),
181   fTrackArrayTPCPar(0x0),
182   fTrackArrayMC(0x0),
183   fVZEROArray(0x0),
184   fVtxCut(10.0),  
185   fEtaCut(0.9),
186   fEtaCutStack(1.2),    
187   fMinPt(0.1),
188   fMinPtV0(0.1),
189   fMinCent(0.0),
190   fMaxCent(100.0),
191   fMassCut(0.1),//
192   fLowPtFraction(0.01),
193   fTreeOption(0),
194   fRequireRecV0(kTRUE), //
195   fStoreMcIn(kFALSE),//
196   fMcProcessType(-999),
197   fTriggeredEventMB(-999),
198   fVtxStatus(-999),
199   fZvtx(-999),
200   fZvtxMC(-999),
201   fRun(-999),
202   fEventId(-999),
203   fListOfObjects(0), 
204   fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),
205   fTree(0x0),
206   fn1(0),
207   fn2(0),
208   fcent(0)
209
210 {
211
212   if(fTree)fTree=0;
213   // Output slot #1 writes into a TList
214   DefineOutput(1, TList::Class());
215
216 }
217
218 //_____________________________________________________________________________
219 AliAnalysisTaskHighPtDeDx::~AliAnalysisTaskHighPtDeDx()
220 {
221   // Destructor
222   // histograms are in the output list and deleted when the output
223   // list is deleted by the TSelector dtor
224   if (fListOfObjects && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
225     delete fListOfObjects;
226     fListOfObjects = 0;
227   }
228   if (fRandom) delete fRandom;
229   fRandom=0;
230   
231   // //for proof running; I cannot create tree do to memory limitations -> work with THnSparse 
232   // if (fListOfObjects  && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutputList;
233   
234   
235   
236 }
237
238 //______________________________________________________________________________
239 void AliAnalysisTaskHighPtDeDx::UserCreateOutputObjects()
240
241  
242   // This method is called once per worker node
243   // Here we define the output: histograms and debug tree if requested 
244   // We also create the random generator here so it might get different seeds...
245   fRandom = new TRandom(0); // 0 means random seed
246
247   OpenFile(1);
248   fListOfObjects = new TList();
249   fListOfObjects->SetOwner();
250   
251   //
252   // Histograms
253   //  
254   fEvents = new TH1I("fEvents","Number of analyzed events; Events; Counts", 3, 0, 3);
255   fListOfObjects->Add(fEvents);
256
257   fn1=new TH1F("fn1","fn1",5001,-1,5000);
258   fListOfObjects->Add(fn1);
259
260   fn2=new TH1F("fn2","fn2",5001,-1,5000);
261   fListOfObjects->Add(fn2);
262
263   fcent=new TH1F("fcent","fcent",104,-2,102);
264   fListOfObjects->Add(fcent);
265
266   fVtx = new TH1I("fVtx","Vtx info (0=no, 1=yes); Vtx; Counts", 2, -0.5, 1.5);
267   fListOfObjects->Add(fVtx);
268
269   fVtxBeforeCuts = new TH1F("fVtxBeforeCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
270   fListOfObjects->Add(fVtxBeforeCuts);
271   
272   fVtxAfterCuts = new TH1F("fVtxAfterCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
273   fListOfObjects->Add(fVtxAfterCuts);
274
275   if (fAnalysisMC) {    
276     fVtxMC = new TH1I("fVtxMC","Vtx info - no trigger cut (0=no, 1=yes); Vtx; Counts", 2, -0.5, 1.5);
277     fListOfObjects->Add(fVtxMC);
278   }
279
280   if (fTreeOption) {
281
282     fTree = new TTree("tree","Event data");
283     fEvent = new DeDxEvent();
284     fTree->Branch("event", &fEvent);
285     //array with tracks
286     fTrackArrayGlobalPar = new TClonesArray("DeDxTrack", 1000);
287     fTree->Bronch("trackGlobalPar", "TClonesArray", &fTrackArrayGlobalPar);
288     //array with v0's
289     fV0ArrayGlobalPar = new TClonesArray("DeDxV0", 1000);
290     fTree->Bronch("v0GlobalPar", "TClonesArray", &fV0ArrayGlobalPar);
291
292     if(fTPCBranch){
293       fTrackArrayTPCPar = new TClonesArray("DeDxTrack", 1000);
294       fTree->Bronch("trackTPCPar", "TClonesArray", &fTrackArrayTPCPar);
295       fV0ArrayTPCPar = new TClonesArray("DeDxV0", 1000);
296       fTree->Bronch("v0TPCPar", "TClonesArray", &fV0ArrayTPCPar);
297     }
298     if(fVZEROBranch){
299       fVZEROArray = new TClonesArray("VZEROCell", 1000);
300       fTree->Bronch("cellVZERO", "TClonesArray", &fVZEROArray);
301     }
302
303
304     if (fAnalysisMC && fStoreMcIn) {    
305       fTrackArrayMC = new TClonesArray("DeDxTrackMC", 1000);
306       fTree->Bronch("trackMC", "TClonesArray", &fTrackArrayMC);
307     }
308
309
310
311     fTree->SetDirectory(0);
312
313     fListOfObjects->Add(fTree);
314
315   }
316
317   // Post output data.
318   PostData(1, fListOfObjects);
319 }
320
321 //______________________________________________________________________________
322 void AliAnalysisTaskHighPtDeDx::UserExec(Option_t *) 
323 {
324   // Main loop
325
326   //
327   // First we make sure that we have valid input(s)!
328   //
329
330
331
332   AliVEvent *event = InputEvent();
333   if (!event) {
334      Error("UserExec", "Could not retrieve event");
335      return;
336   }
337
338
339   if (fAnalysisType == "ESD"){
340     fESD = dynamic_cast<AliESDEvent*>(event);
341     if(!fESD){
342       Printf("%s:%d ESDEvent not found in Input Manager",(char*)__FILE__,__LINE__);
343       this->Dump();
344       return;
345     }    
346   } else {
347     fAOD = dynamic_cast<AliAODEvent*>(event);
348     if(!fAOD){
349       Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
350       this->Dump();
351       return;
352     }    
353
354   }
355
356  
357   
358   if (fAnalysisMC) {
359
360     if (fAnalysisType == "ESD"){
361       fMC = dynamic_cast<AliMCEvent*>(MCEvent());
362       if(!fMC){
363         Printf("%s:%d MCEvent not found in Input Manager",(char*)__FILE__,__LINE__);
364         this->Dump();
365         return;
366       }    
367
368       fMCStack = fMC->Stack();
369       
370       if(!fMCStack){
371         Printf("%s:%d MCStack not found in Input Manager",(char*)__FILE__,__LINE__);
372         this->Dump();
373         return;
374       }    
375     } else { // AOD
376
377       fMC = dynamic_cast<AliMCEvent*>(MCEvent());
378       if(fMC)
379         fMC->Dump();
380
381       fMCArray = (TClonesArray*)fAOD->FindListObject("mcparticles");
382       if(!fMCArray){
383         Printf("%s:%d AOD MC array not found in Input Manager",(char*)__FILE__,__LINE__);
384         this->Dump();
385         return;
386       }    
387     }
388   }
389
390   
391   // Get trigger decision
392   fTriggeredEventMB = 0; //init
393
394   // cout << " trigger " << ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler())) ->IsEventSelected() << "cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc"<<endl;
395     if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
396      ->IsEventSelected() & ftrigBit1 ){
397     fn1->Fill(1);
398     fTriggeredEventMB = 1;  //event triggered as minimum bias
399   }
400   if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
401      ->IsEventSelected() & ftrigBit2 ){
402     // From AliVEvent:
403     //    kINT7         = BIT(1), // V0AND trigger, offline V0 selection
404     fTriggeredEventMB += 2;  
405     fn2->Fill(1);
406   }
407  
408   //if(fTriggeredEventMB == 0) fTriggeredEventMB = 1; //ignore trigger //hljunggr
409   //std::cout << "trigger " << fTriggeredEventMB << " ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc" << std::endl;
410   // Get process type for MC
411   fMcProcessType = 0; // -1=invalid, 0=data, 1=ND, 2=SD, 3=DD
412
413   // real data that are not triggered we skip
414   if(!fAnalysisMC && !fTriggeredEventMB)
415     return; 
416
417
418   //   std::cout << "trigger Martin: " << (bitset<20>) ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() << std::endl; print trigger
419
420
421   if (fAnalysisMC) {
422     
423
424
425     if (fAnalysisType == "ESD"){
426
427   
428
429       AliHeader* headerMC = fMC->Header();
430       if (headerMC) {
431         
432         AliGenEventHeader* genHeader = headerMC->GenEventHeader();
433         TArrayF vtxMC(3); // primary vertex  MC 
434         vtxMC[0]=9999; vtxMC[1]=9999;  vtxMC[2]=9999; //initialize with dummy
435         if (genHeader) {
436           genHeader->PrimaryVertex(vtxMC);
437         }
438         fZvtxMC = vtxMC[2];
439         
440         // PYTHIA:
441         AliGenPythiaEventHeader* pythiaGenHeader =
442           dynamic_cast<AliGenPythiaEventHeader*>(headerMC->GenEventHeader());
443         if (pythiaGenHeader) {  //works only for pythia
444           fMcProcessType =  GetPythiaEventProcessType(pythiaGenHeader->ProcessType());
445         }
446         // PHOJET:
447         AliGenDPMjetEventHeader* dpmJetGenHeader =
448           dynamic_cast<AliGenDPMjetEventHeader*>(headerMC->GenEventHeader());
449         if (dpmJetGenHeader) {
450           fMcProcessType = GetDPMjetEventProcessType(dpmJetGenHeader->ProcessType());
451         }
452       }
453     } else { // AOD
454       
455   
456
457       AliAODMCHeader* mcHeader = dynamic_cast<AliAODMCHeader*>(fAOD->FindListObject("mcHeader")); 
458
459
460       if(mcHeader) {
461         fZvtxMC = mcHeader->GetVtxZ();
462         
463
464
465         if(strstr(mcHeader->GetGeneratorName(), "Pythia")) {
466           fMcProcessType =  GetPythiaEventProcessType(mcHeader->GetEventType());
467         } else {
468           fMcProcessType =  GetDPMjetEventProcessType(mcHeader->GetEventType());
469         }
470       }
471     }
472   }
473   
474   // There are 3 cases
475   // Vertex: NO  - status -1
476   // SDP Vertex wit bad quality: status -2
477   // Vertex: YES : outside cut - status 0
478   //             : inside cut  - status 1
479   // We have to be careful how we normalize because we probably want to
480   // normalize to:
481   // Nevents=(No vertex + outside + inside)/(out + in)*in
482   
483
484   if (fAnalysisType == "ESD"){
485     
486     const AliESDVertex *vtxESD = fESD->GetPrimaryVertexTracks();
487     if(vtxESD->GetNContributors()<1) {
488       // SPD vertex
489       vtxESD = fESD->GetPrimaryVertexSPD();
490       /* quality checks on SPD-vertex */
491       TString vertexType = vtxESD->GetTitle();
492       if (vertexType.Contains("vertexer: Z") && (vtxESD->GetDispersion() > 0.04 || vtxESD->GetZRes() > 0.25))  
493         fZvtx  = -1599; //vertex = 0x0; //
494       else if (vtxESD->GetNContributors()<1) 
495         fZvtx  = -999; //vertex = 0x0; //
496       else
497         fZvtx = vtxESD->GetZ();
498     }  
499     else
500       fZvtx = vtxESD->GetZ();
501     
502   }
503   else // AOD
504     fZvtx = GetVertex(fAOD);
505     
506   fVtxStatus = -999;
507   
508
509
510   if(fZvtx<-1500) {
511     
512     fVtxStatus = -2;
513  
514   } 
515   else if(fZvtx<-990&&fZvtx>-1500) {
516     
517     fVtxStatus = -1;
518     if(fTriggeredEventMB)
519       fVtx->Fill(0);
520     if(fAnalysisMC)
521       fVtxMC->Fill(0);
522   } else {
523     
524     if(fTriggeredEventMB)
525       fVtx->Fill(1);
526     if(fAnalysisMC)
527       fVtxMC->Fill(1);
528     fVtxBeforeCuts->Fill(fZvtx);
529     fVtxStatus = 0;
530     if (TMath::Abs(fZvtx) < fVtxCut) {  
531       fVtxAfterCuts->Fill(fZvtx);
532       fVtxStatus = 1;
533     }
534   }
535   
536   // store MC event data no matter what
537   //if(fAnalysisMC && fStoreMcIn) {
538   if(fAnalysisMC) {
539
540     if (fAnalysisType == "ESD") {
541       ProcessMCTruthESD();
542     } else { // AOD
543
544       ProcessMCTruthAOD();
545   
546     }
547   }      
548   
549
550
551   Float_t centralityV0M = -10;
552   Float_t centralityV0A = -10;
553   Float_t centralityZNA = -10;
554   Float_t centralityCL1 = -10;
555
556   // only analyze triggered events
557   //   cout << " fTriggeredEventMB " << fTriggeredEventMB << "fAnalysisType " << fAnalysisType << " fAnalysisPbPb " <<fAnalysisPbPb << endl;
558   if(fTriggeredEventMB) {
559     //cout << "inside triggered !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! " << endl;
560     if (fAnalysisType == "ESD"){
561       if(fAnalysisPbPb){
562         AliCentrality *centObject = fESD->GetCentrality();
563         centralityV0M = centObject->GetCentralityPercentile("V0M");
564         centralityV0A = centObject->GetCentralityPercentile("V0A");
565         centralityZNA = centObject->GetCentralityPercentile("ZNA");
566         centralityCL1 = centObject->GetCentralityPercentile("CL1");
567      
568         //if(  !(fMinCent>=20 && fMaxCent <=40) ) return; //take out //hljunggr
569         if((centralityV0M>fMaxCent)||(centralityV0M<fMinCent))return; //hljunggr comment out, put back in!!
570       }
571       fcent->Fill(centralityV0A);
572
573       AnalyzeESD(fESD);
574     } else { // AOD
575       if(fAnalysisPbPb){
576         AliCentrality *centObject = fAOD->GetCentrality();
577         if(centObject){
578           centralityV0M = centObject->GetCentralityPercentile("V0M"); 
579           centralityV0A = centObject->GetCentralityPercentile("V0A");
580           centralityZNA = centObject->GetCentralityPercentile("ZNA");
581           centralityCL1 = centObject->GetCentralityPercentile("CL1");
582         }
583
584         //if(  !(fMinCent>=20 && fMaxCent <=40) ) return; //take out //hljunggr
585         if((centralityV0M>fMaxCent)||(centralityV0M<fMinCent))return; //hljunggr comment out, put back in!!
586       }
587       fcent->Fill(centralityV0M);
588       AnalyzeAOD(fAOD);
589     }
590   }
591
592   if( fTreeOption) {
593     
594     fEvent->process = fMcProcessType;
595     fEvent->trig    = fTriggeredEventMB;
596     fEvent->zvtxMC  = fZvtxMC;
597     fEvent->cent      = centralityV0M;
598     fEvent->centV0A      = centralityV0A;
599     fEvent->centZNA      = centralityZNA;
600     fEvent->centCL1      = centralityCL1;
601
602
603     fTree->Fill();
604     if(fTrackArrayGlobalPar)fTrackArrayGlobalPar->Clear();
605     if(fV0ArrayGlobalPar)fV0ArrayGlobalPar->Clear();
606
607     if(fTPCBranch){
608       if(fTrackArrayTPCPar)fTrackArrayTPCPar->Clear();
609       if(fV0ArrayTPCPar)fV0ArrayTPCPar->Clear();
610     }
611     if(fVZEROArray)
612       fVZEROArray->Clear();
613
614     if (fAnalysisMC)    
615       fTrackArrayMC->Clear();
616       
617
618   }
619   
620   // Post output data.
621   PostData(1, fListOfObjects);
622 }
623
624 //________________________________________________________________________
625 void AliAnalysisTaskHighPtDeDx::AnalyzeESD(AliESDEvent* esdEvent)
626 {
627   fRun  = esdEvent->GetRunNumber();
628   fEventId = 0;
629   if(esdEvent->GetHeader())
630     fEventId = GetEventIdAsLong(esdEvent->GetHeader());
631   
632   Short_t isPileup = esdEvent->IsPileupFromSPD();
633   /*
634   cout"  nTPCtracks="<<esdEvent->nTPCtracks<<"  nITStracks="<<esdEvent->nITStracks<<endl;
635   if(esdEvent->nTPCtracks==0 && esdEvent->nITStracks==0){
636     cout<<"Evento defectuoso!!"<<endl;
637     }*/
638
639
640   //  Int_t     event     = esdEvent->GetEventNumberInFile();
641   UInt_t    time      = esdEvent->GetTimeStamp();
642   //  ULong64_t trigger   = esdEvent->GetTriggerMask();
643   Float_t   magf      = esdEvent->GetMagneticField();
644
645
646
647
648
649   if(fTriggeredEventMB) {// Only MC case can we have not triggered events
650     
651     // accepted event
652     fEvents->Fill(0);
653     
654     
655     //Change, 10/04/13. Now accept all events to do a correct normalization
656     //if(fVtxStatus!=1) return; // accepted vertex
657     Int_t nESDTracks = esdEvent->GetNumberOfTracks();
658     
659     ProduceArrayTrksESD( esdEvent, kGlobalTrk );//produce array with global track parameters
660     ProduceArrayV0ESD( esdEvent, kGlobalTrk );//v0's
661
662     if(fTPCBranch){
663       ProduceArrayTrksESD( esdEvent, kTPCTrk );//array with tpc track parametes
664       ProduceArrayV0ESD( esdEvent, kTPCTrk );
665     }
666     
667     fEvents->Fill(1);
668     
669     if(fVZEROBranch){
670       AliESDVZERO *esdV0 = esdEvent->GetVZEROData();// loop sobre canales del V0 para obtener las multiplicidad
671       for (Int_t iCh=0; iCh<64; ++iCh) { 
672         Float_t multv=esdV0->GetMultiplicity(iCh); 
673         Int_t intexv=iCh;
674         VZEROCell* cellv0 = new((*fVZEROArray)[iCh]) VZEROCell();
675         cellv0->cellindex=intexv;
676         cellv0->cellmult= multv;
677       }   
678     }
679
680
681
682   } // end if triggered
683   
684   if(fTreeOption) {
685
686     fEvent->run       = fRun;
687     fEvent->eventid   = fEventId;
688     fEvent->time      = time;
689     //fEvent->cent      = centrality;
690     fEvent->mag       = magf;
691     fEvent->zvtx      = fZvtx;
692     fEvent->vtxstatus = fVtxStatus;
693     fEvent->pileup    = isPileup;
694
695   }
696
697
698
699
700 }
701
702 //________________________________________________________________________
703 void AliAnalysisTaskHighPtDeDx::AnalyzeAOD(AliAODEvent* aodEvent)
704 {
705   fRun  = aodEvent->GetRunNumber();
706   fEventId = 0;
707   if(aodEvent->GetHeader())
708     fEventId = GetEventIdAsLong(aodEvent->GetHeader());
709    
710   UInt_t    time      = 0; // Missing AOD info? aodEvent->GetTimeStamp();
711   Float_t   magf      = aodEvent->GetMagneticField();
712
713   //Int_t     trackmult = 0; // no pt cuts
714   //Int_t     nadded    = 0;
715   Short_t isPileup = aodEvent->IsPileupFromSPD();
716
717   
718
719   if(fTriggeredEventMB) {// Only MC case can we have not triggered events
720     
721     // accepted event
722     fEvents->Fill(0);
723     
724     //if(fVtxStatus!=1) return; // accepted vertex
725     //Int_t nAODTracks = aodEvent->GetNumberOfTracks();  
726
727     ProduceArrayTrksAOD( aodEvent, kGlobalTrk );
728
729     ProduceArrayV0AOD( aodEvent, kGlobalTrk );
730     if(fTPCBranch){
731       ProduceArrayTrksAOD( aodEvent, kTPCTrk );
732       ProduceArrayV0AOD( aodEvent, kTPCTrk );
733     }
734     fEvents->Fill(1);
735
736     if(fVZEROBranch){
737       AliAODVZERO *esdV0 = aodEvent->GetVZEROData();// loop sobre canales del V0 para obtener las multiplicidad
738       for (Int_t iCh=0; iCh<64; ++iCh) { 
739         Float_t multv=esdV0->GetMultiplicity(iCh); 
740         Int_t intexv=iCh;
741         VZEROCell* cellv0 = new((*fVZEROArray)[iCh]) VZEROCell();
742         cellv0->cellindex=intexv;
743         cellv0->cellmult= multv;
744       }   
745     }
746
747
748
749   } // end if triggered
750   
751   if(fTreeOption) {
752
753     //Sort(fTrackArrayGlobalPar, kFALSE);
754
755     fEvent->run       = fRun;
756     fEvent->eventid   = fEventId;
757     fEvent->time      = time;
758     //fEvent->cent      = centrality;
759     fEvent->mag       = magf;
760     fEvent->zvtx      = fZvtx;
761     fEvent->vtxstatus = fVtxStatus;
762     //fEvent->trackmult = trackmult;
763     //fEvent->n         = nadded;
764     fEvent->pileup    = isPileup;
765   }
766 }
767
768 //_____________________________________________________________________________
769 Float_t AliAnalysisTaskHighPtDeDx::GetVertex(const AliVEvent* event) const
770 {
771   Float_t zvtx = -999;
772   
773   const AliVVertex* primaryVertex = event->GetPrimaryVertex(); 
774   
775   if(primaryVertex->GetNContributors()>0)
776     zvtx = primaryVertex->GetZ();
777
778   return zvtx;
779 }
780
781 //_____________________________________________________________________________
782 Short_t AliAnalysisTaskHighPtDeDx::GetPidCode(Int_t pdgCode) const 
783 {
784   // return our internal code for pions, kaons, and protons
785   
786   Short_t pidCode = 6;
787   
788   switch (TMath::Abs(pdgCode)) {
789   case 211:
790     pidCode = 1; // pion
791     break;
792   case 321:
793     pidCode = 2; // kaon
794     break;
795   case 2212:
796     pidCode = 3; // proton
797     break;
798   case 11:
799     pidCode = 4; // electron
800     break;
801   case 13:
802     pidCode = 5; // muon
803     break;
804   default:
805     pidCode = 6;  // something else?
806   };
807   
808   return pidCode;
809 }
810
811
812 //_____________________________________________________________________________
813 void AliAnalysisTaskHighPtDeDx::ProcessMCTruthESD() 
814 {
815   // Fill the special MC histogram with the MC truth info
816   
817  
818   //cout << "mc truth " << " fAnalysisMC " << fAnalysisMC << " fstoremcin " << fStoreMcIn << endl;
819   Short_t trackmult = 0;
820   Short_t nadded    = 0;
821   const Int_t nTracksMC = fMCStack->GetNtrack();
822
823   Float_t  sphericityMC=GetSphericityTrue(fMCStack, 0.8, 0.5);
824   Float_t  spherocityMC=GetSpherocityTrue(fMCStack, 0.8, 0.5);
825
826
827   for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
828     
829    
830
831     //Cuts
832     if(!(fMCStack->IsPhysicalPrimary(iTracks)))
833       continue;
834     
835     TParticle* trackMC = fMCStack->Particle(iTracks);
836     
837     TParticlePDG* pdgPart = trackMC ->GetPDG();
838     Double_t chargeMC = pdgPart->Charge();
839
840     
841     if (TMath::Abs(trackMC->Eta()) > fEtaCutStack )
842       continue;
843     
844     trackmult++;
845     
846     //   "charge:pt:p:eta:phi:pidCode"
847     Float_t ptMC      = trackMC->Pt();
848     Float_t pMC       = trackMC->P();
849     Float_t etaMC     = trackMC->Eta();
850     Float_t phiMC     = trackMC->Phi();
851     Float_t yMC       = trackMC->Y();
852
853     Int_t pdgCode = trackMC->GetPdgCode();
854     Short_t pidCodeMC = 0;
855     pidCodeMC = GetPidCode(pdgCode);
856     
857     // Here we want to add some of the MC histograms!
858     
859     Bool_t lIsStrangeness = kFALSE; 
860     if ( TMath::Abs(pdgCode)==310 || TMath::Abs(pdgCode)==3122 || TMath::Abs(pdgCode)==3312 || TMath::Abs(pdgCode)==3334 ) lIsStrangeness = kTRUE; 
861     
862     // And therefore we first cut here!
863     if (trackMC->Pt() < fMinPt && !lIsStrangeness) {
864       // Keep small fraction of low pT tracks
865       if(fRandom->Rndm() > fLowPtFraction)      continue; 
866     } // else {
867     // Here we want to add the high pt part of the MC histograms!
868     //    }
869     
870     // And therefore we first cut here!
871     if (trackMC->Pt() < fMinPtV0 && lIsStrangeness) {
872       // Keep small fraction of low pT tracks
873       if(fRandom->Rndm() > fLowPtFraction)      continue; 
874     } // else {
875     // Here we want to add the high pt part of the MC histograms!
876     //    }
877     
878     if(fTreeOption) {
879       
880       DeDxTrackMC* track = new((*fTrackArrayMC)[nadded]) DeDxTrackMC();
881       nadded++;
882       
883       track->pMC   = pMC;
884       track->ptMC  = ptMC;
885       track->etaMC = etaMC;
886       track->phiMC = phiMC;
887       track->yMC   = yMC;
888       track->qMC   = Short_t(chargeMC);
889       track->pidMC = pidCodeMC;
890       track->pdgMC = pdgCode;
891     }
892     
893   }//MC track loop
894   
895   if(fTreeOption) {
896     
897     Sort(fTrackArrayMC, kTRUE);
898
899     fEvent->trackmultMC = trackmult;
900     fEvent->nMC         = nadded;
901     fEvent->sphericityMC         = sphericityMC;
902     fEvent->spherocityMC         = spherocityMC;
903
904   }
905   
906 }
907
908 //_____________________________________________________________________________
909 void AliAnalysisTaskHighPtDeDx::ProcessMCTruthAOD() 
910 {
911   // Fill the special MC histogram with the MC truth info
912   //cout << " debug 1 " << endl;
913   Short_t trackmult = 0;
914   Short_t nadded    = 0;
915   const Int_t nTracksMC = fMCArray->GetEntriesFast();
916   //cout << " debug 2 " << endl;
917
918   //cout << " debug 3 " << endl;
919   for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
920     //cout << " debug 31 " << endl;
921     AliAODMCParticle* trackMC = dynamic_cast<AliAODMCParticle*>(fMCArray->At(iTracks));
922     //cout << " debug 32 " << endl;
923     //Cuts
924     if(!(trackMC->IsPhysicalPrimary()))
925       continue;
926     //cout << " debug 33 " << endl;
927     Double_t chargeMC = trackMC->Charge();
928
929     
930     if (TMath::Abs(trackMC->Eta()) > fEtaCutStack )
931       continue;
932     //cout << " debug 34 " << endl;
933     trackmult++;
934     
935     //   "charge:pt:p:eta:phi:pidCode"
936     Float_t ptMC      = trackMC->Pt();
937     Float_t pMC       = trackMC->P();
938     Float_t etaMC     = trackMC->Eta();
939     Float_t phiMC     = trackMC->Phi();
940     Float_t yMC       = trackMC->Y();
941     //cout << " debug 35 " << endl;
942     Int_t pdgCode = trackMC->PdgCode();
943     Short_t pidCodeMC = 0;
944     pidCodeMC = GetPidCode(pdgCode);
945    
946
947  
948  
949     // Here we want to add some of the MC histograms!
950         
951     Bool_t lIsStrangeness = kFALSE; 
952     if ( TMath::Abs(pdgCode)==310 || TMath::Abs(pdgCode)==3122 || TMath::Abs(pdgCode)==3312 || TMath::Abs(pdgCode)==3334 ) lIsStrangeness = kTRUE; 
953
954     // And therefore we first cut here!
955     if (trackMC->Pt() < fMinPt && !lIsStrangeness) {
956       // Keep small fraction of low pT tracks
957       if(fRandom->Rndm() > fLowPtFraction)      continue; 
958     } // else {
959     // Here we want to add the high pt part of the MC histograms!
960     //    }
961     
962     // And therefore we first cut here!
963     if (trackMC->Pt() < fMinPtV0 && lIsStrangeness) {
964       // Keep small fraction of low pT tracks
965       if(fRandom->Rndm() > fLowPtFraction)      continue; 
966     } // else {
967     // Here we want to add the high pt part of the MC histograms!
968     //    }
969     
970     //cout << " debug 36 " << endl;
971  
972     //cout << "fTreeOption " << fTreeOption << endl;
973     if(fTreeOption) {
974       
975
976       //cout << " debug 361, nadded =  " << nadded << " fTrackArrayMC" << fTrackArrayMC <<  endl;
977       DeDxTrackMC* track = new((*fTrackArrayMC)[nadded]) DeDxTrackMC();
978       //cout << " debug 3611, nadded =  " << nadded <<  endl;
979       nadded++;
980       //cout << " debug 3612, nadded =  " << nadded <<  endl;
981       //cout << " debug 362 " << endl;
982       track->pMC   = pMC;
983       //cout << " debug 363 " << endl;
984       track->ptMC  = ptMC;
985       track->etaMC = etaMC;
986       track->phiMC = phiMC;
987       track->yMC   = yMC;
988       track->qMC   = Short_t(chargeMC);
989       track->pidMC = pidCodeMC;
990       track->pdgMC = pdgCode;
991       //cout << " debug 364 " << endl;
992
993     }
994     //cout << " debug 37 " << endl;
995   }//MC track loop
996   //cout << " debug 4 " << endl;
997   if(fTreeOption) {
998     
999     Sort(fTrackArrayMC, kTRUE);
1000
1001     fEvent->trackmultMC = trackmult;
1002     fEvent->nMC         = nadded;
1003     fEvent->sphericityMC         = 0;
1004     fEvent->spherocityMC         = 0;
1005
1006   }
1007   //cout << " debug 5 " << endl;
1008 }
1009
1010 //_____________________________________________________________________________
1011 Short_t AliAnalysisTaskHighPtDeDx::GetPythiaEventProcessType(Int_t pythiaType) {
1012   //
1013   // Get the process type of the event.  PYTHIA
1014   //
1015   // source PWG0   dNdpt 
1016
1017   Short_t globalType = -1; //init
1018       
1019   if(pythiaType==92||pythiaType==93){
1020     globalType = 2; //single diffractive
1021   }
1022   else if(pythiaType==94){
1023     globalType = 3; //double diffractive
1024   }
1025   //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB??
1026   else {
1027     globalType = 1;  //non diffractive
1028   }
1029   return globalType;
1030 }
1031
1032 //_____________________________________________________________________________
1033 Short_t AliAnalysisTaskHighPtDeDx::GetDPMjetEventProcessType(Int_t dpmJetType) {
1034   //
1035   // get the process type of the event.  PHOJET
1036   //
1037   //source PWG0   dNdpt 
1038   // can only read pythia headers, either directly or from cocktalil header
1039   Short_t globalType = -1;
1040   
1041   if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction
1042     globalType = 1;
1043   }
1044   else if (dpmJetType==5 || dpmJetType==6) {
1045     globalType = 2;
1046   }
1047   else if (dpmJetType==7) {
1048     globalType = 3;
1049   }
1050   return globalType;
1051 }
1052
1053 //_____________________________________________________________________________
1054 ULong64_t AliAnalysisTaskHighPtDeDx::GetEventIdAsLong(AliVHeader* header) const
1055 {
1056   // To have a unique id for each event in a run!
1057   // Modified from AliRawReader.h
1058   return ((ULong64_t)header->GetBunchCrossNumber()+
1059           (ULong64_t)header->GetOrbitNumber()*3564+
1060           (ULong64_t)header->GetPeriodNumber()*16777215*3564);
1061 }
1062
1063
1064 //____________________________________________________________________
1065 TParticle* AliAnalysisTaskHighPtDeDx::FindPrimaryMother(AliStack* stack, Int_t label)
1066 {
1067   //
1068   // Finds the first mother among the primary particles of the particle identified by <label>,
1069   // i.e. the primary that "caused" this particle
1070   //
1071   // Taken from AliPWG0Helper class
1072   //
1073
1074   Int_t motherLabel = FindPrimaryMotherLabel(stack, label);
1075   if (motherLabel < 0)
1076     return 0;
1077
1078   return stack->Particle(motherLabel);
1079 }
1080
1081 //____________________________________________________________________
1082 Int_t AliAnalysisTaskHighPtDeDx::FindPrimaryMotherLabel(AliStack* stack, Int_t label)
1083 {
1084   //
1085   // Finds the first mother among the primary particles of the particle identified by <label>,
1086   // i.e. the primary that "caused" this particle
1087   //
1088   // returns its label
1089   //
1090   // Taken from AliPWG0Helper class
1091   //
1092   const Int_t nPrim  = stack->GetNprimary();
1093   
1094   while (label >= nPrim) {
1095
1096     //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
1097
1098     TParticle* particle = stack->Particle(label);
1099     if (!particle) {
1100       
1101       AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label));
1102       return -1;
1103     }
1104     
1105     // find mother
1106     if (particle->GetMother(0) < 0) {
1107
1108       AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
1109       return -1;
1110     }
1111     
1112     label = particle->GetMother(0);
1113   }
1114   
1115   return label;
1116 }
1117
1118 //____________________________________________________________________
1119 AliAODMCParticle* AliAnalysisTaskHighPtDeDx::FindPrimaryMotherAOD(AliAODMCParticle* startParticle)
1120 {
1121   //
1122   // Finds the first mother among the primary particles of the particle identified by <label>,
1123   // i.e. the primary that "caused" this particle
1124   //
1125   // Taken from AliPWG0Helper class
1126   //
1127
1128   AliAODMCParticle* mcPart = startParticle;
1129
1130   while (mcPart)
1131     {
1132       
1133       if(mcPart->IsPrimary())
1134         return mcPart;
1135       
1136       Int_t mother = mcPart->GetMother();
1137
1138       mcPart = dynamic_cast<AliAODMCParticle*>(fMCArray->At(mother));
1139     }
1140
1141   return 0;
1142 }
1143
1144
1145 //V0______________________________________
1146 //____________________________________________________________________
1147 TParticle* AliAnalysisTaskHighPtDeDx::FindPrimaryMotherV0(AliStack* stack, Int_t label)
1148 {
1149   //
1150   // Finds the first mother among the primary particles of the particle identified by <label>,
1151   // i.e. the primary that "caused" this particle
1152   //
1153   // Taken from AliPWG0Helper class
1154   //
1155
1156   Int_t nSteps = 0;
1157
1158   Int_t motherLabel = FindPrimaryMotherLabelV0(stack, label, nSteps);
1159   if (motherLabel < 0)
1160     return 0;
1161
1162   return stack->Particle(motherLabel);
1163 }
1164
1165 //____________________________________________________________________
1166 Int_t AliAnalysisTaskHighPtDeDx::FindPrimaryMotherLabelV0(AliStack* stack, Int_t label, Int_t& nSteps)
1167 {
1168   //
1169   // Finds the first mother among the primary particles of the particle identified by <label>,
1170   // i.e. the primary that "caused" this particle
1171   //
1172   // returns its label
1173   //
1174   // Taken from AliPWG0Helper class
1175   //
1176   nSteps = 0;
1177   const Int_t nPrim  = stack->GetNprimary();
1178   
1179   while (label >= nPrim) {
1180
1181     //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
1182
1183     nSteps++; // 1 level down
1184     
1185     TParticle* particle = stack->Particle(label);
1186     if (!particle) {
1187       
1188       AliDebugGeneral("FindPrimaryMotherLabelV0", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label));
1189       return -1;
1190     }
1191     
1192     // find mother
1193     if (particle->GetMother(0) < 0) {
1194
1195       AliDebugGeneral("FindPrimaryMotherLabelV0", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
1196       return -1;
1197     }
1198     
1199     label = particle->GetMother(0);
1200   }
1201   
1202   return label;
1203 }
1204
1205 //____________________________________________________________________
1206 AliAODMCParticle* AliAnalysisTaskHighPtDeDx::FindPrimaryMotherAODV0(AliAODMCParticle* startParticle, Int_t& nSteps)
1207 {
1208   //
1209   // Finds the first mother among the primary particles of the particle identified by <label>,
1210   // i.e. the primary that "caused" this particle
1211   //
1212   // Taken from AliPWG0Helper class
1213   //
1214
1215   nSteps = 0;
1216
1217   AliAODMCParticle* mcPart = startParticle;
1218
1219   while (mcPart)
1220     {
1221       
1222       if(mcPart->IsPrimary())
1223         return mcPart;
1224       
1225       Int_t mother = mcPart->GetMother();
1226
1227       mcPart = dynamic_cast<AliAODMCParticle*>(fMCArray->At(mother));
1228       nSteps++; // 1 level down
1229     }
1230
1231   return 0;
1232 }
1233
1234
1235
1236
1237  
1238 //_____________________________________________________________________________
1239 void AliAnalysisTaskHighPtDeDx::Sort(TClonesArray* array, Bool_t isMC) 
1240 {
1241   const Int_t n = array->GetEntriesFast(); 
1242   if(n==0) {
1243
1244     if(isMC) 
1245       fEvent->ptmaxMC = 0;
1246     else
1247       fEvent->ptmax   = 0;
1248       
1249     return;
1250   }
1251
1252   Float_t ptArray[n];
1253   Int_t   indexArray[n];
1254   
1255   for(Int_t i = 0; i < n; i++) {
1256
1257     Float_t pt = 0;
1258     if(isMC) {
1259       DeDxTrackMC* track = (DeDxTrackMC*)array->At(i);
1260       pt = track->ptMC;
1261     } else {
1262       DeDxTrack* track = (DeDxTrack*)array->At(i);
1263       pt = track->pt;
1264     }
1265     ptArray[i]    = pt;
1266     indexArray[i] = i;
1267   }
1268
1269   TMath::Sort(n, ptArray, indexArray);
1270   
1271   // set max pt
1272   if(isMC) {
1273     DeDxTrackMC* track = (DeDxTrackMC*)array->At(indexArray[0]);
1274     fEvent->ptmaxMC = track->ptMC;
1275   } else {
1276     DeDxTrack* track = (DeDxTrack*)array->At(indexArray[0]);
1277     fEvent->ptmax   = track->pt;
1278   }
1279   
1280   // set order of each track
1281   for(Int_t i = 0; i < n; i++) {
1282     
1283     if(isMC) {
1284       DeDxTrackMC* track = (DeDxTrackMC*)array->At(indexArray[i]);
1285       track->orderMC = i;
1286     } else {
1287       DeDxTrack* track = (DeDxTrack*)array->At(indexArray[i]);
1288       track->order = i;
1289     }
1290   }
1291 }
1292 //__________________________________________________________________
1293 void AliAnalysisTaskHighPtDeDx::ProduceArrayTrksESD( AliESDEvent *ESDevent, AnalysisMode analysisMode ){
1294   
1295   const Int_t nESDTracks = ESDevent->GetNumberOfTracks();
1296   Int_t trackmult=0;
1297
1298
1299   Int_t nadded=0;
1300   if( analysisMode == kGlobalTrk ){
1301     if(fTrackArrayGlobalPar)
1302       fTrackArrayGlobalPar->Clear();
1303   } else if( analysisMode == kTPCTrk ){
1304     if(fTrackArrayTPCPar)
1305       fTrackArrayTPCPar->Clear();
1306   }
1307
1308   //Fix, for pPb LHC13b pass2 it does not work
1309   
1310   AliESDtrackCuts* esdTrackCutsGolden = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE,1);
1311   esdTrackCutsGolden->SetMinNCrossedRowsTPC(120);
1312   esdTrackCutsGolden->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
1313   esdTrackCutsGolden->SetMaxChi2PerClusterITS(36);
1314   esdTrackCutsGolden->SetMaxFractionSharedTPCClusters(0.4);
1315   esdTrackCutsGolden->SetMaxChi2TPCConstrainedGlobal(36);
1316   
1317   trackmult=esdTrackCutsGolden->GetReferenceMultiplicity(ESDevent, AliESDtrackCuts::MultEstTrackType(0), 0.5);
1318   
1319   //trackmult=AliESDtrackCuts::GetReferenceMultiplicity(ESDevent, AliESDtrackCuts::kTrackletsITSTPC, 0.5);
1320
1321
1322   Float_t spherocity=-1;
1323   Float_t sphericity=-1;
1324   Float_t spherocityTPC=-1;
1325   Float_t sphericityTPC=-1;
1326
1327
1328   spherocity=GetSpherocity(ESDevent,fTrackFilterGolden,0.8, 0.5, kFALSE);
1329   sphericity=GetSphericity(ESDevent,fTrackFilterGolden,0.8, 0.5, kFALSE);
1330   //Lot of memory consumption solved by reducing the number of files per job
1331   //spherocityTPC=GetSpherocity(ESDevent,fTrackFilterGolden,0.8, 0.5, kTRUE);
1332   //sphericityTPC=GetSphericity(ESDevent,fTrackFilterGolden,0.8, 0.5, kTRUE);
1333   
1334   spherocityTPC=spherocity;
1335   sphericityTPC=sphericity;
1336
1337
1338   if( analysisMode==kGlobalTrk ){  
1339     
1340     for(Int_t iT = 0; iT < nESDTracks; iT++) {
1341       
1342       AliESDtrack* esdTrack = ESDevent->GetTrack(iT);
1343       
1344       
1345       if(TMath::Abs(esdTrack->Eta()) > fEtaCut)
1346         continue;
1347       
1348       UShort_t filterFlag = 0;
1349       
1350       UInt_t selectDebug = 0;
1351       if (fTrackFilterGolden) {
1352         selectDebug = fTrackFilterGolden->IsSelected(esdTrack);
1353         if (selectDebug) {
1354           filterFlag +=1;
1355         }
1356       }
1357       
1358       if (fTrackFilterTPC) {
1359         
1360         selectDebug = fTrackFilterTPC->IsSelected(esdTrack);
1361         if (selectDebug){//only tracks which pass the TPC-only track cuts
1362           filterFlag +=2;
1363           
1364         }
1365         
1366       }
1367       
1368       if (fTrackFilter) {
1369         selectDebug = fTrackFilter->IsSelected(esdTrack);
1370         if (selectDebug) {
1371           filterFlag +=4;
1372         }
1373       }
1374       
1375      
1376       if(filterFlag==0)
1377         continue;
1378       
1379       Int_t sharedtpcclusters = esdTrack->GetTPCnclsS();
1380    
1381
1382       //
1383       // Track was accepted
1384       //      
1385       
1386       // Here we want to add histograms!
1387       
1388       if (esdTrack->Pt() < fMinPt) {
1389         
1390         // Keep small fraction of low pT tracks
1391         if(fRandom->Rndm() > fLowPtFraction)
1392           continue;
1393       } // else {
1394       // Here we want to add the high pt part of the histograms!
1395       //    }
1396     
1397       Short_t charge  = esdTrack->Charge();
1398       Float_t pt      = esdTrack->Pt();
1399       Float_t p       = esdTrack->P(); 
1400       Float_t eta     = esdTrack->Eta();
1401       Float_t phi     = esdTrack->Phi();
1402       Short_t ncl     = esdTrack->GetTPCsignalN();
1403       Short_t neff    = Short_t(esdTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1404       //          Short_t nclf    = esdTrack->GetTPCNclsF();
1405       Float_t dedx    = esdTrack->GetTPCsignal();
1406       Float_t tpcchi  = 0;
1407       if(esdTrack->GetTPCNcls() > 0)
1408         tpcchi = esdTrack->GetTPCchi2()/Float_t(esdTrack->GetTPCNcls());
1409       Float_t b[2];
1410       Float_t bCov[3];
1411       esdTrack->GetImpactParameters(b,bCov);
1412       Float_t dcaxy   = b[0];
1413       Float_t dcaz    = b[1];
1414       Double_t p_con[3] = {0, 0, 0};
1415       esdTrack->GetConstrainedPxPyPz(p_con);
1416       
1417       
1418       //        Float_t pt_con = (Float_t)TMath::Sqrt(p_con[0]*p_con[0] + p_con[1]*p_con[1]);
1419       // const AliExternalTrackParam* tpcParam = esdTrack->GetTPCInnerParam();
1420       // Float_t pttpc   = tpcParam->Pt();
1421       // Float_t ptpc    = tpcParam->P();
1422       
1423       Float_t ptMC        = 0;
1424       Short_t pidCode     = 0; // 0 = real data / no mc track!
1425       Short_t primaryFlag = 0; // 0 = real data / not primary mc track  
1426       Int_t   pdgMother   = 0;
1427       
1428       
1429       //with Globaltrack parameters????
1430       if(fAnalysisMC) {
1431         
1432         const Int_t label = TMath::Abs(esdTrack->GetLabel());
1433         TParticle* mcTrack = fMCStack->Particle(label);     
1434         if (mcTrack){
1435           
1436           if(fMCStack->IsPhysicalPrimary(label))
1437             primaryFlag = 1;
1438           
1439           //10/01/13. Add a flag to see if it is from material or WD
1440
1441           if(fMCStack->IsSecondaryFromWeakDecay(label))
1442             primaryFlag = 2;
1443
1444           if(fMCStack->IsSecondaryFromMaterial(label))
1445             primaryFlag = 3;
1446
1447
1448
1449
1450           Int_t pdgCode = mcTrack->GetPdgCode();
1451           pidCode = GetPidCode(pdgCode);
1452           
1453           ptMC      = mcTrack->Pt();
1454           
1455           TParticle* mother = FindPrimaryMother(fMCStack, label);
1456           pdgMother = mother->GetPdgCode();
1457         }
1458       }
1459       
1460     
1461       //TOF
1462       /*
1463         Float_t beta = -99;
1464         if (esdTrack->GetStatus()&AliESDtrack::kTOFpid){
1465         if ((esdTrack->GetIntegratedLength() != 0) && (esdTrack->GetTOFsignal() != 0))
1466           beta = esdTrack->GetIntegratedLength()/esdTrack->GetTOFsignal()/fgkClight;
1467           }
1468       */
1469       Bool_t IsTOFout=kFALSE;
1470       if ((esdTrack->GetStatus()&AliESDtrack::kTOFout)==0)
1471         IsTOFout=kTRUE;
1472       
1473       Bool_t HasTOFTime=kTRUE;
1474       if ((esdTrack->GetStatus()&AliESDtrack::kTIME)==0)
1475         HasTOFTime=kFALSE;
1476       
1477       Bool_t IsMatched=kFALSE;
1478       if (!(esdTrack->GetStatus()&AliESDtrack::kTOFmismatch)==0)
1479         IsMatched=kFALSE;
1480       else
1481         IsMatched=kTRUE;
1482
1483       Float_t lengthtrack=esdTrack->GetIntegratedLength();
1484
1485       Float_t timeTOF=esdTrack->GetTOFsignal();
1486
1487       Double_t inttime[5]={0,0,0,0,0}; 
1488       esdTrack->GetIntegratedTimes(inttime);// Returns the array with integrated times for each particle hypothesis
1489       Float_t fexptimeel=inttime[0];
1490       Float_t fexptimemu=inttime[1];
1491       Float_t fexptimepi=inttime[2];
1492       Float_t fexptimeka=inttime[3];
1493       Float_t fexptimepr=inttime[4];  
1494
1495       
1496       if(fTreeOption) {
1497         
1498         DeDxTrack* track = new((*fTrackArrayGlobalPar)[nadded]) DeDxTrack();
1499         nadded++;
1500         
1501         track->p          = p;
1502         track->pt         = pt;
1503         //        track->ptcon   = pt_con;
1504         //        track->tpcchi  = tpcchi;
1505         track->eta        = eta;
1506         track->phi        = phi;
1507         track->q          = charge;
1508         track->filter     = filterFlag;
1509         track->ncl        = ncl;
1510         track->neff       = neff;
1511         track->dedx       = dedx;
1512
1513         track->isTOFout   = IsTOFout;
1514         track->hasTOFtime = HasTOFTime;
1515         track->isTOFmatched = IsMatched;
1516         track->flength    = lengthtrack;
1517         track->ftimetof   = timeTOF;
1518         track->exptoftimeel = fexptimeel;
1519         track->exptoftimemu = fexptimemu;
1520         track->exptoftimepi = fexptimepi;
1521         track->exptoftimeka = fexptimeka;
1522         track->exptoftimepr = fexptimepr;
1523
1524         track->dcaxy      = dcaxy;
1525         track->dcaz       = dcaz;
1526         track->pid        = pidCode;
1527         track->primary    = primaryFlag;
1528         track->pttrue     = ptMC;
1529         track->mother     = pdgMother;
1530         track->tpcnclS    = sharedtpcclusters;
1531         
1532       }
1533     }//end of track loop
1534
1535   }//end global
1536   
1537   else if( analysisMode==kTPCTrk ){  
1538     const AliESDVertex *vtxSPD = ESDevent->GetPrimaryVertexSPD();
1539     if( vtxSPD->GetNContributors() < 1 || TMath::Abs(vtxSPD->GetZ()) > 10.0 ) return;
1540  
1541
1542     for(Int_t iT = 0; iT < nESDTracks; iT++) {
1543       
1544       AliESDtrack* esdTrack = ESDevent->GetTrack(iT);
1545
1546       AliESDtrack *trackTPC = AliESDtrackCuts::GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(ESDevent),esdTrack->GetID());
1547       
1548       if(!trackTPC) continue;
1549     
1550       UInt_t selectDebug = 0;
1551       selectDebug = fTrackFilterTPC->IsSelected(trackTPC);
1552       if(selectDebug==0) continue;
1553     
1554       if(trackTPC->Pt()>0.){
1555         // only constrain tracks above threshold
1556         AliExternalTrackParam exParam;
1557         // take the B-field from the ESD, no 3D fieldMap available at this point
1558         Bool_t relate = false;
1559         relate = trackTPC->RelateToVertexTPC(vtxSPD,ESDevent->GetMagneticField(),
1560                                            kVeryBig,&exParam);
1561         if(!relate){
1562           delete trackTPC;
1563           continue;
1564         }
1565         trackTPC->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),
1566                       exParam.GetCovariance());
1567       }
1568       else continue;     
1569       
1570  
1571       Int_t sharedtpcclusters = trackTPC->GetTPCnclsS();
1572
1573       if(TMath::Abs(trackTPC->Eta()) > fEtaCut)
1574         continue;
1575       
1576       //
1577       // Track was accepted
1578       //      
1579       
1580       // Here we want to add histograms!
1581       
1582       if (trackTPC->Pt() < fMinPt) {
1583         
1584         // Keep small fraction of low pT tracks
1585         if(fRandom->Rndm() > fLowPtFraction)
1586           continue;
1587       } // else {
1588       // Here we want to add the high pt part of the histograms!
1589       //    }
1590       
1591       Short_t charge  = trackTPC->Charge();
1592       Float_t pt      = trackTPC->Pt();
1593       Float_t p       = trackTPC->P(); 
1594       Float_t eta     = trackTPC->Eta();
1595       Float_t phi     = trackTPC->Phi();
1596       Short_t ncl     = trackTPC->GetTPCsignalN();
1597       Short_t neff    = Short_t(trackTPC->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1598       //          Short_t nclf    = esdTrack->GetTPCNclsF();
1599       Float_t dedx    = trackTPC->GetTPCsignal();
1600       Float_t tpcchi  = 0;
1601       if(trackTPC->GetTPCNcls() > 0)
1602         tpcchi = trackTPC->GetTPCchi2()/Float_t(trackTPC->GetTPCNcls());
1603       //Float_t b[2];
1604       //Float_t bCov[3];
1605       //trackTPC->GetImpactParameters(b,bCov);
1606       //Float_t dcaxy   = b[0];
1607       //Float_t dcaz    = b[1];
1608
1609       Float_t dcaxy = 0.;
1610       Float_t dcaz = 0.;
1611       trackTPC->GetImpactParameters(dcaxy,dcaz);
1612
1613
1614       Double_t p_con[3] = {0, 0, 0};
1615       trackTPC->GetConstrainedPxPyPz(p_con);
1616       
1617     
1618       // Float_t pt_con = (Float_t)TMath::Sqrt(p_con[0]*p_con[0] + p_con[1]*p_con[1]);
1619       // const AliExternalTrackParam* tpcParam = esdTrack->GetTPCInnerParam();
1620       // Float_t pttpc   = tpcParam->Pt();
1621       // Float_t ptpc    = tpcParam->P();
1622       
1623       Float_t ptMC        = 0;
1624       Short_t pidCode     = 0; // 0 = real data / no mc track!
1625       Short_t primaryFlag = 0; // 0 = real data / not primary mc track  
1626       Int_t   pdgMother   = 0;
1627       
1628       
1629       //with Globaltrack parameters????
1630       if(fAnalysisMC) {
1631         
1632         const Int_t label = TMath::Abs(esdTrack->GetLabel());
1633         TParticle* mcTrack = fMCStack->Particle(label);     
1634         if (mcTrack){
1635           
1636           if(fMCStack->IsPhysicalPrimary(label))
1637             primaryFlag = 1;
1638           
1639           
1640           //10/01/13. Add a flag to see if it is from material or WD
1641           
1642           if(fMCStack->IsSecondaryFromWeakDecay(label))
1643             primaryFlag = 2;
1644           
1645           if(fMCStack->IsSecondaryFromMaterial(label))
1646             primaryFlag = 3;
1647
1648
1649           Int_t pdgCode = mcTrack->GetPdgCode();
1650           pidCode = GetPidCode(pdgCode);
1651           
1652           ptMC      = mcTrack->Pt();
1653           
1654           TParticle* mother = FindPrimaryMother(fMCStack, label);
1655           pdgMother = mother->GetPdgCode();
1656         }
1657       }
1658     
1659     
1660       //TOF
1661       /*
1662       Float_t beta = -99;
1663       if (esdTrack->GetStatus()&AliESDtrack::kTOFpid){
1664         if ((esdTrack->GetIntegratedLength() != 0) && (esdTrack->GetTOFsignal() != 0))
1665           beta = esdTrack->GetIntegratedLength()/esdTrack->GetTOFsignal()/fgkClight;
1666       }
1667       */
1668       if(fTreeOption) {
1669         
1670         DeDxTrack* tracktpc = new((*fTrackArrayTPCPar)[nadded]) DeDxTrack();
1671         nadded++;
1672         
1673         tracktpc->p          = p;
1674         tracktpc->pt         = pt;
1675         //        track->ptcon   = pt_con;
1676         //        track->tpcchi  = tpcchi;
1677         tracktpc->eta        = eta;
1678         tracktpc->phi        = phi;
1679         tracktpc->q          = charge;
1680         tracktpc->filter     = 1;
1681         tracktpc->ncl        = ncl;
1682         tracktpc->neff       = neff;
1683         tracktpc->dedx       = dedx;
1684
1685         tracktpc->isTOFout   = 0;
1686         tracktpc->hasTOFtime = 0;
1687         tracktpc->isTOFmatched = 0;
1688         tracktpc->flength    = 0;
1689         tracktpc->ftimetof   = 0;
1690         tracktpc->exptoftimeel = 0;
1691         tracktpc->exptoftimemu = 0;
1692         tracktpc->exptoftimepi = 0;
1693         tracktpc->exptoftimeka = 0;
1694         tracktpc->exptoftimepr = 0;
1695
1696
1697         tracktpc->dcaxy      = dcaxy;
1698         tracktpc->dcaz       = dcaz;
1699         tracktpc->pid        = pidCode;
1700         tracktpc->primary    = primaryFlag;
1701         tracktpc->pttrue     = ptMC;
1702         tracktpc->mother     = pdgMother;
1703         tracktpc->tpcnclS    = sharedtpcclusters;
1704       }
1705
1706
1707       if(trackTPC)
1708         delete trackTPC;
1709
1710
1711     }//end of track loop
1712  
1713   }//end of: if isglobal==kFALSE
1714
1715
1716   if(fTreeOption) {
1717
1718     if( analysisMode==kGlobalTrk ){
1719       Sort(fTrackArrayGlobalPar, kFALSE);
1720       
1721       fEvent->trackmult = trackmult;
1722       fEvent->n         = nadded;
1723
1724       fEvent->sphericity         = sphericity;
1725       fEvent->spherocity         = spherocity;
1726       fEvent->sphericityTPC      = sphericityTPC;
1727       fEvent->spherocityTPC      = spherocityTPC;
1728
1729     }
1730     else if( analysisMode==kTPCTrk ){
1731       Sort(fTrackArrayTPCPar, kFALSE);
1732     }
1733
1734   }
1735
1736
1737 }
1738 //__________________________________________________________________
1739 void AliAnalysisTaskHighPtDeDx::ProduceArrayTrksAOD( AliAODEvent *AODevent, AnalysisMode analysisMode ){
1740
1741
1742   Int_t     trackmult = 0; // no pt cuts
1743   Int_t     nadded    = 0;
1744   Int_t nAODTracks = AODevent->GetNumberOfTracks();
1745   if( analysisMode == kGlobalTrk ){
1746     
1747     if(fTrackArrayGlobalPar)
1748       fTrackArrayGlobalPar->Clear();
1749     
1750   } 
1751   if( analysisMode == kTPCTrk ){
1752     if(fTrackArrayTPCPar)
1753       fTrackArrayTPCPar->Clear();
1754   }
1755
1756
1757
1758   //const AliAODVertex* vertexSPD= (AliAODVertex*)AODevent->GetPrimaryVertexSPD();//GetPrimaryVertex()
1759   //if( vertexSPD->GetNContributors() < 1 || TMath::Abs( vertexSPD->GetZ()) > 10.0 ) return;
1760
1761
1762   
1763   if( analysisMode==kGlobalTrk ){  
1764  
1765     
1766      const AliAODVertex* vertex = AODevent->GetPrimaryVertex();
1767     for(Int_t iT = 0; iT < nAODTracks; iT++) {
1768       
1769       AliAODTrack* aodTrack = dynamic_cast<AliAODTrack*>(AODevent->GetTrack(iT));
1770       if(!aodTrack) AliFatal("Not a standard AOD");
1771       
1772       //FOR AOD068, filterCut_Set2=KTRUE WHEN THE TRACK SATISFIES THE CUTS FROM JET ANALYSIS
1773       
1774       UShort_t filterFlag = 0;
1775
1776       
1777       if (fTrackFilterGolden) {
1778         
1779         // "Global track RAA analysis QM2011 + Chi2ITS<36"; bit 1024
1780         if(aodTrack->TestFilterBit(1024)) {
1781           filterFlag +=1;
1782         }
1783       }
1784       
1785       
1786       if (fTrackFilterTPC) {
1787         // TPC only cuts is bit 1 according to above macro
1788         // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX, 
1789         if(aodTrack->TestFilterBit(1)){
1790           filterFlag +=2;
1791         }
1792       }
1793       
1794     
1795       
1796       if (fTrackFilter) {
1797         /*
1798         // tighter cuts on primary particles for high pT tracks
1799         // take the standard cuts, which include already 
1800         // ITSrefit and use only primaries...
1801         bit 32
1802         */
1803         if(aodTrack->TestFilterBit(32)) {
1804           filterFlag +=4;
1805         }
1806
1807       }
1808       if(aodTrack->TestFilterBit(768)) {
1809         filterFlag +=8; // 272 for 2012 hljunggr, included to get hybrid tracks
1810         if(!aodTrack->IsHybridGlobalConstrainedGlobal()) filterFlag+=128; //just to test the selection, remove this line (/Martin)
1811       }
1812       
1813       
1814       // if(aodTrack->IsOn(AliAODTrack::kTPCrefit)) { //hljunggr: Tuvas flag for tpcrefit
1815       //        // Minimum number of clusters
1816       //        Float_t nCrossedRowsTPC = aodTrack->GetTPCClusterInfo(2,1);
1817       //   if (nCrossedRowsTPC >= 70) {
1818       //          Int_t findable=aodTrack->GetTPCNclsF();
1819       //          if (findable > 0){ 
1820       //            if (nCrossedRowsTPC/findable >= 0.8) filterFlag += 16;
1821       //          }
1822       //        }
1823         
1824
1825
1826       // }
1827       
1828
1829       
1830       if(filterFlag==0)
1831         continue;
1832       
1833       
1834       Double_t b[2], cov[3];
1835       if (!(aodTrack->PropagateToDCA(vertex, AODevent->GetMagneticField(), kVeryBig, b, cov)))
1836         filterFlag = 32; // propagation failed!!!!!;
1837       Float_t dcaxy   = b[0];
1838       Float_t dcaz    = b[1];
1839       
1840       TBits sharedTPC=aodTrack->GetTPCSharedMap();
1841       Int_t sharedtpcclusters=sharedTPC.CountBits(0)-sharedTPC.CountBits(159);
1842  
1843
1844
1845       // As I understand this routine recalculates the momentum so it should be called first!
1846       //Double_t b[2], cov[3];
1847   
1848       
1849       //if(!aodTrack->PropagateToDCA(vertex, AODevent->GetMagneticField(), kVeryBig, b, cov))
1850       //        filterFlag = 32; // propagation failed!!!!!
1851       
1852       if(TMath::Abs(aodTrack->Eta()) > fEtaCut)
1853         continue;
1854       if (aodTrack->Pt() < fMinPt) {
1855         
1856         // Keep small fraction of low pT tracks
1857         if(fRandom->Rndm() > fLowPtFraction)
1858           continue;
1859       } // else {
1860       // Here we want to add the high pt part of the histograms!
1861       //    }
1862      
1863      
1864
1865       
1866       Short_t charge  = aodTrack->Charge();
1867       Float_t pt      = aodTrack->Pt();
1868       Float_t p       = aodTrack->P(); 
1869       Float_t eta     = aodTrack->Eta();
1870       Float_t phi     = aodTrack->Phi();
1871       AliAODPid* aodPid = aodTrack->GetDetPid();
1872       Short_t ncl     = -10;
1873       Short_t neff    = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1874       //          Short_t nclf    = aodTrack->GetTPCNclsF();
1875       Float_t dedx    = -10;
1876       Bool_t IsTOFout=kFALSE;
1877       Bool_t HasTOFTime=kTRUE;
1878       Bool_t IsMatched=kFALSE;
1879       Float_t lengthtrack=-999;//in aod we do not have that information, beta must be: beta=inttime/timeTOF 
1880         
1881       Float_t timeTOF=-999;
1882       Float_t fexptimeel=-999;
1883       Float_t fexptimemu=-999;
1884       Float_t fexptimepi=-999;
1885       Float_t fexptimeka=-999;
1886       Float_t fexptimepr=-999;  
1887
1888
1889       //Float_t beta = -99;
1890       if(aodPid) {
1891         ncl     = aodPid->GetTPCsignalN();
1892         dedx    = aodPid->GetTPCsignal();
1893         //TOF
1894         /*
1895         if (aodTrack->GetStatus()&AliESDtrack::kTOFpid){
1896           Double_t tof[5];
1897           aodPid->GetIntegratedTimes(tof);
1898           beta = tof[0]/aodPid->GetTOFsignal();
1899           }*/
1900
1901
1902
1903         if ((aodTrack->GetStatus()&AliESDtrack::kTOFout)==0)
1904           IsTOFout=kTRUE;
1905         
1906
1907         if ((aodTrack->GetStatus()&AliESDtrack::kTIME)==0)
1908           HasTOFTime=kFALSE;
1909         
1910
1911         if (!(aodTrack->GetStatus()&AliESDtrack::kTOFmismatch)==0)
1912           IsMatched=kFALSE;
1913         else
1914           IsMatched=kTRUE;
1915         
1916         lengthtrack=-999;//in aod we do not have that information, beta must be: beta=inttime/timeTOF 
1917         
1918         timeTOF=aodPid->GetTOFsignal();
1919         
1920         Double_t inttime[5]={0,0,0,0,0}; 
1921         aodPid->GetIntegratedTimes(inttime);// Returns the array with integrated times for each particle hypothesis
1922         fexptimeel=inttime[0];
1923         fexptimemu=inttime[1];
1924         fexptimepi=inttime[2];
1925         fexptimeka=inttime[3];
1926         fexptimepr=inttime[4];  
1927         
1928
1929
1930       }
1931       //        Float_t tpcchi = aodTrack->Chi2perNDF();
1932       
1933       // Double_t p_con[3] = {0, 0, 0};
1934       // aodTrack->GetConstrainedPxPyPz(p_con);
1935       //        Float_t pt_con = 0; // This is not there! (Float_t)TMath::Sqrt(p_con[0]*p_con[0] + p_con[1]*p_con[1]);
1936       // const AliExternalTrackParam* tpcParam = aodTrack->GetTPCInnerParam();
1937       // Float_t pttpc   = tpcParam->Pt();
1938       // Float_t ptpc    = tpcParam->P();
1939       
1940       Float_t ptMC        = 0;
1941       Short_t pidCode     = 0; // 0 = real data / no mc track!
1942       Short_t primaryFlag = 0; // 0 = real data / not primary mc track  
1943       Int_t   pdgMother   = 0;
1944       
1945       if(fAnalysisMC) {
1946         
1947         const Int_t label = TMath::Abs(aodTrack->GetLabel());
1948         
1949         AliAODMCParticle* mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(label));
1950         if (mcTrack){
1951           
1952           if(mcTrack->IsPhysicalPrimary())
1953             primaryFlag = 1;
1954           
1955
1956           Int_t pdgCode = mcTrack->GetPdgCode();
1957           pidCode = GetPidCode(pdgCode);
1958           
1959           ptMC      = mcTrack->Pt();
1960           
1961           AliAODMCParticle* mother = FindPrimaryMotherAOD(mcTrack);
1962           pdgMother = mother->GetPdgCode();         
1963         }
1964       }
1965     
1966       
1967       if(fTreeOption) {
1968         
1969         DeDxTrack* track = new((*fTrackArrayGlobalPar)[nadded]) DeDxTrack();
1970         nadded++;
1971         
1972         track->p          = p;
1973         track->pt         = pt;
1974         //        track->ptcon   = pt_con;
1975         //        track->tpcchi  = tpcchi;
1976         track->eta        = eta;
1977         track->phi        = phi;
1978         track->q          = charge;
1979         track->filter     = filterFlag;
1980         track->ncl        = ncl;
1981         track->neff       = neff;
1982         track->dedx       = dedx;
1983         //track->beta       = beta;
1984         track->isTOFout   = IsTOFout;
1985         track->hasTOFtime = HasTOFTime;
1986         track->isTOFmatched = IsMatched;
1987         track->flength    = lengthtrack;
1988         track->ftimetof   = timeTOF;
1989         track->exptoftimeel = fexptimeel;
1990         track->exptoftimemu = fexptimemu;
1991         track->exptoftimepi = fexptimepi;
1992         track->exptoftimeka = fexptimeka;
1993         track->exptoftimepr = fexptimepr;
1994
1995         track->dcaxy      = dcaxy;
1996         track->dcaz       = dcaz;
1997         track->pid        = pidCode;
1998         track->primary    = primaryFlag;
1999         track->pttrue     = ptMC;
2000         track->mother     = pdgMother;
2001         track->tpcnclS    = sharedtpcclusters;
2002       }
2003     }//end of track loop
2004
2005  
2006
2007   }//end of global track analysis
2008   
2009
2010
2011   else if( analysisMode==kTPCTrk ){  
2012  
2013     const AliAODVertex* vertexSPD= (AliAODVertex*)AODevent->GetPrimaryVertexSPD();//GetPrimaryVertex()
2014     if( vertexSPD->GetNContributors() < 1 || TMath::Abs(vertexSPD->GetZ()) > 10.0 ) return;
2015
2016
2017     //const AliAODVertex* vertex = AODevent->GetPrimaryVertex();
2018     for(Int_t iT = 0; iT < nAODTracks; iT++) {
2019       
2020       AliAODTrack* aodTrack = dynamic_cast<AliAODTrack*>(AODevent->GetTrack(iT));
2021       if(!aodTrack) AliFatal("Not a standard AOD");
2022       
2023       UShort_t filterFlag = 0;
2024   
2025           
2026       // TPC only cuts is bit 1 according to above macro
2027       // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX, 
2028       if(aodTrack->TestFilterBit(128)){
2029         filterFlag +=2;
2030         trackmult++;
2031       }
2032       
2033       
2034       if(filterFlag==0)
2035         continue;
2036       
2037       TBits sharedTPC=aodTrack->GetTPCSharedMap();
2038       Int_t sharedtpcclusters=sharedTPC.CountBits(0)-sharedTPC.CountBits(159);
2039
2040       
2041       Double_t b[2], cov[3];
2042       //AliAODVertex* vertex = AODevent->GetPrimaryVertex();
2043       //Double_t b[2] = {-99., -99.};
2044       //Double_t bCov[3] = {-99., -99., -99.};
2045       //aodTrack->PropagateToDCA(aod->GetPrimaryVertex(), aod->GetMagneticField(), 100., b, bCov);
2046       if (!(aodTrack->PropagateToDCA(vertexSPD, AODevent->GetMagneticField(), kVeryBig, b, cov)))
2047         filterFlag = 32; // propagation failed!!!!!;
2048       Float_t dcaxy   = b[0];
2049       Float_t dcaz    = b[1];
2050
2051
2052
2053
2054       // As I understand this routine recalculates the momentum so it should be called first!
2055       //Double_t b[2], cov[3];
2056   
2057       
2058       //if(!aodTrack->PropagateToDCA(vertex, AODevent->GetMagneticField(), kVeryBig, b, cov))
2059       //        filterFlag = 32; // propagation failed!!!!!
2060       
2061       if(TMath::Abs(aodTrack->Eta()) > fEtaCut)
2062         continue;
2063       if (aodTrack->Pt() < fMinPt) {
2064         
2065         // Keep small fraction of low pT tracks
2066         if(fRandom->Rndm() > fLowPtFraction)
2067           continue;
2068       } // else {
2069       // Here we want to add the high pt part of the histograms!
2070       //    }
2071      
2072      
2073
2074       
2075       Short_t charge  = aodTrack->Charge();
2076       Float_t pt      = aodTrack->Pt();
2077       Float_t p       = aodTrack->P(); 
2078       Float_t eta     = aodTrack->Eta();
2079       Float_t phi     = aodTrack->Phi();
2080       AliAODPid* aodPid = aodTrack->GetDetPid();
2081       Short_t ncl     = -10;
2082       Short_t neff    = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
2083       //          Short_t nclf    = aodTrack->GetTPCNclsF();
2084       Float_t dedx    = -10;
2085       //Float_t beta = -99;
2086       if(aodPid) {
2087         ncl     = aodPid->GetTPCsignalN();
2088         dedx    = aodPid->GetTPCsignal();
2089         //TOF
2090         /*
2091         if (aodTrack->GetStatus()&AliESDtrack::kTOFpid){
2092           Double_t tof[5];
2093           aodPid->GetIntegratedTimes(tof);
2094           beta = tof[0]/aodPid->GetTOFsignal();
2095           }*/
2096       }
2097       //        Float_t tpcchi = aodTrack->Chi2perNDF();
2098       
2099       // Double_t p_con[3] = {0, 0, 0};
2100       // aodTrack->GetConstrainedPxPyPz(p_con);
2101       //        Float_t pt_con = 0; // This is not there! (Float_t)TMath::Sqrt(p_con[0]*p_con[0] + p_con[1]*p_con[1]);
2102       // const AliExternalTrackParam* tpcParam = aodTrack->GetTPCInnerParam();
2103       // Float_t pttpc   = tpcParam->Pt();
2104       // Float_t ptpc    = tpcParam->P();
2105       
2106       Float_t ptMC        = 0;
2107       Short_t pidCode     = 0; // 0 = real data / no mc track!
2108       Short_t primaryFlag = 0; // 0 = real data / not primary mc track  
2109       Int_t   pdgMother   = 0;
2110       
2111       if(fAnalysisMC) {
2112         
2113         const Int_t label = TMath::Abs(aodTrack->GetLabel());
2114         
2115         AliAODMCParticle* mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(label));
2116         if (mcTrack){
2117           
2118           if(mcTrack->IsPhysicalPrimary())
2119             primaryFlag = 1;
2120
2121           
2122           Int_t pdgCode = mcTrack->GetPdgCode();
2123           pidCode = GetPidCode(pdgCode);
2124           
2125           ptMC      = mcTrack->Pt();
2126           
2127           AliAODMCParticle* mother = FindPrimaryMotherAOD(mcTrack);
2128           pdgMother = mother->GetPdgCode();         
2129         }
2130       }
2131     
2132     
2133       if(fTreeOption) {
2134         
2135         DeDxTrack* tracktpc = new((*fTrackArrayTPCPar)[nadded]) DeDxTrack();
2136         nadded++;
2137         
2138         tracktpc->p          = p;
2139         tracktpc->pt         = pt;
2140         //        track->ptcon   = pt_con;
2141         //        track->tpcchi  = tpcchi;
2142         tracktpc->eta        = eta;
2143         tracktpc->phi        = phi;
2144         tracktpc->q          = charge;
2145         tracktpc->filter     = filterFlag;
2146         tracktpc->ncl        = ncl;
2147         tracktpc->neff       = neff;
2148         tracktpc->dedx       = dedx;
2149         //tracktpc->beta       = beta;
2150
2151         tracktpc->isTOFout   = 0;
2152         tracktpc->hasTOFtime = 0;
2153         tracktpc->isTOFmatched = 0;
2154         tracktpc->flength    = 0;
2155         tracktpc->ftimetof   = 0;
2156         tracktpc->exptoftimeel = 0;
2157         tracktpc->exptoftimemu = 0;
2158         tracktpc->exptoftimepi = 0;
2159         tracktpc->exptoftimeka = 0;
2160         tracktpc->exptoftimepr = 0;
2161
2162
2163         tracktpc->dcaxy      = dcaxy;
2164         tracktpc->dcaz       = dcaz;
2165         tracktpc->pid        = pidCode;
2166         tracktpc->primary    = primaryFlag;
2167         tracktpc->pttrue     = ptMC;
2168         tracktpc->mother     = pdgMother;
2169         tracktpc->tpcnclS       = sharedtpcclusters;
2170       }
2171     }//end of track loop
2172
2173
2174   }//end of global track analysis
2175
2176
2177   if(fTreeOption) {
2178     
2179     if( analysisMode==kGlobalTrk ){
2180       Sort(fTrackArrayGlobalPar, kFALSE);
2181       
2182       fEvent->trackmult = trackmult;
2183       fEvent->n         = nadded;
2184       fEvent->sphericity         = 0;
2185       fEvent->spherocity         = 0;
2186
2187     }
2188     if( analysisMode==kTPCTrk ){
2189       Sort(fTrackArrayTPCPar, kFALSE);
2190     }
2191     
2192   }
2193   
2194 }
2195 //__________________________________________________
2196 void AliAnalysisTaskHighPtDeDx::ProduceArrayV0ESD( AliESDEvent *ESDevent, AnalysisMode analysisMode ){
2197   Int_t nv0s = ESDevent->GetNumberOfV0s();
2198   if(nv0s<1)return;
2199   //Int_t     trackmult = 0; // no pt cuts
2200   Int_t     nadded    = 0;
2201   const AliESDVertex *myBestPrimaryVertex = ESDevent->GetPrimaryVertex();
2202   if (!myBestPrimaryVertex) return;
2203   if (!(myBestPrimaryVertex->GetStatus())) return;
2204   Double_t  lPrimaryVtxPosition[3];
2205   myBestPrimaryVertex->GetXYZ(lPrimaryVtxPosition);
2206   
2207   Double_t  lPrimaryVtxCov[6];
2208   myBestPrimaryVertex->GetCovMatrix(lPrimaryVtxCov);
2209   Double_t  lPrimaryVtxChi2 = myBestPrimaryVertex->GetChi2toNDF();
2210   
2211   AliAODVertex* myPrimaryVertex = new AliAODVertex(lPrimaryVtxPosition, lPrimaryVtxCov, lPrimaryVtxChi2, NULL, -1, AliAODVertex::kPrimary);
2212
2213
2214   if( analysisMode == kGlobalTrk ){
2215     if(fV0ArrayGlobalPar)
2216       fV0ArrayGlobalPar->Clear();
2217   } else if( analysisMode == kTPCTrk ){
2218     if(fV0ArrayTPCPar)
2219       fV0ArrayTPCPar->Clear();
2220   }
2221   
2222   if( analysisMode==kGlobalTrk ){  
2223
2224
2225
2226     
2227     // ################################
2228     // #### BEGINNING OF V0 CODE ######
2229     // ################################
2230
2231     
2232     for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
2233       
2234       // This is the begining of the V0 loop  
2235       AliESDv0 *esdV0 = ESDevent->GetV0(iV0);
2236       if (!esdV0) continue;
2237       
2238       // AliESDTrack (V0 Daughters)
2239       UInt_t lKeyPos = (UInt_t)TMath::Abs(esdV0->GetPindex());
2240       UInt_t lKeyNeg = (UInt_t)TMath::Abs(esdV0->GetNindex());
2241       
2242       AliESDtrack *pTrack = ESDevent->GetTrack(lKeyPos);
2243       AliESDtrack *nTrack = ESDevent->GetTrack(lKeyNeg);
2244       if (!pTrack || !nTrack) {
2245         Printf("ERROR: Could not retreive one of the daughter track");
2246         continue;
2247       }
2248       
2249       // Remove like-sign
2250       if (pTrack->GetSign() == nTrack->GetSign()) {
2251         //cout<< "like sign, continue"<< endl;
2252         continue;
2253       } 
2254       
2255       // Eta cut on decay products
2256       if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)
2257         continue;
2258       
2259         //Pre-selection to reduce output size
2260         // Pt cut on decay products
2261         if (esdV0->Pt() < fMinPtV0) continue;
2262         // No point in keeping low cospa values...
2263         if (esdV0->GetV0CosineOfPointingAngle() < 0.996 ) continue;
2264         //Reject on-the-fly tracks too
2265         if (esdV0->GetOnFlyStatus() != 0 ) continue;
2266         
2267
2268       //filter for positive track
2269       UShort_t filterFlag_p = 0;
2270       
2271       UInt_t selectDebug_p = 0;
2272       if (fTrackFilterGolden) {
2273         selectDebug_p = fTrackFilterGolden->IsSelected(pTrack);
2274         if (selectDebug_p) {
2275           filterFlag_p +=1;
2276         }
2277       }
2278       
2279       if (fTrackFilterTPC) {
2280         
2281         selectDebug_p = fTrackFilterTPC->IsSelected(pTrack);
2282         if (selectDebug_p){//only tracks which pass the TPC-only track cuts
2283           filterFlag_p +=2;
2284         }
2285         
2286       }
2287       
2288       if (fTrackFilter) {
2289         selectDebug_p = fTrackFilter->IsSelected(pTrack);
2290         if (selectDebug_p) {
2291           filterFlag_p +=4;
2292         }
2293       }
2294       
2295
2296     
2297       
2298
2299       //filter for negative track
2300       UShort_t filterFlag_n = 0;
2301       
2302       UInt_t selectDebug_n = 0;
2303       if (fTrackFilterGolden) {
2304         selectDebug_n = fTrackFilterGolden->IsSelected(nTrack);
2305         if (selectDebug_n) {
2306           filterFlag_n +=1;
2307         }
2308       }
2309       
2310       if (fTrackFilterTPC) {
2311         
2312         selectDebug_n = fTrackFilterTPC->IsSelected(nTrack);
2313         if (selectDebug_n){//only tracks which pass the TPC-only track cuts
2314           filterFlag_n +=2;
2315         }
2316         
2317       }
2318       
2319       if (fTrackFilter) {
2320         selectDebug_n = fTrackFilter->IsSelected(nTrack);
2321         if (selectDebug_n) {
2322           filterFlag_n +=4;
2323         }
2324       }
2325       
2326
2327
2328     
2329
2330
2331       // Check if switch does anything!
2332       Bool_t isSwitched = kFALSE;
2333       if (pTrack->GetSign() < 0) { // switch
2334         
2335         isSwitched = kTRUE;
2336         AliESDtrack* helpTrack = nTrack;
2337         nTrack = pTrack;
2338         pTrack = helpTrack;
2339       } 
2340       
2341       Float_t alpha = esdV0->AlphaV0();
2342       Float_t ptarm = esdV0->PtArmV0();
2343       // Double_t pVtxPos= v0->PrimaryVtxPosition();      
2344       
2345       Double_t  lV0Position[3];
2346       esdV0->GetXYZ(lV0Position[0], lV0Position[1], lV0Position[2]);
2347       
2348       Double_t lV0Radius      = TMath::Sqrt(lV0Position[0]*lV0Position[0]+lV0Position[1]*lV0Position[1]);
2349       Double_t lV0DecayLength = TMath::Sqrt(TMath::Power(lV0Position[0] - lPrimaryVtxPosition[0],2) +
2350                                             TMath::Power(lV0Position[1] - lPrimaryVtxPosition[1],2) +
2351                                             TMath::Power(lV0Position[2] - lPrimaryVtxPosition[2],2 ));
2352       AliKFVertex primaryVtxKF( *myPrimaryVertex );
2353       AliKFParticle::SetField(ESDevent->GetMagneticField());
2354       
2355       // Also implement switch here!!!!!!
2356       AliKFParticle* negEKF  = 0; // e-
2357       AliKFParticle* posEKF  = 0; // e+
2358       AliKFParticle* negPiKF = 0; // pi -
2359       AliKFParticle* posPiKF = 0; // pi +
2360       AliKFParticle* posPKF  = 0; // p
2361       AliKFParticle* negAPKF = 0; // p-bar
2362       
2363       if(!isSwitched) {
2364         negEKF  = new AliKFParticle( *(esdV0->GetParamN()) , 11);
2365         posEKF  = new AliKFParticle( *(esdV0->GetParamP()) ,-11);
2366         negPiKF = new AliKFParticle( *(esdV0->GetParamN()) ,-211);
2367         posPiKF = new AliKFParticle( *(esdV0->GetParamP()) , 211);
2368         posPKF  = new AliKFParticle( *(esdV0->GetParamP()) , 2212);
2369         negAPKF = new AliKFParticle( *(esdV0->GetParamN()) ,-2212);
2370       } else { // switch + and - 
2371         negEKF  = new AliKFParticle( *(esdV0->GetParamP()) , 11);
2372         posEKF  = new AliKFParticle( *(esdV0->GetParamN()) ,-11);
2373         negPiKF = new AliKFParticle( *(esdV0->GetParamP()) ,-211);
2374         posPiKF = new AliKFParticle( *(esdV0->GetParamN()) , 211);
2375         posPKF  = new AliKFParticle( *(esdV0->GetParamN()) , 2212);
2376         negAPKF = new AliKFParticle( *(esdV0->GetParamP()) ,-2212);
2377       }
2378       
2379       AliKFParticle v0GKF;  // Gamma e.g. from pi0
2380       v0GKF+=(*negEKF);
2381       v0GKF+=(*posEKF);
2382       v0GKF.SetProductionVertex(primaryVtxKF);
2383       
2384       AliKFParticle v0K0sKF; // K0 short
2385       v0K0sKF+=(*negPiKF);
2386       v0K0sKF+=(*posPiKF);
2387       v0K0sKF.SetProductionVertex(primaryVtxKF);
2388       
2389       AliKFParticle v0LambdaKF; // Lambda
2390       v0LambdaKF+=(*negPiKF);
2391       v0LambdaKF+=(*posPKF);    
2392       v0LambdaKF.SetProductionVertex(primaryVtxKF);
2393       
2394       AliKFParticle v0AntiLambdaKF; // Lambda-bar
2395       v0AntiLambdaKF+=(*posPiKF);
2396       v0AntiLambdaKF+=(*negAPKF);
2397       v0AntiLambdaKF.SetProductionVertex(primaryVtxKF);
2398       
2399       Double_t deltaInvMassG     = v0GKF.GetMass();
2400       Double_t deltaInvMassK0s   = v0K0sKF.GetMass()-0.498;
2401       Double_t deltaInvMassL     = v0LambdaKF.GetMass()-1.116;
2402       Double_t deltaInvMassAntiL = v0AntiLambdaKF.GetMass()-1.116;
2403
2404
2405       if(TMath::Abs(deltaInvMassK0s) > fMassCut &&
2406          TMath::Abs(deltaInvMassL) > fMassCut &&
2407          TMath::Abs(deltaInvMassAntiL) > fMassCut &&
2408         TMath::Abs(deltaInvMassG) > fMassCut )
2409         continue;
2410
2411       
2412       // Extract track information
2413       
2414       Short_t pcharge  = pTrack->Charge();
2415       Float_t ppt      = pTrack->Pt();
2416       Float_t pp       = pTrack->P(); 
2417       Float_t peta     = pTrack->Eta();
2418       Float_t pphi     = pTrack->Phi();
2419       Short_t pncl     = pTrack->GetTPCsignalN();
2420       Short_t pneff    = Short_t(pTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
2421       Float_t pdedx    = pTrack->GetTPCsignal();
2422       Int_t   psharedtpcclusters = pTrack->GetTPCnclsS();
2423       
2424       Float_t ptpcchi  = 0;
2425       if(pTrack->GetTPCNcls() > 0)
2426         ptpcchi = pTrack->GetTPCchi2()/Float_t(pTrack->GetTPCNcls());
2427       
2428
2429       Bool_t pIsTOFout=kFALSE;
2430       if ((pTrack->GetStatus()&AliESDtrack::kTOFout)==0)
2431         pIsTOFout=kTRUE;
2432       
2433       Bool_t pHasTOFTime=kTRUE;
2434       if ((pTrack->GetStatus()&AliESDtrack::kTIME)==0)
2435         pHasTOFTime=kFALSE;
2436       
2437       Bool_t pIsMatched=kFALSE;
2438       if (!(pTrack->GetStatus()&AliESDtrack::kTOFmismatch)==0)
2439         pIsMatched=kFALSE;
2440       else
2441         pIsMatched=kTRUE;
2442
2443       Float_t plengthtrack=pTrack->GetIntegratedLength();
2444
2445       Float_t ptimeTOF=pTrack->GetTOFsignal();
2446
2447       Double_t pinttime[5]={0,0,0,0,0}; 
2448       pTrack->GetIntegratedTimes(pinttime);// Returns the array with integrated times for each particle hypothesis
2449       Float_t pfexptimeel=pinttime[0];
2450       Float_t pfexptimemu=pinttime[1];
2451       Float_t pfexptimepi=pinttime[2];
2452       Float_t pfexptimeka=pinttime[3];
2453       Float_t pfexptimepr=pinttime[4];  
2454
2455
2456
2457
2458       Short_t ncharge  = nTrack->Charge();
2459       Float_t npt      = nTrack->Pt();
2460       Float_t np       = nTrack->P(); 
2461       Float_t neta     = nTrack->Eta();
2462       Float_t nphi     = nTrack->Phi();
2463       Short_t nncl     = nTrack->GetTPCsignalN();
2464       Short_t nneff    = Short_t(nTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
2465       Float_t ndedx    = nTrack->GetTPCsignal();
2466       Int_t   nsharedtpcclusters = nTrack->GetTPCnclsS();
2467
2468       Float_t ntpcchi  = 0;
2469       if(nTrack->GetTPCNcls() > 0)
2470         ntpcchi = nTrack->GetTPCchi2()/Float_t(nTrack->GetTPCNcls());
2471       
2472
2473
2474       Bool_t nIsTOFout=kFALSE;
2475       if ((nTrack->GetStatus()&AliESDtrack::kTOFout)==0)
2476         nIsTOFout=kTRUE;
2477       
2478       Bool_t nHasTOFTime=kTRUE;
2479       if ((nTrack->GetStatus()&AliESDtrack::kTIME)==0)
2480         nHasTOFTime=kFALSE;
2481       
2482       Bool_t nIsMatched=kFALSE;
2483       if (!(nTrack->GetStatus()&AliESDtrack::kTOFmismatch)==0)
2484         nIsMatched=kFALSE;
2485       else
2486         nIsMatched=kTRUE;
2487
2488       Float_t nlengthtrack=nTrack->GetIntegratedLength();
2489
2490       Float_t ntimeTOF=nTrack->GetTOFsignal();
2491
2492       Double_t ninttime[5]={0,0,0,0,0}; 
2493       nTrack->GetIntegratedTimes(ninttime);// Returns the array with integrated times for each particle hypothesis
2494       Float_t nfexptimeel=ninttime[0];
2495       Float_t nfexptimemu=ninttime[1];
2496       Float_t nfexptimepi=ninttime[2];
2497       Float_t nfexptimeka=ninttime[3];
2498       Float_t nfexptimepr=ninttime[4];  
2499
2500
2501
2502       Float_t b[2];
2503       Float_t bCov[3];
2504       pTrack->GetImpactParameters(b,bCov);
2505       Float_t pdcaxy   = b[0];
2506       Float_t pdcaz    = b[1];
2507       nTrack->GetImpactParameters(b,bCov);
2508       Float_t ndcaxy   = b[0];
2509       Float_t ndcaz    = b[1];
2510       
2511       Int_t   primaryV0     = 0; // 0 means that the tracks are not both daughters of a primary particle (1 means they are)
2512       Int_t   pdgV0         = 0; // 0 means that they don't have same origin for MC (1 means they have the same original mother)
2513       Float_t p_ptMC        = 0;
2514       Short_t p_pidCode     = 0; // 0 = real data / no mc track!
2515       Short_t p_primaryFlag = 0; // 0 = real data / not primary mc track  
2516       Int_t   p_pdgMother   = 0;
2517       Float_t n_ptMC        = 0;
2518       Short_t n_pidCode     = 0; // 0 = real data / no mc track!
2519       Short_t n_primaryFlag = 0; // 0 = real data / not primary mc track  
2520       Int_t   n_pdgMother   = 0;
2521       if(fAnalysisMC) {
2522         
2523         Int_t p_mother_label = 0;
2524         Int_t p_mother_steps = 0;
2525         Int_t n_mother_label = 0;
2526         Int_t n_mother_steps = 0;
2527         
2528         // positive track
2529         const Int_t p_label = TMath::Abs(pTrack->GetLabel());
2530         TParticle* p_mcTrack = fMCStack->Particle(p_label);         
2531         if (p_mcTrack){
2532           
2533           if(fMCStack->IsPhysicalPrimary(p_label))
2534             p_primaryFlag = 1;
2535           
2536           
2537           //10/01/13. Add a flag to see if it is from material or WD
2538
2539           if(fMCStack->IsSecondaryFromWeakDecay(p_label))
2540             p_primaryFlag = 2;
2541
2542           if(fMCStack->IsSecondaryFromMaterial(p_label))
2543             p_primaryFlag = 3;
2544
2545           Int_t p_pdgCode = p_mcTrack->GetPdgCode();
2546           p_pidCode = GetPidCode(p_pdgCode);
2547           
2548           p_ptMC      = p_mcTrack->Pt();
2549           
2550           p_mother_label = FindPrimaryMotherLabelV0(fMCStack, p_label, 
2551                                                   p_mother_steps);
2552           if(p_mother_label>0) {
2553             TParticle* p_mother = fMCStack->Particle(p_mother_label);
2554             p_pdgMother = p_mother->GetPdgCode();
2555           }
2556         }
2557         
2558         // negative track
2559         const Int_t n_label = TMath::Abs(nTrack->GetLabel());
2560         TParticle* n_mcTrack = fMCStack->Particle(n_label);         
2561         if (n_mcTrack){
2562           
2563           if(fMCStack->IsPhysicalPrimary(n_label))
2564             n_primaryFlag = 1;
2565           
2566           //10/01/13. Add a flag to see if it is from material or WD
2567
2568           if(fMCStack->IsSecondaryFromWeakDecay(n_label))
2569             n_primaryFlag = 2;
2570
2571           if(fMCStack->IsSecondaryFromMaterial(n_label))
2572             n_primaryFlag = 3;
2573
2574           Int_t n_pdgCode = n_mcTrack->GetPdgCode();
2575           n_pidCode = GetPidCode(n_pdgCode);
2576           
2577           n_ptMC      = n_mcTrack->Pt();
2578           
2579           n_mother_label = FindPrimaryMotherLabelV0(fMCStack, n_label, 
2580                                                   n_mother_steps);
2581           if(n_mother_label>0) {
2582             TParticle* n_mother = fMCStack->Particle(n_mother_label);
2583             n_pdgMother = n_mother->GetPdgCode();
2584           }
2585         }
2586         
2587         // Check if V0 is primary = first and the same mother of both partciles
2588         if(p_mother_label>0 && n_mother_label>0 && p_mother_label == n_mother_label) {
2589           pdgV0 = p_pdgMother;
2590           if(p_mother_steps == 1 && n_mother_steps == 1) {
2591             primaryV0 = 1;
2592           }
2593         }
2594       }
2595
2596  
2597     
2598       if(fTreeOption) {
2599         
2600
2601
2602         DeDxV0* v0data = new((*fV0ArrayGlobalPar)[nadded]) DeDxV0();
2603         nadded++;
2604         
2605         // v0 data
2606         v0data->p       = esdV0->P();
2607         v0data->pt      = esdV0->Pt();
2608         v0data->eta     = esdV0->Eta();
2609         v0data->phi     = esdV0->Phi();
2610         v0data->pdca    = TMath::Sqrt(pdcaxy*pdcaxy + pdcaz*pdcaz);
2611         v0data->ndca    = TMath::Sqrt(ndcaxy*ndcaxy + ndcaz*ndcaz);
2612         v0data->dmassG  = deltaInvMassG;
2613         v0data->dmassK0 = deltaInvMassK0s;
2614         v0data->dmassL  = deltaInvMassL;
2615         v0data->dmassAL = deltaInvMassAntiL;
2616         v0data->alpha   = alpha;
2617         v0data->ptarm   = ptarm;
2618         v0data->decayr  = lV0Radius;
2619         v0data->decayl  = lV0DecayLength;
2620         
2621         // New parameters
2622         v0data->status  = esdV0->GetOnFlyStatus();
2623         v0data->chi2    = esdV0->GetChi2V0();
2624         v0data->cospt   = esdV0->GetV0CosineOfPointingAngle(); 
2625         // cospt: as I understand this means that the pointing to the vertex
2626         // is fine so I remove the dcaxy and dcaz for the V= class
2627         v0data->dcadaughters = esdV0->GetDcaV0Daughters();
2628         v0data->primary = primaryV0;
2629         v0data->pdg     = pdgV0;
2630         
2631         // positive track
2632         v0data->ptrack.p       = pp;
2633         v0data->ptrack.pt      = ppt;
2634         //        v0data->ptrack.ptcon   = ppt_con;
2635         //        v0data->ptrack.tpcchi  = ptpcchi;
2636         v0data->ptrack.eta     = peta;
2637         v0data->ptrack.phi     = pphi;
2638         v0data->ptrack.q       = pcharge;
2639         v0data->ptrack.ncl     = pncl;
2640         v0data->ptrack.neff    = pneff;
2641         v0data->ptrack.dedx    = pdedx;
2642
2643
2644         v0data->ptrack.isTOFout   = pIsTOFout;
2645         v0data->ptrack.hasTOFtime = pHasTOFTime;
2646         v0data->ptrack.isTOFmatched = pIsMatched;
2647         v0data->ptrack.flength    =   plengthtrack;
2648         v0data->ptrack.ftimetof   = ptimeTOF;
2649         v0data->ptrack.exptoftimeel = pfexptimeel;
2650         v0data->ptrack.exptoftimemu = pfexptimemu;
2651         v0data->ptrack.exptoftimepi = pfexptimepi;
2652         v0data->ptrack.exptoftimeka = pfexptimeka;
2653         v0data->ptrack.exptoftimepr = pfexptimepr;
2654
2655         v0data->ptrack.dcaxy   = pdcaxy;
2656         v0data->ptrack.dcaz    = pdcaz;
2657         v0data->ptrack.pid     = p_pidCode;
2658         v0data->ptrack.primary = p_primaryFlag;
2659         v0data->ptrack.pttrue  = p_ptMC;
2660         v0data->ptrack.mother  = p_pdgMother;
2661         v0data->ptrack.filter  = filterFlag_p;
2662         v0data->ptrack.tpcnclS    = psharedtpcclusters;
2663
2664         // negative track
2665         v0data->ntrack.p       = np;
2666         v0data->ntrack.pt      = npt;
2667         //        v0data->ntrack.ptcon   = npt_con;
2668         //        v0data->ntrack.tpcchi  = ntpcchi;
2669         v0data->ntrack.eta     = neta;
2670         v0data->ntrack.phi     = nphi;
2671         v0data->ntrack.q       = ncharge;
2672         v0data->ntrack.ncl     = nncl;
2673         v0data->ntrack.neff    = nneff;
2674         v0data->ntrack.dedx    = ndedx;
2675
2676         v0data->ntrack.isTOFout   = nIsTOFout;
2677         v0data->ntrack.hasTOFtime = nHasTOFTime;
2678         v0data->ntrack.isTOFmatched = nIsMatched;
2679         v0data->ntrack.flength    =   nlengthtrack;
2680         v0data->ntrack.ftimetof   = ntimeTOF;
2681         v0data->ntrack.exptoftimeel = nfexptimeel;
2682         v0data->ntrack.exptoftimemu = nfexptimemu;
2683         v0data->ntrack.exptoftimepi = nfexptimepi;
2684         v0data->ntrack.exptoftimeka = nfexptimeka;
2685         v0data->ntrack.exptoftimepr = nfexptimepr;
2686
2687
2688         v0data->ntrack.dcaxy   = ndcaxy;
2689         v0data->ntrack.dcaz    = ndcaz;
2690         v0data->ntrack.pid     = n_pidCode;
2691         v0data->ntrack.primary = n_primaryFlag;
2692         v0data->ntrack.pttrue  = n_ptMC;
2693         v0data->ntrack.mother  = n_pdgMother;
2694         v0data->ntrack.filter  = filterFlag_n;
2695         v0data->ntrack.tpcnclS    = nsharedtpcclusters;
2696
2697       }
2698       
2699       // clean up loop over v0
2700
2701       delete negPiKF;
2702       delete posPiKF;
2703       delete posPKF;
2704       delete negAPKF;
2705   
2706
2707
2708     }
2709   
2710     // clean up event
2711     //delete myPrimaryVertex;
2712   }
2713   else if( analysisMode==kTPCTrk ){  
2714  
2715     const AliESDVertex *vtxSPD = ESDevent->GetPrimaryVertexSPD();
2716     if( vtxSPD->GetNContributors() < 1 || TMath::Abs(vtxSPD->GetZ()) > 10.0 ) return;
2717     
2718     
2719     // ################################
2720     // #### BEGINNING OF V0 CODE ######
2721     // ################################
2722
2723     
2724     for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
2725       
2726       // This is the begining of the V0 loop  
2727       AliESDv0 *esdV0 = ESDevent->GetV0(iV0);
2728       if (!esdV0) continue;
2729       
2730       // AliESDTrack (V0 Daughters)
2731       UInt_t lKeyPos = (UInt_t)TMath::Abs(esdV0->GetPindex());
2732       UInt_t lKeyNeg = (UInt_t)TMath::Abs(esdV0->GetNindex());
2733       
2734       AliESDtrack *pTrack = ESDevent->GetTrack(lKeyPos);
2735       AliESDtrack *nTrack = ESDevent->GetTrack(lKeyNeg);
2736       if (!pTrack || !nTrack) {
2737         Printf("ERROR: Could not retreive one of the daughter track");
2738         continue;
2739       }
2740
2741       // Remove like-sign
2742       if (pTrack->GetSign() == nTrack->GetSign()) {
2743         //cout<< "like sign, continue"<< endl;
2744         continue;
2745       } 
2746
2747       AliESDtrack *pTrackTPC = AliESDtrackCuts::GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(ESDevent),pTrack->GetID());
2748       AliESDtrack *nTrackTPC = AliESDtrackCuts::GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(ESDevent),nTrack->GetID());
2749
2750       if (!pTrackTPC || !nTrackTPC) {
2751         Printf("ERROR: Could not retreive one of the daughter TPC track");
2752         continue;
2753       }
2754
2755       //filter for positive track
2756       UInt_t selectDebug_p = 0;
2757       UShort_t filterFlag_p = 0;
2758       selectDebug_p = fTrackFilterTPC->IsSelected(pTrackTPC);
2759       //if(selectDebug_p==0) continue;
2760
2761       //filter for negative track
2762       UInt_t selectDebug_n = 0;
2763       UShort_t filterFlag_n = 0;
2764       selectDebug_n = fTrackFilterTPC->IsSelected(nTrackTPC);
2765       //if(selectDebug_n==0 || selectDebug_p==0) continue;
2766
2767       if(selectDebug_p)
2768         filterFlag_p += 1;
2769
2770       if(pTrackTPC->Pt()>0.){
2771         // only constrain tracks above threshold
2772         AliExternalTrackParam exParamp;
2773         // take the B-field from the ESD, no 3D fieldMap available at this point
2774         Bool_t relate_p = false;
2775         relate_p = pTrackTPC->RelateToVertexTPC(vtxSPD,ESDevent->GetMagneticField(),
2776                                                 kVeryBig,&exParamp);
2777         if(!relate_p){
2778           delete pTrackTPC;
2779           continue;
2780         }
2781         pTrackTPC->Set(exParamp.GetX(),exParamp.GetAlpha(),exParamp.GetParameter(),
2782                        exParamp.GetCovariance());
2783       }
2784       else continue;     
2785       
2786
2787       //filter for negative track
2788  
2789       if(selectDebug_n)
2790         filterFlag_n += 1;
2791
2792   
2793       if(nTrackTPC->Pt()>0.){
2794         // only constrain tracks above threshold
2795         AliExternalTrackParam exParamn;
2796         // take the B-field from the ESD, no 3D fieldMap available at this point
2797         Bool_t relate_n = false;
2798         relate_n = nTrackTPC->RelateToVertexTPC(vtxSPD,ESDevent->GetMagneticField(),
2799                                                 kVeryBig,&exParamn);
2800         if(!relate_n){
2801           delete nTrackTPC;
2802           continue;
2803         }
2804         nTrackTPC->Set(exParamn.GetX(),exParamn.GetAlpha(),exParamn.GetParameter(),
2805                        exParamn.GetCovariance());
2806       }
2807       else continue;  
2808       
2809       
2810       // Eta cut on decay products
2811       if(TMath::Abs(pTrackTPC->Eta()) > fEtaCut || TMath::Abs(nTrackTPC->Eta()) > fEtaCut)
2812         continue;
2813       
2814         //Pre-selection to reduce output size
2815         // Pt cut on decay products
2816         if (esdV0->Pt() < fMinPtV0) continue;
2817         // No point in keeping low cospa values...
2818         if (esdV0->GetV0CosineOfPointingAngle() < 0.996 ) continue;
2819         //Reject on-the-fly tracks too
2820         if (esdV0->GetOnFlyStatus() != 0 ) continue;
2821  
2822       
2823       // Check if switch does anything!
2824       Bool_t isSwitched = kFALSE;
2825       if (pTrackTPC->GetSign() < 0) { // switch
2826         
2827         isSwitched = kTRUE;
2828         AliESDtrack* helpTrack = nTrack;
2829         nTrackTPC = pTrackTPC;
2830         pTrackTPC = helpTrack;
2831       } 
2832       
2833
2834       // Extract track information
2835       
2836       Short_t pcharge  = pTrackTPC->Charge();
2837       Float_t ppt      = pTrackTPC->Pt();
2838       Float_t pp       = pTrackTPC->P(); 
2839       Float_t peta     = pTrackTPC->Eta();
2840       Float_t pphi     = pTrackTPC->Phi();
2841       Short_t pncl     = pTrackTPC->GetTPCsignalN();
2842       Short_t pneff    = Short_t(pTrackTPC->GetTPCClusterInfo(2, 1)); // effective track length for pT res
2843       Float_t pdedx    = pTrackTPC->GetTPCsignal();
2844       Int_t   psharedtpcclusters = pTrackTPC->GetTPCnclsS();
2845
2846       Float_t ptpcchi  = 0;
2847       if(pTrackTPC->GetTPCNcls() > 0)
2848         ptpcchi = pTrackTPC->GetTPCchi2()/Float_t(pTrackTPC->GetTPCNcls());
2849       
2850       Short_t ncharge  = nTrackTPC->Charge();
2851       Float_t npt      = nTrackTPC->Pt();
2852       Float_t np       = nTrackTPC->P(); 
2853       Float_t neta     = nTrackTPC->Eta();
2854       Float_t nphi     = nTrackTPC->Phi();
2855       Short_t nncl     = nTrackTPC->GetTPCsignalN();
2856       Short_t nneff    = Short_t(nTrackTPC->GetTPCClusterInfo(2, 1)); // effective track length for pT res
2857       Float_t ndedx    = nTrackTPC->GetTPCsignal();
2858       Int_t   nsharedtpcclusters = nTrackTPC->GetTPCnclsS();
2859
2860       Float_t ntpcchi  = 0;
2861       if(nTrackTPC->GetTPCNcls() > 0)
2862         ntpcchi = nTrackTPC->GetTPCchi2()/Float_t(nTrackTPC->GetTPCNcls());
2863       
2864       Float_t bp[2]={0,0};
2865       Float_t bCovp[3]={0,0,0};
2866       pTrackTPC->GetImpactParameters(bp,bCovp);
2867       Float_t pdcaxy   = bp[0];
2868       Float_t pdcaz    = bp[1];
2869
2870       Float_t bn[2]={0,0};
2871       Float_t bCovn[3]={0,0,0};
2872       nTrackTPC->GetImpactParameters(bn,bCovn);
2873       Float_t ndcaxy   = bn[0];
2874       Float_t ndcaz    = bn[1];
2875
2876
2877       Float_t alpha = esdV0->AlphaV0();
2878       Float_t ptarm = esdV0->PtArmV0();
2879       // Double_t pVtxPos= v0->PrimaryVtxPosition();      
2880       
2881       Double_t  lV0Position[3];
2882       esdV0->GetXYZ(lV0Position[0], lV0Position[1], lV0Position[2]);
2883       
2884       Double_t lV0Radius      = TMath::Sqrt(lV0Position[0]*lV0Position[0]+lV0Position[1]*lV0Position[1]);
2885       Double_t lV0DecayLength = TMath::Sqrt(TMath::Power(lV0Position[0] - lPrimaryVtxPosition[0],2) +
2886                                             TMath::Power(lV0Position[1] - lPrimaryVtxPosition[1],2) +
2887                                             TMath::Power(lV0Position[2] - lPrimaryVtxPosition[2],2 ));
2888       AliKFVertex primaryVtxKF( *myPrimaryVertex );
2889       AliKFParticle::SetField(ESDevent->GetMagneticField());
2890       
2891       // Also implement switch here!!!!!!
2892       AliKFParticle* negEKF  = 0; // e-
2893       AliKFParticle* posEKF  = 0; // e+
2894       AliKFParticle* negPiKF = 0; // pi -
2895       AliKFParticle* posPiKF = 0; // pi +
2896       AliKFParticle* posPKF  = 0; // p
2897       AliKFParticle* negAPKF = 0; // p-bar
2898       
2899       
2900       if(!isSwitched) {
2901         /*      
2902                 negEKF  = new AliKFParticle( *(esdV0->GetParamN()) , 11);
2903                 posEKF  = new AliKFParticle( *(esdV0->GetParamP()) ,-11);
2904                 negPiKF = new AliKFParticle( *(esdV0->GetParamN()) ,-211);
2905                 posPiKF = new AliKFParticle( *(esdV0->GetParamP()) , 211);
2906                 posPKF  = new AliKFParticle( *(esdV0->GetParamP()) , 2212);
2907                 negAPKF = new AliKFParticle( *(esdV0->GetParamN()) ,-2212);
2908         */
2909
2910         negEKF  = new AliKFParticle( *(dynamic_cast<AliVTrack*>(nTrackTPC)) , 11);
2911         posEKF  = new AliKFParticle( *(dynamic_cast<AliVTrack*>(pTrackTPC)) ,-11);
2912         negPiKF = new AliKFParticle( *(dynamic_cast<AliVTrack*>(nTrackTPC)) ,-211);
2913         posPiKF = new AliKFParticle( *(dynamic_cast<AliVTrack*>(pTrackTPC)) , 211);
2914         posPKF  = new AliKFParticle( *(dynamic_cast<AliVTrack*>(pTrackTPC)) , 2212);
2915         negAPKF = new AliKFParticle( *(dynamic_cast<AliVTrack*>(nTrackTPC)) ,-2212);
2916
2917
2918
2919       } else { // switch + and - 
2920         /*
2921         negEKF  = new AliKFParticle( *(esdV0->GetParamP()) , 11);
2922         posEKF  = new AliKFParticle( *(esdV0->GetParamN()) ,-11);
2923         negPiKF = new AliKFParticle( *(esdV0->GetParamP()) ,-211);
2924         posPiKF = new AliKFParticle( *(esdV0->GetParamN()) , 211);
2925         posPKF  = new AliKFParticle( *(esdV0->GetParamN()) , 2212);
2926         negAPKF = new AliKFParticle( *(esdV0->GetParamP()) ,-2212);
2927         */
2928
2929         negEKF  = new AliKFParticle(  *(dynamic_cast<AliVTrack*>(pTrackTPC)), 11);
2930         posEKF  = new AliKFParticle(  *(dynamic_cast<AliVTrack*>(nTrackTPC)),-11);
2931         negPiKF = new AliKFParticle(  *(dynamic_cast<AliVTrack*>(pTrackTPC)),-211);
2932         posPiKF = new AliKFParticle(  *(dynamic_cast<AliVTrack*>(nTrackTPC)), 211);
2933         posPKF  = new AliKFParticle(  *(dynamic_cast<AliVTrack*>(nTrackTPC)), 2212);
2934         negAPKF = new AliKFParticle(  *(dynamic_cast<AliVTrack*>(pTrackTPC)),-2212);
2935
2936
2937       }
2938    
2939
2940
2941
2942       
2943       AliKFParticle v0GKF;  // Gamma e.g. from pi0
2944       v0GKF+=(*negEKF);
2945       v0GKF+=(*posEKF);
2946       v0GKF.SetProductionVertex(primaryVtxKF);
2947       
2948       AliKFParticle v0K0sKF; // K0 short
2949       v0K0sKF+=(*negPiKF);
2950       v0K0sKF+=(*posPiKF);
2951       v0K0sKF.SetProductionVertex(primaryVtxKF);
2952       
2953       AliKFParticle v0LambdaKF; // Lambda
2954       v0LambdaKF+=(*negPiKF);
2955       v0LambdaKF+=(*posPKF);    
2956       v0LambdaKF.SetProductionVertex(primaryVtxKF);
2957       
2958       AliKFParticle v0AntiLambdaKF; // Lambda-bar
2959       v0AntiLambdaKF+=(*posPiKF);
2960       v0AntiLambdaKF+=(*negAPKF);
2961       v0AntiLambdaKF.SetProductionVertex(primaryVtxKF);
2962       
2963       Double_t deltaInvMassG     = v0GKF.GetMass();
2964       Double_t deltaInvMassK0s   = v0K0sKF.GetMass()-0.498;
2965       Double_t deltaInvMassL     = v0LambdaKF.GetMass()-1.116;
2966       Double_t deltaInvMassAntiL = v0AntiLambdaKF.GetMass()-1.116;
2967       
2968       
2969       if(TMath::Abs(deltaInvMassK0s) > fMassCut &&
2970          TMath::Abs(deltaInvMassL) > fMassCut &&
2971          TMath::Abs(deltaInvMassAntiL) > fMassCut &&
2972          TMath::Abs(deltaInvMassG) > fMassCut )
2973         continue;
2974       
2975
2976
2977       // TODO: Whe should these be different? Different mass hypothesis = energy loss
2978       // This is not important for us as we focus on the decay products!
2979       // Double_t ptK0s        = v0K0sKF.GetPt(); 
2980       // Double_t ptL          = v0LambdaKF.GetPt();
2981       // Double_t ptAntiL      = v0AntiLambdaKF.GetPt();     
2982       
2983       Int_t   primaryV0     = 0; // 0 means that the tracks are not both daughters of a primary particle (1 means they are)
2984       Int_t   pdgV0         = 0; // 0 means that they don't have same origin for MC (1 means they have the same original mother)
2985       Float_t p_ptMC        = 0;
2986       Short_t p_pidCode     = 0; // 0 = real data / no mc track!
2987       Short_t p_primaryFlag = 0; // 0 = real data / not primary mc track  
2988       Int_t   p_pdgMother   = 0;
2989       Float_t n_ptMC        = 0;
2990       Short_t n_pidCode     = 0; // 0 = real data / no mc track!
2991       Short_t n_primaryFlag = 0; // 0 = real data / not primary mc track  
2992       Int_t   n_pdgMother   = 0;
2993       if(fAnalysisMC) {
2994         
2995         Int_t p_mother_label = 0;
2996         Int_t p_mother_steps = 0;
2997         Int_t n_mother_label = 0;
2998         Int_t n_mother_steps = 0;
2999         
3000         // positive track
3001         const Int_t p_label = TMath::Abs(pTrackTPC->GetLabel());
3002         TParticle* p_mcTrack = fMCStack->Particle(p_label);         
3003         if (p_mcTrack){
3004           
3005           if(fMCStack->IsPhysicalPrimary(p_label))
3006             p_primaryFlag = 1;
3007
3008           //10/01/13. Add a flag to see if it is from material or WD
3009           
3010           if(fMCStack->IsSecondaryFromWeakDecay(p_label))
3011             p_primaryFlag = 2;
3012           
3013           if(fMCStack->IsSecondaryFromMaterial(p_label))
3014             p_primaryFlag = 3;
3015
3016
3017           
3018           Int_t p_pdgCode = p_mcTrack->GetPdgCode();
3019           p_pidCode = GetPidCode(p_pdgCode);
3020           
3021           p_ptMC      = p_mcTrack->Pt();
3022           
3023           p_mother_label = FindPrimaryMotherLabelV0(fMCStack, p_label, 
3024                                                   p_mother_steps);
3025           if(p_mother_label>0) {
3026             TParticle* p_mother = fMCStack->Particle(p_mother_label);
3027             p_pdgMother = p_mother->GetPdgCode();
3028           }
3029         }
3030         
3031         // negative track
3032         const Int_t n_label = TMath::Abs(nTrackTPC->GetLabel());
3033         TParticle* n_mcTrack = fMCStack->Particle(n_label);         
3034         if (n_mcTrack){
3035           
3036           if(fMCStack->IsPhysicalPrimary(n_label))
3037             n_primaryFlag = 1;
3038
3039           //10/01/13. Add a flag to see if it is from material or WD
3040           
3041           if(fMCStack->IsSecondaryFromWeakDecay(n_label))
3042             n_primaryFlag = 2;
3043           
3044           if(fMCStack->IsSecondaryFromMaterial(n_label))
3045             n_primaryFlag = 3;
3046
3047           
3048           Int_t n_pdgCode = n_mcTrack->GetPdgCode();
3049           n_pidCode = GetPidCode(n_pdgCode);
3050           
3051           n_ptMC      = n_mcTrack->Pt();
3052           
3053           n_mother_label = FindPrimaryMotherLabelV0(fMCStack, n_label, 
3054                                                   n_mother_steps);
3055           if(n_mother_label>0) {
3056             TParticle* n_mother = fMCStack->Particle(n_mother_label);
3057             n_pdgMother = n_mother->GetPdgCode();
3058           }
3059         }
3060         
3061         // Check if V0 is primary = first and the same mother of both partciles
3062         if(p_mother_label>0 && n_mother_label>0 && p_mother_label == n_mother_label) {
3063           pdgV0 = p_pdgMother;
3064           if(p_mother_steps == 1 && n_mother_steps == 1) {
3065             primaryV0 = 1;
3066           }
3067         }
3068       }
3069       
3070
3071       if(fTreeOption) {
3072         
3073         DeDxV0* v0datatpc = new((*fV0ArrayTPCPar)[nadded]) DeDxV0();
3074         nadded++;
3075         
3076         // v0 data
3077         v0datatpc->p       = esdV0->P();
3078         v0datatpc->pt      = esdV0->Pt();
3079         v0datatpc->eta     = esdV0->Eta();
3080         v0datatpc->phi     = esdV0->Phi();
3081         v0datatpc->pdca    = TMath::Sqrt(pdcaxy*pdcaxy + pdcaz*pdcaz);
3082         v0datatpc->ndca    = TMath::Sqrt(ndcaxy*ndcaxy + ndcaz*ndcaz);
3083         v0datatpc->dmassG  = deltaInvMassG;
3084         v0datatpc->dmassK0 = deltaInvMassK0s;
3085         v0datatpc->dmassL  = deltaInvMassL;
3086         v0datatpc->dmassAL = deltaInvMassAntiL;
3087         v0datatpc->alpha   = alpha;
3088         v0datatpc->ptarm   = ptarm;
3089         v0datatpc->decayr  = lV0Radius;
3090         v0datatpc->decayl  = lV0DecayLength;
3091         
3092         // New parameters
3093         v0datatpc->status  = esdV0->GetOnFlyStatus();
3094         v0datatpc->chi2    = esdV0->GetChi2V0();
3095         v0datatpc->cospt   = esdV0->GetV0CosineOfPointingAngle(); 
3096         // cospt: as I understand this means that the pointing to the vertex
3097         // is fine so I remove the dcaxy and dcaz for the V= class
3098         v0datatpc->dcadaughters = esdV0->GetDcaV0Daughters();
3099         v0datatpc->primary = primaryV0;
3100         v0datatpc->pdg     = pdgV0;
3101         
3102         // positive track
3103         v0datatpc->ptrack.p       = pp;
3104         v0datatpc->ptrack.pt      = ppt;
3105         //        v0data->ptrack.ptcon   = ppt_con;
3106         //        v0data->ptrack.tpcchi  = ptpcchi;
3107         v0datatpc->ptrack.eta     = peta;
3108         v0datatpc->ptrack.phi     = pphi;
3109         v0datatpc->ptrack.q       = pcharge;
3110         v0datatpc->ptrack.ncl     = pncl;
3111         v0datatpc->ptrack.neff    = pneff;
3112         v0datatpc->ptrack.dedx    = pdedx;
3113
3114         v0datatpc->ptrack.isTOFout   = 0;
3115         v0datatpc->ptrack.hasTOFtime = 0;
3116         v0datatpc->ptrack.isTOFmatched = 0;
3117         v0datatpc->ptrack.flength    =   0;
3118         v0datatpc->ptrack.ftimetof   = 0;
3119         v0datatpc->ptrack.exptoftimeel = 0;
3120         v0datatpc->ptrack.exptoftimemu = 0;
3121         v0datatpc->ptrack.exptoftimepi = 0;
3122         v0datatpc->ptrack.exptoftimeka = 0;
3123         v0datatpc->ptrack.exptoftimepr = 0;
3124
3125
3126         v0datatpc->ptrack.dcaxy   = pdcaxy;
3127         v0datatpc->ptrack.dcaz    = pdcaz;
3128         v0datatpc->ptrack.pid     = p_pidCode;
3129         v0datatpc->ptrack.primary = p_primaryFlag;
3130         v0datatpc->ptrack.pttrue  = p_ptMC;
3131         v0datatpc->ptrack.mother  = p_pdgMother;
3132         v0datatpc->ptrack.filter  = filterFlag_p;
3133         v0datatpc->ptrack.tpcnclS    = psharedtpcclusters;
3134
3135
3136         // negative track
3137         v0datatpc->ntrack.p       = np;
3138         v0datatpc->ntrack.pt      = npt;
3139         //        v0data->ntrack.ptcon   = npt_con;
3140         //        v0data->ntrack.tpcchi  = ntpcchi;
3141         v0datatpc->ntrack.eta     = neta;
3142         v0datatpc->ntrack.phi     = nphi;
3143         v0datatpc->ntrack.q       = ncharge;
3144         v0datatpc->ntrack.ncl     = nncl;
3145         v0datatpc->ntrack.neff    = nneff;
3146         v0datatpc->ntrack.dedx    = ndedx;
3147
3148         v0datatpc->ntrack.isTOFout   = 0;
3149         v0datatpc->ntrack.hasTOFtime = 0;
3150         v0datatpc->ntrack.isTOFmatched = 0;
3151         v0datatpc->ntrack.flength    =   0;
3152         v0datatpc->ntrack.ftimetof   = 0;
3153         v0datatpc->ntrack.exptoftimeel = 0;
3154         v0datatpc->ntrack.exptoftimemu = 0;
3155         v0datatpc->ntrack.exptoftimepi = 0;
3156         v0datatpc->ntrack.exptoftimeka = 0;
3157         v0datatpc->ntrack.exptoftimepr = 0;
3158
3159         v0datatpc->ntrack.dcaxy   = ndcaxy;
3160         v0datatpc->ntrack.dcaz    = ndcaz;
3161         v0datatpc->ntrack.pid     = n_pidCode;
3162         v0datatpc->ntrack.primary = n_primaryFlag;
3163         v0datatpc->ntrack.pttrue  = n_ptMC;
3164         v0datatpc->ntrack.mother  = n_pdgMother;
3165         v0datatpc->ntrack.filter  = filterFlag_n;
3166         v0datatpc->ntrack.tpcnclS    = nsharedtpcclusters;
3167       }
3168       
3169       // clean up loop over v0
3170       
3171       delete negPiKF;
3172       delete posPiKF;
3173       delete posPKF;
3174       delete negAPKF;
3175  
3176
3177
3178     }
3179
3180  
3181  
3182
3183
3184   }
3185   delete myPrimaryVertex;
3186
3187   
3188 }
3189 //__________________________________________________________________________
3190 void AliAnalysisTaskHighPtDeDx::ProduceArrayV0AOD( AliAODEvent *AODevent, AnalysisMode analysisMode ){
3191   Int_t nv0s = AODevent->GetNumberOfV0s();
3192   if(nv0s<1)return;
3193   //Int_t     trackmult = 0; // no pt cuts
3194   Int_t     nadded    = 0;
3195   
3196   AliAODVertex *myBestPrimaryVertex = AODevent->GetPrimaryVertex();
3197   if (!myBestPrimaryVertex) return;
3198   
3199   
3200   if( analysisMode == kGlobalTrk ){
3201     if(fV0ArrayGlobalPar)
3202       fV0ArrayGlobalPar->Clear();
3203   } else if( analysisMode == kTPCTrk ){
3204     if(fV0ArrayTPCPar)
3205       fV0ArrayTPCPar->Clear();
3206   }
3207   
3208   if( analysisMode==kGlobalTrk ){  
3209     
3210     
3211     // ################################
3212     // #### BEGINNING OF V0 CODE ######
3213     // ################################
3214     // This is the begining of the V0 loop  
3215     for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
3216       AliAODv0 *aodV0 = AODevent->GetV0(iV0);
3217       if (!aodV0) continue;
3218       
3219       // common part
3220       
3221       // AliAODTrack (V0 Daughters)
3222       AliAODVertex* vertex = aodV0->GetSecondaryVtx();
3223       if (!vertex) {
3224         Printf("ERROR: Could not retrieve vertex");
3225         continue;
3226       }
3227       
3228       AliAODTrack *pTrack = (AliAODTrack*)vertex->GetDaughter(0);
3229       AliAODTrack *nTrack = (AliAODTrack*)vertex->GetDaughter(1);
3230       if (!pTrack || !nTrack) {
3231         Printf("ERROR: Could not retrieve one of the daughter track");
3232         continue;
3233       }
3234       
3235       // Remove like-sign
3236       if (pTrack->Charge() == nTrack->Charge()) {
3237         //cout<< "like sign, continue"<< endl;
3238         continue;
3239       } 
3240       
3241       // Make sure charge ordering is ok
3242       if (pTrack->Charge() < 0) {
3243         AliAODTrack* helpTrack = pTrack;
3244         pTrack = nTrack;
3245         nTrack = helpTrack;
3246       } 
3247       
3248       // Eta cut on decay products
3249       if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)
3250         continue;
3251       
3252         //Pre-selection to reduce output size
3253         // Pt cut on decay products
3254         if (aodV0->Pt() < fMinPtV0) continue;
3255         // No point in keeping low cospa values...
3256         if (aodV0->CosPointingAngle(myBestPrimaryVertex) < 0.996 ) continue;
3257         //Reject on-the-fly tracks too
3258         if (aodV0->GetOnFlyStatus() != 0 ) continue;
3259       
3260
3261
3262       //check positive tracks
3263       UShort_t filterFlag_p = 0;
3264       
3265       if (fTrackFilterGolden) {  
3266         // For AOD068 
3267         
3268         if(pTrack->TestFilterBit(1024)) {
3269           filterFlag_p +=1;
3270         }
3271       }
3272       
3273       
3274       if (fTrackFilterTPC) {
3275         // TPC only cuts is bit 1 according to above macro
3276         // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX, 
3277         if(pTrack->TestFilterBit(1)){
3278           filterFlag_p +=2;       
3279         }
3280       }
3281
3282
3283       // if(pTrack->IsOn(AliAODTrack::kTPCrefit)) { //hljunggr: Tuvas flag for tpcrefit
3284       //        // Minimum number of clusters
3285       //   Float_t nCrossedRowsTPC = pTrack->GetTPCClusterInfo(2,1);
3286       //   if (nCrossedRowsTPC >= 70) {
3287       //          Int_t findable=pTrack->GetTPCNclsF();
3288       //          if (findable > 0){ 
3289       //            if (nCrossedRowsTPC/findable >= 0.8) 
3290       //              filterFlag_p += 16;
3291       //          }
3292       //        }
3293         
3294       // }
3295       
3296
3297       
3298       //check negative tracks
3299       UShort_t filterFlag_n = 0;
3300
3301       
3302       if (fTrackFilterGolden) {
3303         
3304         // ITSTPC2010 cuts is bit 32 according to above macro, new versions of aliroot includes the golden cuts
3305         if(nTrack->TestFilterBit(1024)) {
3306           filterFlag_n +=1;
3307         }
3308       }
3309       
3310       
3311       if (fTrackFilterTPC) {
3312         // TPC only cuts is bit 1 according to above macro
3313         // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX, 
3314         if(nTrack->TestFilterBit(1)){
3315           filterFlag_n +=2;
3316         }
3317       }
3318       
3319
3320       // if(nTrack->IsOn(AliAODTrack::kTPCrefit)) { //hljunggr: Tuvas flag for tpcrefit
3321       //        // Minimum number of clusters
3322       //   Float_t nCrossedRowsTPC = pTrack->GetTPCClusterInfo(2,1);
3323       //   if (nCrossedRowsTPC >= 70) {
3324       //          Int_t findable=nTrack->GetTPCNclsF();
3325       //          if (findable > 0){ 
3326       //            if (nCrossedRowsTPC/findable >= 0.8) filterFlag_n += 16;
3327       //          }
3328       //        }
3329       // }
3330       
3331       
3332       Float_t alpha = aodV0->AlphaV0();
3333       Float_t ptarm = aodV0->PtArmV0();
3334       // Double_t pVtxPos= v0->PrimaryVtxPosition();      
3335       
3336       Double_t lV0Radius      = aodV0->RadiusV0();
3337       Double_t lV0DecayLength = aodV0->DecayLength(myBestPrimaryVertex);
3338       
3339       Double_t deltaInvMassG     = aodV0->InvMass2Prongs(0,1,11,11);
3340       Double_t deltaInvMassK0s   = aodV0->MassK0Short()-0.498;
3341       Double_t deltaInvMassL     = aodV0->MassLambda()-1.116;
3342       Double_t deltaInvMassAntiL = aodV0->MassAntiLambda()-1.116;
3343       
3344       if(TMath::Abs(deltaInvMassK0s) > fMassCut &&
3345          TMath::Abs(deltaInvMassL) > fMassCut &&
3346          TMath::Abs(deltaInvMassAntiL) > fMassCut &&
3347          TMath::Abs(deltaInvMassG) > fMassCut )
3348         continue;
3349
3350
3351       // TODO: Why should these be different? Different mass hypothesis = energy loss
3352       // This is not important for us as we focus on the decay products!
3353       // Double_t ptK0s        = v0K0sKF.GetPt(); 
3354       // Double_t ptL          = v0LambdaKF.GetPt();
3355       // Double_t ptAntiL      = v0AntiLambdaKF.GetPt();     
3356       
3357       // Extract track information
3358       
3359       Double_t b[2], cov[3];
3360       if(!pTrack->PropagateToDCA(myBestPrimaryVertex, AODevent->GetMagneticField(), kVeryBig, b, cov))
3361         filterFlag_p += 32; // propagation failed!!!!!
3362       
3363       Float_t pdcaxy   = b[0];
3364       Float_t pdcaz    = b[1];
3365       if(!nTrack->PropagateToDCA(myBestPrimaryVertex, AODevent->GetMagneticField(), kVeryBig, b, cov))
3366         filterFlag_n += 32; // propagation failed!!!!!
3367       Float_t ndcaxy   = b[0];
3368       Float_t ndcaz    = b[1];
3369       
3370       Short_t pcharge  = pTrack->Charge();
3371       Float_t ppt      = pTrack->Pt();
3372       Float_t pp       = pTrack->P(); 
3373       Float_t peta     = pTrack->Eta();
3374       Float_t pphi     = pTrack->Phi();
3375       //        Float_t ptpcchi  = pTrack->Chi2perNDF();
3376       
3377       AliAODPid* pPid = pTrack->GetDetPid();
3378       Short_t pncl     = -10;
3379       Short_t pneff    = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
3380       Float_t pdedx    = -10;
3381    
3382       Bool_t pIsTOFout=kFALSE;
3383       Bool_t pHasTOFTime=kTRUE;
3384       Bool_t pIsMatched=kFALSE;
3385       Float_t plengthtrack=-999;//in aod we do not have that information, beta must be: beta=inttime/timeTOF 
3386       Float_t ptimeTOF=-999;
3387       Float_t pfexptimeel=-999;
3388       Float_t pfexptimemu=-999;
3389       Float_t pfexptimepi=-999;
3390       Float_t pfexptimeka=-999;
3391       Float_t pfexptimepr=-999;  
3392
3393
3394       if(pPid) {
3395         pncl     = pPid->GetTPCsignalN();
3396         pdedx    = pPid->GetTPCsignal();
3397         //TOF
3398
3399         /*
3400         if (pTrack->GetStatus()&AliESDtrack::kTOFpid){
3401           Double_t tof[5];
3402           pPid->GetIntegratedTimes(tof);
3403           pbeta = tof[0]/pPid->GetTOFsignal();
3404         }
3405         */
3406
3407
3408         if ((pTrack->GetStatus()&AliESDtrack::kTOFout)==0)
3409           pIsTOFout=kTRUE;
3410         
3411
3412         if ((pTrack->GetStatus()&AliESDtrack::kTIME)==0)
3413           pHasTOFTime=kFALSE;
3414         
3415
3416         if (!(pTrack->GetStatus()&AliESDtrack::kTOFmismatch)==0)
3417           pIsMatched=kFALSE;
3418         else
3419           pIsMatched=kTRUE;
3420         
3421         ptimeTOF=pPid->GetTOFsignal();
3422         
3423         Double_t pinttime[5]={0,0,0,0,0}; 
3424         pPid->GetIntegratedTimes(pinttime);// Returns the array with integrated times for each particle hypothesis
3425         pfexptimeel=pinttime[0];
3426         pfexptimemu=pinttime[1];
3427         pfexptimepi=pinttime[2];
3428         pfexptimeka=pinttime[3];
3429         pfexptimepr=pinttime[4];  
3430
3431
3432
3433
3434
3435       }
3436       TBits psharedTPC=pTrack->GetTPCSharedMap();
3437       Int_t psharedtpcclusters=psharedTPC.CountBits(0)-psharedTPC.CountBits(159);
3438
3439       TBits nsharedTPC=nTrack->GetTPCSharedMap();
3440       Int_t nsharedtpcclusters=nsharedTPC.CountBits(0)-nsharedTPC.CountBits(159);
3441
3442
3443       
3444       Short_t ncharge  = nTrack->Charge();
3445       Float_t npt      = nTrack->Pt();
3446       Float_t np       = nTrack->P(); 
3447       Float_t neta     = nTrack->Eta();
3448       Float_t nphi     = nTrack->Phi();
3449       //        Float_t ntpcchi  = nTrack->Chi2perNDF();
3450       
3451       AliAODPid* nPid = nTrack->GetDetPid();
3452       Short_t nncl     = -10;
3453       Short_t nneff    = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
3454       Float_t ndedx    = -10;
3455       //Float_t nbeta = -99;
3456
3457       Bool_t nIsTOFout=kFALSE;
3458       Bool_t nHasTOFTime=kTRUE;
3459       Bool_t nIsMatched=kFALSE;
3460       Float_t nlengthtrack=-999;//in aod we do not have that information, beta must be: beta=inttime/timeTOF 
3461       Float_t ntimeTOF=-999;
3462       Float_t nfexptimeel=-999;
3463       Float_t nfexptimemu=-999;
3464       Float_t nfexptimepi=-999;
3465       Float_t nfexptimeka=-999;
3466       Float_t nfexptimepr=-999;  
3467
3468
3469       if(pPid) {
3470         nncl     = nPid->GetTPCsignalN();
3471         ndedx    = nPid->GetTPCsignal();
3472         //TOF
3473         /*
3474         if (nTrack->GetStatus()&AliESDtrack::kTOFpid){
3475           Double_t tof[5];
3476           nPid->GetIntegratedTimes(tof);
3477           nbeta = tof[0]/nPid->GetTOFsignal();
3478           }*/
3479
3480
3481         if ((nTrack->GetStatus()&AliESDtrack::kTOFout)==0)
3482           nIsTOFout=kTRUE;
3483         
3484
3485         if ((nTrack->GetStatus()&AliESDtrack::kTIME)==0)
3486           nHasTOFTime=kFALSE;
3487         
3488         if (!(nTrack->GetStatus()&AliESDtrack::kTOFmismatch)==0)
3489           nIsMatched=kFALSE;
3490         else
3491           nIsMatched=kTRUE;
3492        
3493         ntimeTOF=nPid->GetTOFsignal();
3494         
3495         Double_t ninttime[5]={0,0,0,0,0}; 
3496         nPid->GetIntegratedTimes(ninttime);// Returns the array with integrated times for each particle hypothesis
3497         nfexptimeel=ninttime[0];
3498         nfexptimemu=ninttime[1];
3499         nfexptimepi=ninttime[2];
3500         nfexptimeka=ninttime[3];
3501         nfexptimepr=ninttime[4];  
3502
3503       }
3504       
3505       Int_t   primaryV0     = 0; // 0 means that the tracks are not both daughters of a primary particle (1 means they are)
3506       Int_t   pdgV0         = 0; // 0 means that they don't have same origin for MC (1 means they have the same original mother)
3507       Float_t p_ptMC        = 0;
3508       Short_t p_pidCode     = 0; // 0 = real data / no mc track!
3509       Short_t p_primaryFlag = 0; // 0 = real data / not primary mc track  
3510       Int_t   p_pdgMother   = 0;
3511       Float_t n_ptMC        = 0;
3512       Short_t n_pidCode     = 0; // 0 = real data / no mc track!
3513       Short_t n_primaryFlag = 0; // 0 = real data / not primary mc track  
3514       Int_t   n_pdgMother   = 0;
3515       if(fAnalysisMC) {
3516         
3517         AliAODMCParticle* p_mother = 0;
3518         Int_t p_mother_steps = 0;
3519         AliAODMCParticle* n_mother = 0;
3520         Int_t n_mother_steps = 0;
3521         
3522         // positive track
3523         const Int_t p_label = TMath::Abs(pTrack->GetLabel());
3524         
3525         AliAODMCParticle* p_mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(p_label));
3526         if (p_mcTrack){
3527           
3528           if(p_mcTrack->IsPhysicalPrimary())
3529             p_primaryFlag = 1;
3530           
3531
3532