Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / grid / AliAnalysisTaskHighPtDeDxV0.cxx
1 #include "AliAnalysisTaskHighPtDeDxV0.h"
2
3 // ROOT includes
4 #include <TList.h>
5 #include <TTree.h>
6 #include <TMath.h>
7 #include <TH1.h>
8 #include <TParticle.h>
9 #include <TFile.h>
10
11 // AliRoot includes
12 #include <AliAnalysisManager.h>
13 #include <AliAnalysisFilter.h>
14 #include <AliESDInputHandler.h>
15 #include <AliESDEvent.h>
16 #include <AliESDVertex.h>
17 #include <AliLog.h>
18 #include <AliExternalTrackParam.h>
19
20 #include <AliMCEventHandler.h>
21 #include <AliMCEvent.h>
22 #include <AliStack.h>
23
24 #include <AliHeader.h>
25 #include <AliGenPythiaEventHeader.h>
26 #include <AliGenDPMjetEventHeader.h>
27
28 #include <AliAODVertex.h>
29 #include <AliESDv0.h>
30 #include <AliKFVertex.h>
31
32 #include "AliCentrality.h" 
33 #include "AliESDtrackCuts.h"
34
35 #include <AliAODTrack.h> 
36 #include <AliAODPid.h> 
37 #include <AliAODMCHeader.h> 
38
39 // STL includes
40 #include <iostream>
41 using namespace std;
42
43 //_____________________________________________________________________________
44 //
45 // Responsible:
46 // Peter Christiansen (Lund)
47 //
48 //_____________________________________________________________________________
49
50
51 ClassImp(AliAnalysisTaskHighPtDeDxV0)
52
53 const Double_t AliAnalysisTaskHighPtDeDxV0::fgkClight = 2.99792458e-2;
54
55 //_____________________________________________________________________________
56 AliAnalysisTaskHighPtDeDxV0::AliAnalysisTaskHighPtDeDxV0():
57   AliAnalysisTaskSE(),
58   fESD(0x0),
59   fAOD(0x0),
60   fMC(0x0),
61   fMCStack(0x0),
62   fMCArray(0x0),
63   fTrackFilter(0x0),
64   fTrackFilterGolden(0x0),
65   fTrackFilterTPC(0x0),
66   fAnalysisType("ESD"),
67   fAnalysisMC(kFALSE),
68   fAnalysisPbPb(kFALSE),
69   ftrigBit1(0x0),
70   ftrigBit2(0x0),
71   fEvent(0x0),
72   fV0ArrayGlobalPar(0x0),
73   fV0ArrayTPCPar(0x0),
74   fTrackArrayMC(0x0),
75   fVtxCut(10.0),  
76   fEtaCut(0.9),  
77   fMinPt(0.1),
78   fMinCent(0.0),
79   fMaxCent(100.0),
80   fMassCut(0.1),
81   fTreeOption(0),
82   fRequireRecV0(kTRUE),
83   fStoreMcIn(kFALSE),
84   fMcProcessType(-999),
85   fTriggeredEventMB(-999),
86   fVtxStatus(-999),
87   fZvtx(-999),
88   fZvtxMC(-999),
89   fRun(-999),
90   fEventId(-999),
91   fListOfObjects(0), 
92   fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),
93   fTree(0x0),
94   fnv0(0)
95 {
96   // Default constructor (should not be used)
97
98 }
99
100 //______________________________________________________________________________
101 AliAnalysisTaskHighPtDeDxV0::AliAnalysisTaskHighPtDeDxV0(const char *name):
102   AliAnalysisTaskSE(name),
103   fESD(0x0),
104   fAOD(0x0),
105   fMC(0x0),
106   fMCStack(0x0),
107   fMCArray(0x0),
108   fTrackFilter(0x0),
109   fTrackFilterGolden(0x0),
110   fTrackFilterTPC(0x0),
111   fAnalysisType("ESD"),
112   fAnalysisMC(kFALSE),
113   fAnalysisPbPb(kFALSE),
114   ftrigBit1(0x0),
115   ftrigBit2(0x0),
116   fEvent(0x0),
117   fV0ArrayGlobalPar(0x0),
118   fV0ArrayTPCPar(0x0),
119   fTrackArrayMC(0x0),
120   fVtxCut(10.0),  
121   fEtaCut(0.9),  
122   fMinPt(0.1),
123   fMinCent(0.0),
124   fMaxCent(100.0),
125   fMassCut(0.1),
126   fTreeOption(0),
127   fRequireRecV0(kTRUE),
128   fStoreMcIn(kFALSE),
129   fMcProcessType(-999),
130   fTriggeredEventMB(-999),
131   fVtxStatus(-999),
132   fZvtx(-999),
133   fZvtxMC(-999),
134   fRun(-999),
135   fEventId(-999),
136   fListOfObjects(0), 
137   fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),
138   fTree(0x0),
139   fnv0(0)
140 {
141   // Output slot #1 writes into a TList
142   DefineOutput(1, TList::Class());
143 }
144
145 //_____________________________________________________________________________
146 AliAnalysisTaskHighPtDeDxV0::~AliAnalysisTaskHighPtDeDxV0()
147 {
148   // Destructor
149   // histograms are in the output list and deleted when the output
150   // list is deleted by the TSelector dtor
151   if (fListOfObjects && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
152     delete fListOfObjects;
153     fListOfObjects = 0;
154   }
155   
156   // //for proof running; I cannot create tree do to memory limitations -> work with THnSparse 
157   // if (fListOfObjects  && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutputList;
158   
159   
160   
161 }
162
163 //______________________________________________________________________________
164 void AliAnalysisTaskHighPtDeDxV0::UserCreateOutputObjects()
165
166   // This method is called once per worker node
167   // Here we define the output: histograms and debug tree if requested 
168
169   OpenFile(1);
170   fListOfObjects = new TList();
171   fListOfObjects->SetOwner();
172   
173   //
174   // Histograms
175   //  
176   fEvents = new TH1I("fEvents","Number of analyzed events; Events; Counts", 1, 0, 1);
177   fListOfObjects->Add(fEvents);
178   
179   fVtx = new TH1I("fVtx","Vtx info (0=no, 1=yes); Vtx; Counts", 2, -0.5, 1.5);
180   fListOfObjects->Add(fVtx);
181
182   fnv0 = new TH1F("fnv0","# V0's", 10000,0,10000);
183   fListOfObjects->Add(fnv0);
184
185   fVtxBeforeCuts = new TH1F("fVtxBeforeCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
186   fListOfObjects->Add(fVtxBeforeCuts);
187   
188   fVtxAfterCuts = new TH1F("fVtxAfterCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
189   fListOfObjects->Add(fVtxAfterCuts);
190
191   if (fAnalysisMC) {    
192     fVtxMC = new TH1I("fVtxMC","Vtx info - no trigger cut (0=no, 1=yes); Vtx; Counts", 2, -0.5, 1.5);
193     fListOfObjects->Add(fVtxMC);
194   
195   }
196
197   if (fTreeOption) {
198
199     fTree = new TTree("tree","Event data");
200     fEvent = new DeDxEvent();
201     fTree->Branch("event", &fEvent);
202
203     fV0ArrayGlobalPar = new TClonesArray("DeDxV0", 1000);
204     fTree->Bronch("v0GlobalPar", "TClonesArray", &fV0ArrayGlobalPar);
205
206     fV0ArrayTPCPar = new TClonesArray("DeDxV0", 1000);
207     fTree->Bronch("v0TPCPar", "TClonesArray", &fV0ArrayTPCPar);
208
209     if (fAnalysisMC && fStoreMcIn) {    
210       fTrackArrayMC = new TClonesArray("DeDxTrackMC", 1000);
211       fTree->Bronch("trackMC", "TClonesArray", &fTrackArrayMC);
212     }
213
214     fTree->SetDirectory(0);
215
216     fListOfObjects->Add(fTree);
217
218   }
219
220   // Post output data.
221   PostData(1, fListOfObjects);
222 }
223
224 //______________________________________________________________________________
225 void AliAnalysisTaskHighPtDeDxV0::UserExec(Option_t *) 
226 {
227   // Main loop
228   
229   //
230   // Comment: This method matches completely the same method for the high pT
231   // tracks
232   //
233
234
235   //
236   // First we make sure that we have valid input(s)!
237   //
238   if (fAnalysisType == "ESD"){
239     fESD = dynamic_cast<AliESDEvent*>(InputEvent());
240     if(!fESD){
241       Printf("%s:%d ESDEvent not found in Input Manager",(char*)__FILE__,__LINE__);
242       this->Dump();
243       return;
244     }    
245   } else {
246     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
247     if(!fAOD){
248       Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
249       this->Dump();
250       return;
251     }    
252   }
253
254   if (fAnalysisMC) {
255
256     if (fAnalysisType == "ESD"){
257       fMC = dynamic_cast<AliMCEvent*>(MCEvent());
258       if(!fMC){
259         Printf("%s:%d MCEvent not found in Input Manager",(char*)__FILE__,__LINE__);
260         this->Dump();
261         return;
262       }    
263
264       fMCStack = fMC->Stack();
265       
266       if(!fMCStack){
267         Printf("%s:%d MCStack not found in Input Manager",(char*)__FILE__,__LINE__);
268         this->Dump();
269         return;
270       }    
271     } else { // AOD
272
273       fMC = dynamic_cast<AliMCEvent*>(MCEvent());
274       if(fMC)
275         fMC->Dump();
276
277       fMCArray = (TClonesArray*)fAOD->FindListObject("mcparticles");
278       if(!fMCArray){
279         Printf("%s:%d AOD MC array not found in Input Manager",(char*)__FILE__,__LINE__);
280         this->Dump();
281         return;
282       }    
283     }
284   }
285
286     // Get trigger decision
287   fTriggeredEventMB = 0; //init
288  
289
290
291   if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
292      ->IsEventSelected() & ftrigBit1 ){
293     fTriggeredEventMB = 1;  //event triggered as minimum bias
294   }
295   if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
296      ->IsEventSelected() & ftrigBit2 ){
297     // From AliVEvent:
298     //    kINT7         = BIT(1), // V0AND trigger, offline V0 selection
299     fTriggeredEventMB += 2;  
300   }
301
302   
303   // Get process type for MC
304   fMcProcessType = 0; // -1=invalid, 0=data, 1=ND, 2=SD, 3=DD
305
306   // real data that are not triggered we skip
307   if(!fAnalysisMC && !fTriggeredEventMB)
308     return; 
309   
310   if (fAnalysisMC) {
311     
312     if (fAnalysisType == "ESD"){
313
314       AliHeader* headerMC = fMC->Header();
315       if (headerMC) {
316         
317         AliGenEventHeader* genHeader = headerMC->GenEventHeader();
318         TArrayF vtxMC(3); // primary vertex  MC 
319         vtxMC[0]=9999; vtxMC[1]=9999;  vtxMC[2]=9999; //initialize with dummy
320         if (genHeader) {
321           genHeader->PrimaryVertex(vtxMC);
322         }
323         fZvtxMC = vtxMC[2];
324         
325         // PYTHIA:
326         AliGenPythiaEventHeader* pythiaGenHeader =
327           dynamic_cast<AliGenPythiaEventHeader*>(headerMC->GenEventHeader());
328         if (pythiaGenHeader) {  //works only for pythia
329           fMcProcessType =  GetPythiaEventProcessType(pythiaGenHeader->ProcessType());
330         }
331         // PHOJET:
332         AliGenDPMjetEventHeader* dpmJetGenHeader =
333           dynamic_cast<AliGenDPMjetEventHeader*>(headerMC->GenEventHeader());
334         if (dpmJetGenHeader) {
335           fMcProcessType = GetDPMjetEventProcessType(dpmJetGenHeader->ProcessType());
336         }
337       }
338     } else { // AOD
339       
340       AliAODMCHeader* mcHeader = dynamic_cast<AliAODMCHeader*>(fAOD->FindListObject("mcHeader")); 
341       if(mcHeader) {
342         fZvtxMC = mcHeader->GetVtxZ();
343         
344         if(strstr(mcHeader->GetGeneratorName(), "Pythia")) {
345           fMcProcessType =  GetPythiaEventProcessType(mcHeader->GetEventType());
346         } else {
347           fMcProcessType =  GetDPMjetEventProcessType(mcHeader->GetEventType());
348         }
349       }
350     }
351   }
352   
353   // There are 3 cases
354   // Vertex: NO  - status -1
355   // Vertex: YES : outside cut - status 0
356   //             : inside cut  - status 1
357   // We have to be careful how we normalize because we probably want to
358   // normalize to:
359   // Nevents=(No vertex + outside + inside)/(out + in)*in
360   
361   if (fAnalysisType == "ESD")
362     fZvtx = GetVertex(fESD);
363   else // AOD
364     fZvtx = GetVertex(fAOD);
365     
366   fVtxStatus = -999;
367   
368   if(fZvtx<-990) {
369     
370     fVtxStatus = -1;
371     if(fTriggeredEventMB)
372       fVtx->Fill(0);
373     if(fAnalysisMC)
374       fVtxMC->Fill(0);
375   } else {
376     
377     if(fTriggeredEventMB)
378       fVtx->Fill(1);
379     if(fAnalysisMC)
380       fVtxMC->Fill(1);
381     fVtxBeforeCuts->Fill(fZvtx);
382     fVtxStatus = 0;
383     if (TMath::Abs(fZvtx) < fVtxCut) {  
384       fVtxAfterCuts->Fill(fZvtx);
385       fVtxStatus = 1;
386     }
387   }
388   
389   // store MC event data no matter what
390   if(fAnalysisMC && fStoreMcIn) {
391
392     if (fAnalysisType == "ESD") {
393       ProcessMCTruthESD();
394     } else { // AOD
395       ProcessMCTruthAOD();
396     }
397   }      
398   Float_t centrality = -10;
399   // only analyze triggered events
400   if(fTriggeredEventMB) {
401     
402     if (fAnalysisType == "ESD"){
403       if(fAnalysisPbPb){
404         AliCentrality *centObject = fESD->GetCentrality();
405         centrality = centObject->GetCentralityPercentile("V0M"); 
406         if((centrality>fMaxCent)||(centrality<fMinCent))return;
407       }
408       Int_t nv0test=fESD->GetNumberOfV0s();
409       cout<<"&&&&&&&&&&&&&&&&   hello world"<<endl;
410       fnv0->Fill(nv0test);
411       AnalyzeESD(fESD);
412       
413     } else { // AOD
414       if(fAnalysisPbPb){
415         AliCentrality *centObject = fAOD->GetCentrality();
416         centrality = centObject->GetCentralityPercentile("V0M"); 
417         if((centrality>fMaxCent)||(centrality<fMinCent))return;
418       }
419       Int_t nv0test=fAOD->GetNumberOfV0s();
420       fnv0->Fill(nv0test);
421       AnalyzeAOD(fAOD);
422       
423     }
424   }
425   
426
427
428   if( fTreeOption) {
429
430     fEvent->process = fMcProcessType;
431     fEvent->trig    = fTriggeredEventMB;
432     fEvent->zvtxMC  = fZvtxMC;
433     fEvent->cent      = centrality;
434
435     if(!fRequireRecV0 || fEvent->n > 0) // only fill if there are accepted V0s
436       fTree->Fill();
437     fV0ArrayGlobalPar->Clear();
438     fV0ArrayTPCPar->Clear();
439     if (fAnalysisMC && fStoreMcIn)
440       fTrackArrayMC->Clear();
441   }
442   
443   // Post output data.
444   PostData(1, fListOfObjects);
445 }
446
447 //________________________________________________________________________
448 void AliAnalysisTaskHighPtDeDxV0::AnalyzeESD(AliESDEvent* esdEvent)
449 {
450   fRun  = esdEvent->GetRunNumber();
451   fEventId = 0;
452   if(esdEvent->GetHeader())
453     fEventId = GetEventIdAsLong(esdEvent->GetHeader());
454
455   //  Int_t     event     = esdEvent->GetEventNumberInFile();
456   UInt_t    time      = esdEvent->GetTimeStamp();
457   //  ULong64_t trigger   = esdEvent->GetTriggerMask();
458   Float_t   magf      = esdEvent->GetMagneticField();
459
460
461
462   Short_t isPileup = esdEvent->IsPileupFromSPD();
463
464   // centrality
465   Float_t centrality = 120;
466   AliCentrality *centObject = esdEvent->GetCentrality();
467   if(centObject)
468     centrality = centObject->GetCentralityPercentile("V0M"); 
469
470   if(fTriggeredEventMB) {// Only MC case can we have not triggered events
471     
472     if(fVtxStatus==1) { // accepted vertex
473
474       // accepted event
475       fEvents->Fill(0);
476       
477       ProduceArrayTrksESD( esdEvent, kGlobalTrk );//produce array with global track parameters
478       ProduceArrayTrksESD( esdEvent, kTPCTrk );//array with tpc track parametes
479
480
481     } // end if vtx status
482   } // end if triggered
483   
484   if(fTreeOption) {
485
486     fEvent->run       = fRun;
487     fEvent->eventid   = fEventId;
488     fEvent->time      = time;
489     //fEvent->cent      = centrality;
490     fEvent->mag       = magf;
491     fEvent->zvtx      = fZvtx;
492     fEvent->vtxstatus = fVtxStatus;
493     //fEvent->trackmult = trackmult;
494     //fEvent->n         = nadded;
495     fEvent->pileup    = isPileup;
496   }
497 }
498
499 //________________________________________________________________________
500 void AliAnalysisTaskHighPtDeDxV0::AnalyzeAOD(AliAODEvent* aodEvent)
501 {
502   fRun  = aodEvent->GetRunNumber();
503   fEventId = 0;
504   if(aodEvent->GetHeader())
505     fEventId = GetEventIdAsLong(aodEvent->GetHeader());
506
507   //  Int_t     event     = aodEvent->GetEventNumberInFile();
508   UInt_t    time      = 0; // Missing AOD info? aodEvent->GetTimeStamp();
509   //  ULong64_t trigger   = aodEvent->GetTriggerMask();
510   Float_t   magf      = aodEvent->GetMagneticField();
511
512   Int_t     trackmult = 0; // no pt cuts
513   Int_t     nadded    = 0;
514
515   Short_t isPileup = aodEvent->IsPileupFromSPD();
516
517   // centrality
518   Float_t centrality = 120;
519   AliCentrality *centObject = aodEvent->GetCentrality();
520   if(centObject)
521     centrality = centObject->GetCentralityPercentile("V0M"); 
522
523   if(fTriggeredEventMB) {// Only MC case can we have not triggered events
524     
525     if(fVtxStatus==1) { // accepted vertex
526
527       // accepted event
528       fEvents->Fill(0);
529       
530       ProduceArrayTrksAOD( aodEvent, kGlobalTrk );//produce array with global track parameters
531       ProduceArrayTrksAOD( aodEvent, kTPCTrk );//array with tpc track parametes
532
533
534     } // end if vtx status
535   } // end if triggered
536   
537   if(fTreeOption) {
538
539     fEvent->run       = fRun;
540     fEvent->eventid   = fEventId;
541     fEvent->time      = time;
542     //fEvent->cent      = centrality;
543     fEvent->mag       = magf;
544     fEvent->zvtx      = fZvtx;
545     fEvent->vtxstatus = fVtxStatus;
546     //fEvent->trackmult = trackmult;
547     //fEvent->n         = nadded;
548     fEvent->pileup    = isPileup;
549   }
550
551
552
553
554
555
556
557 }
558
559 //_____________________________________________________________________________
560 Float_t AliAnalysisTaskHighPtDeDxV0::GetVertex(const AliVEvent* event) const
561 {
562   Float_t zvtx = -999;
563   
564   const AliVVertex* primaryVertex = event->GetPrimaryVertex(); 
565   
566   if(primaryVertex->GetNContributors()>0)
567     zvtx = primaryVertex->GetZ();
568
569   return zvtx;
570 }
571
572 //_____________________________________________________________________________
573 Short_t AliAnalysisTaskHighPtDeDxV0::GetPidCode(Int_t pdgCode) const 
574 {
575   // return our internal code for pions, kaons, and protons
576   
577   Short_t pidCode = 6;
578   
579   switch (TMath::Abs(pdgCode)) {
580   case 211:
581     pidCode = 1; // pion
582     break;
583   case 321:
584     pidCode = 2; // kaon
585     break;
586   case 2212:
587     pidCode = 3; // proton
588     break;
589   case 11:
590     pidCode = 4; // electron
591     break;
592   case 13:
593     pidCode = 5; // proton
594     break;
595   default:
596     pidCode = 6;  // something else?
597   };
598   
599   return pidCode;
600 }
601
602
603 //_____________________________________________________________________________
604 void AliAnalysisTaskHighPtDeDxV0::ProcessMCTruthESD() 
605 {
606   // Fill the special MC histogram with the MC truth info
607   
608   Short_t trackmult = 0;
609   Short_t nadded    = 0;
610   const Int_t nTracksMC = fMCStack->GetNtrack();
611   
612   for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
613     
614     //Cuts
615     if(!(fMCStack->IsPhysicalPrimary(iTracks)))
616       continue;
617     
618     TParticle* trackMC = fMCStack->Particle(iTracks);
619     
620     Double_t chargeMC = trackMC->GetPDG()->Charge();
621     if (chargeMC != 0) // select charge 0 particles!
622       continue;
623     
624     if (TMath::Abs(trackMC->Eta()) > fEtaCut )
625       continue;
626     
627     trackmult++;
628     
629     //   "charge:pt:p:eta:phi:pidCode"
630     Float_t ptMC      = trackMC->Pt();
631     Float_t pMC       = trackMC->P();
632     Float_t etaMC     = trackMC->Eta();
633     Float_t phiMC     = trackMC->Phi();
634     
635     Int_t pdgCode = trackMC->GetPdgCode();
636     if(!(pdgCode==310 || pdgCode == 3122 || pdgCode == -3122))
637       continue;
638     Short_t pidCodeMC = 0;
639     pidCodeMC = GetPidCode(pdgCode);
640     
641     // Warning: In the data code we cut on the daughters.....
642     if (trackMC->Pt() < fMinPt)
643       continue;
644     
645     if(fTreeOption) {
646       
647       DeDxTrackMC* track = new((*fTrackArrayMC)[nadded]) DeDxTrackMC();
648       nadded++;
649       
650       track->pMC   = pMC;
651       track->ptMC  = ptMC;
652       track->etaMC = etaMC;
653       track->phiMC = phiMC;
654       track->qMC   = Short_t(chargeMC);
655       track->pidMC = pidCodeMC;
656       track->pdgMC = pdgCode;
657     }
658     
659   }//MC track loop
660   
661   if(fTreeOption) {
662     
663     Sort(fTrackArrayMC, kTRUE);
664
665     fEvent->trackmultMC = trackmult;
666     fEvent->nMC         = nadded;
667   }
668   
669 }
670
671 //_____________________________________________________________________________
672 void AliAnalysisTaskHighPtDeDxV0::ProcessMCTruthAOD() 
673 {
674   // Fill the special MC histogram with the MC truth info
675
676   Short_t trackmult = 0;
677   Short_t nadded    = 0;
678   const Int_t nTracksMC = fMCArray->GetEntriesFast();
679   
680   for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
681     
682     AliAODMCParticle* trackMC = dynamic_cast<AliAODMCParticle*>(fMCArray->At(iTracks));
683     
684     //Cuts
685     if(!(trackMC->IsPhysicalPrimary()))
686       continue;
687     
688     
689     Double_t chargeMC = trackMC->Charge();
690     if (chargeMC != 0) // Keep the neutral particles
691       continue;
692     
693     if (TMath::Abs(trackMC->Eta()) > fEtaCut )
694       continue;
695     
696     trackmult++;
697     
698     //   "charge:pt:p:eta:phi:pidCode"
699     Float_t ptMC      = trackMC->Pt();
700     Float_t pMC       = trackMC->P();
701     Float_t etaMC     = trackMC->Eta();
702     Float_t phiMC     = trackMC->Phi();
703     
704     Int_t pdgCode = trackMC->PdgCode();
705     if(!(pdgCode==310 || pdgCode == 3122 || pdgCode == -3122))
706       continue;
707     Short_t pidCodeMC = 0;
708     pidCodeMC = GetPidCode(pdgCode);
709     
710     if (trackMC->Pt() < fMinPt)
711       continue;
712     
713     if(fTreeOption) {
714       
715       DeDxTrackMC* track = new((*fTrackArrayMC)[nadded]) DeDxTrackMC();
716       nadded++;
717       
718       track->pMC   = pMC;
719       track->ptMC  = ptMC;
720       track->etaMC = etaMC;
721       track->phiMC = phiMC;
722       track->qMC   = Short_t(chargeMC);
723       track->pidMC = pidCodeMC;
724       track->pdgMC = pdgCode;
725     }
726     
727   }//MC track loop
728   
729   if(fTreeOption) {
730     
731     Sort(fTrackArrayMC, kTRUE);
732
733     fEvent->trackmultMC = trackmult;
734     fEvent->nMC         = nadded;
735   }
736   
737 }
738
739 //_____________________________________________________________________________
740 Short_t AliAnalysisTaskHighPtDeDxV0::GetPythiaEventProcessType(Int_t pythiaType) {
741   //
742   // Get the process type of the event.  PYTHIA
743   //
744   // source PWG0   dNdpt 
745
746   Short_t globalType = -1; //init
747       
748   if(pythiaType==92||pythiaType==93){
749     globalType = 2; //single diffractive
750   }
751   else if(pythiaType==94){
752     globalType = 3; //double diffractive
753   }
754   //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB??
755   else {
756     globalType = 1;  //non diffractive
757   }
758   return globalType;
759 }
760
761 //_____________________________________________________________________________
762 Short_t AliAnalysisTaskHighPtDeDxV0::GetDPMjetEventProcessType(Int_t dpmJetType) {
763   //
764   // get the process type of the event.  PHOJET
765   //
766   //source PWG0   dNdpt 
767   // can only read pythia headers, either directly or from cocktalil header
768   Short_t globalType = -1;
769   
770   if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction
771     globalType = 1;
772   }
773   else if (dpmJetType==5 || dpmJetType==6) {
774     globalType = 2;
775   }
776   else if (dpmJetType==7) {
777     globalType = 3;
778   }
779   return globalType;
780 }
781
782 //_____________________________________________________________________________
783 ULong64_t AliAnalysisTaskHighPtDeDxV0::GetEventIdAsLong(AliVHeader* header) const
784 {
785   // To have a unique id for each event in a run!
786   // Modified from AliRawReader.h
787   return ((ULong64_t)header->GetBunchCrossNumber()+
788           (ULong64_t)header->GetOrbitNumber()*3564+
789           (ULong64_t)header->GetPeriodNumber()*16777215*3564);
790 }
791
792
793 //____________________________________________________________________
794 TParticle* AliAnalysisTaskHighPtDeDxV0::FindPrimaryMother(AliStack* stack, Int_t label)
795 {
796   //
797   // Finds the first mother among the primary particles of the particle identified by <label>,
798   // i.e. the primary that "caused" this particle
799   //
800   // Taken from AliPWG0Helper class
801   //
802
803   Int_t nSteps = 0;
804
805   Int_t motherLabel = FindPrimaryMotherLabel(stack, label, nSteps);
806   if (motherLabel < 0)
807     return 0;
808
809   return stack->Particle(motherLabel);
810 }
811
812 //____________________________________________________________________
813 Int_t AliAnalysisTaskHighPtDeDxV0::FindPrimaryMotherLabel(AliStack* stack, Int_t label, Int_t& nSteps)
814 {
815   //
816   // Finds the first mother among the primary particles of the particle identified by <label>,
817   // i.e. the primary that "caused" this particle
818   //
819   // returns its label
820   //
821   // Taken from AliPWG0Helper class
822   //
823   nSteps = 0;
824   const Int_t nPrim  = stack->GetNprimary();
825   
826   while (label >= nPrim) {
827
828     //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
829
830     nSteps++; // 1 level down
831     
832     TParticle* particle = stack->Particle(label);
833     if (!particle) {
834       
835       AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label));
836       return -1;
837     }
838     
839     // find mother
840     if (particle->GetMother(0) < 0) {
841
842       AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
843       return -1;
844     }
845     
846     label = particle->GetMother(0);
847   }
848   
849   return label;
850 }
851
852 //____________________________________________________________________
853 AliAODMCParticle* AliAnalysisTaskHighPtDeDxV0::FindPrimaryMotherAOD(AliAODMCParticle* startParticle, Int_t& nSteps)
854 {
855   //
856   // Finds the first mother among the primary particles of the particle identified by <label>,
857   // i.e. the primary that "caused" this particle
858   //
859   // Taken from AliPWG0Helper class
860   //
861
862   nSteps = 0;
863
864   AliAODMCParticle* mcPart = startParticle;
865
866   while (mcPart)
867     {
868       
869       if(mcPart->IsPrimary())
870         return mcPart;
871       
872       Int_t mother = mcPart->GetMother();
873
874       mcPart = dynamic_cast<AliAODMCParticle*>(fMCArray->At(mother));
875       nSteps++; // 1 level down
876     }
877
878   return 0;
879 }
880
881 //_____________________________________________________________________________
882 void AliAnalysisTaskHighPtDeDxV0::Sort(TClonesArray* array, Bool_t isMC) 
883 {
884   const Int_t n = array->GetEntriesFast(); 
885   if(n==0) {
886
887     if(isMC) 
888       fEvent->ptmaxMC = 0;
889     else
890       fEvent->ptmax   = 0;
891       
892     return;
893   }
894
895   Float_t ptArray[n];
896   Int_t   indexArray[n];
897   
898   for(Int_t i = 0; i < n; i++) {
899
900     Float_t pt = 0;
901     if(isMC) {
902       DeDxTrackMC* track = (DeDxTrackMC*)array->At(i);
903       pt = track->ptMC;
904     } else {
905       DeDxTrack* track = (DeDxTrack*)array->At(i);
906       pt = track->pt;
907     }
908     ptArray[i]    = pt;
909     indexArray[i] = i;
910   }
911
912   TMath::Sort(n, ptArray, indexArray);
913   
914   // set max pt
915   if(isMC) {
916     DeDxTrackMC* track = (DeDxTrackMC*)array->At(indexArray[0]);
917     fEvent->ptmaxMC = track->ptMC;
918   } else {
919     DeDxTrack* track = (DeDxTrack*)array->At(indexArray[0]);
920     fEvent->ptmax   = track->pt;
921   }
922   
923   // set order of each track
924   for(Int_t i = 0; i < n; i++) {
925     
926     if(isMC) {
927       DeDxTrackMC* track = (DeDxTrackMC*)array->At(indexArray[i]);
928       track->orderMC = i;
929     } else {
930       DeDxTrack* track = (DeDxTrack*)array->At(indexArray[i]);
931       track->order = i;
932     }
933   }
934 }
935 //__________________________________________________________________
936 void AliAnalysisTaskHighPtDeDxV0::ProduceArrayTrksESD( AliESDEvent *ESDevent, AnalysisMode analysisMode ){
937   Int_t nv0s = ESDevent->GetNumberOfV0s();
938   if(nv0s<1)return;
939   Int_t     trackmult = 0; // no pt cuts
940   Int_t     nadded    = 0;
941   const AliESDVertex *myBestPrimaryVertex = ESDevent->GetPrimaryVertex();
942   if (!myBestPrimaryVertex) return;
943   if (!(myBestPrimaryVertex->GetStatus())) return;
944   Double_t  lPrimaryVtxPosition[3];
945   myBestPrimaryVertex->GetXYZ(lPrimaryVtxPosition);
946   
947   Double_t  lPrimaryVtxCov[6];
948   myBestPrimaryVertex->GetCovMatrix(lPrimaryVtxCov);
949   Double_t  lPrimaryVtxChi2 = myBestPrimaryVertex->GetChi2toNDF();
950   
951   AliAODVertex* myPrimaryVertex = new AliAODVertex(lPrimaryVtxPosition, lPrimaryVtxCov, lPrimaryVtxChi2, NULL, -1, AliAODVertex::kPrimary);
952
953
954   if( analysisMode == kGlobalTrk ){
955     if(fV0ArrayGlobalPar)
956       fV0ArrayGlobalPar->Clear();
957   } else if( analysisMode == kTPCTrk ){
958     if(fV0ArrayTPCPar)
959       fV0ArrayTPCPar->Clear();
960   }
961   
962   if( analysisMode==kGlobalTrk ){  
963
964
965
966     
967     // ################################
968     // #### BEGINNING OF V0 CODE ######
969     // ################################
970
971     
972     for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
973       
974       // This is the begining of the V0 loop  
975       AliESDv0 *esdV0 = ESDevent->GetV0(iV0);
976       if (!esdV0) continue;
977       
978       // AliESDTrack (V0 Daughters)
979       UInt_t lKeyPos = (UInt_t)TMath::Abs(esdV0->GetPindex());
980       UInt_t lKeyNeg = (UInt_t)TMath::Abs(esdV0->GetNindex());
981       
982       AliESDtrack *pTrack = ESDevent->GetTrack(lKeyPos);
983       AliESDtrack *nTrack = ESDevent->GetTrack(lKeyNeg);
984       if (!pTrack || !nTrack) {
985         Printf("ERROR: Could not retreive one of the daughter track");
986         continue;
987       }
988       
989       // Remove like-sign
990       if (pTrack->GetSign() == nTrack->GetSign()) {
991         //cout<< "like sign, continue"<< endl;
992         continue;
993       } 
994       
995       // Eta cut on decay products
996       if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)
997         continue;
998       
999       // Pt cut on decay products
1000       if (esdV0->Pt() < fMinPt)
1001         //      if (pTrack->Pt() < fMinPt && nTrack->Pt() < fMinPt)
1002         continue;
1003
1004       //filter for positive track
1005       UShort_t filterFlag_p = 0;
1006       Bool_t filterCut_Set1_p = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
1007       Bool_t filterCut_Set2_p = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
1008       Bool_t filterCut_Set3_p = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
1009       
1010       UInt_t selectDebug_p = 0;
1011       if (fTrackFilterGolden) {
1012         selectDebug_p = fTrackFilterGolden->IsSelected(pTrack);
1013         if (selectDebug_p) {
1014           filterFlag_p +=1;
1015           filterCut_Set3_p=kTRUE;
1016         }
1017       }
1018       
1019       if (fTrackFilterTPC) {
1020         
1021         selectDebug_p = fTrackFilterTPC->IsSelected(pTrack);
1022         if (selectDebug_p){//only tracks which pass the TPC-only track cuts
1023           filterFlag_p +=2;
1024           filterCut_Set1_p=kTRUE;
1025           
1026         }
1027         
1028       }
1029       
1030       if (fTrackFilter) {
1031         selectDebug_p = fTrackFilter->IsSelected(pTrack);
1032         if (selectDebug_p) {
1033           filterCut_Set2_p=kTRUE;
1034         }
1035       }
1036       
1037       if(filterFlag_p ==0 )
1038         continue;
1039       
1040
1041       //filter for negative track
1042       UShort_t filterFlag_n = 0;
1043       Bool_t filterCut_Set1_n = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
1044       Bool_t filterCut_Set2_n = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
1045       Bool_t filterCut_Set3_n = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
1046       
1047       UInt_t selectDebug_n = 0;
1048       if (fTrackFilterGolden) {
1049         selectDebug_n = fTrackFilterGolden->IsSelected(nTrack);
1050         if (selectDebug_n) {
1051           filterFlag_n +=1;
1052           filterCut_Set3_n=kTRUE;
1053         }
1054       }
1055       
1056       if (fTrackFilterTPC) {
1057         
1058         selectDebug_n = fTrackFilterTPC->IsSelected(nTrack);
1059         if (selectDebug_n){//only tracks which pass the TPC-only track cuts
1060           filterFlag_n +=2;
1061           filterCut_Set1_n=kTRUE;
1062           
1063         }
1064         
1065       }
1066       
1067       if (fTrackFilter) {
1068         selectDebug_n = fTrackFilter->IsSelected(nTrack);
1069         if (selectDebug_n) {
1070           filterCut_Set2_n=kTRUE;
1071         }
1072       }
1073
1074
1075       // Check if switch does anything!
1076       Bool_t isSwitched = kFALSE;
1077       if (pTrack->GetSign() < 0) { // switch
1078         
1079         isSwitched = kTRUE;
1080         AliESDtrack* helpTrack = nTrack;
1081         nTrack = pTrack;
1082         pTrack = helpTrack;
1083       } 
1084       
1085       Float_t alpha = esdV0->AlphaV0();
1086       Float_t ptarm = esdV0->PtArmV0();
1087       // Double_t pVtxPos= v0->PrimaryVtxPosition();      
1088       
1089       Double_t  lV0Position[3];
1090       esdV0->GetXYZ(lV0Position[0], lV0Position[1], lV0Position[2]);
1091       
1092       Double_t lV0Radius      = TMath::Sqrt(lV0Position[0]*lV0Position[0]+lV0Position[1]*lV0Position[1]);
1093       Double_t lV0DecayLength = TMath::Sqrt(TMath::Power(lV0Position[0] - lPrimaryVtxPosition[0],2) +
1094                                             TMath::Power(lV0Position[1] - lPrimaryVtxPosition[1],2) +
1095                                             TMath::Power(lV0Position[2] - lPrimaryVtxPosition[2],2 ));
1096       AliKFVertex primaryVtxKF( *myPrimaryVertex );
1097       AliKFParticle::SetField(ESDevent->GetMagneticField());
1098       
1099       // Also implement switch here!!!!!!
1100       AliKFParticle* negEKF  = 0; // e-
1101       AliKFParticle* posEKF  = 0; // e+
1102       AliKFParticle* negPiKF = 0; // pi -
1103       AliKFParticle* posPiKF = 0; // pi +
1104       AliKFParticle* posPKF  = 0; // p
1105       AliKFParticle* negAPKF = 0; // p-bar
1106       
1107       if(!isSwitched) {
1108         negEKF  = new AliKFParticle( *(esdV0->GetParamN()) , 11);
1109         posEKF  = new AliKFParticle( *(esdV0->GetParamP()) ,-11);
1110         negPiKF = new AliKFParticle( *(esdV0->GetParamN()) ,-211);
1111         posPiKF = new AliKFParticle( *(esdV0->GetParamP()) , 211);
1112         posPKF  = new AliKFParticle( *(esdV0->GetParamP()) , 2212);
1113         negAPKF = new AliKFParticle( *(esdV0->GetParamN()) ,-2212);
1114       } else { // switch + and - 
1115         negEKF  = new AliKFParticle( *(esdV0->GetParamP()) , 11);
1116         posEKF  = new AliKFParticle( *(esdV0->GetParamN()) ,-11);
1117         negPiKF = new AliKFParticle( *(esdV0->GetParamP()) ,-211);
1118         posPiKF = new AliKFParticle( *(esdV0->GetParamN()) , 211);
1119         posPKF  = new AliKFParticle( *(esdV0->GetParamN()) , 2212);
1120         negAPKF = new AliKFParticle( *(esdV0->GetParamP()) ,-2212);
1121       }
1122       
1123       AliKFParticle v0GKF;  // Gamma e.g. from pi0
1124       v0GKF+=(*negEKF);
1125       v0GKF+=(*posEKF);
1126       v0GKF.SetProductionVertex(primaryVtxKF);
1127       
1128       AliKFParticle v0K0sKF; // K0 short
1129       v0K0sKF+=(*negPiKF);
1130       v0K0sKF+=(*posPiKF);
1131       v0K0sKF.SetProductionVertex(primaryVtxKF);
1132       
1133       AliKFParticle v0LambdaKF; // Lambda
1134       v0LambdaKF+=(*negPiKF);
1135       v0LambdaKF+=(*posPKF);    
1136       v0LambdaKF.SetProductionVertex(primaryVtxKF);
1137       
1138       AliKFParticle v0AntiLambdaKF; // Lambda-bar
1139       v0AntiLambdaKF+=(*posPiKF);
1140       v0AntiLambdaKF+=(*negAPKF);
1141       v0AntiLambdaKF.SetProductionVertex(primaryVtxKF);
1142       
1143       Double_t deltaInvMassG     = v0GKF.GetMass();
1144       Double_t deltaInvMassK0s   = v0K0sKF.GetMass()-0.498;
1145       Double_t deltaInvMassL     = v0LambdaKF.GetMass()-1.116;
1146       Double_t deltaInvMassAntiL = v0AntiLambdaKF.GetMass()-1.116;
1147       
1148       if(TMath::Abs(deltaInvMassK0s) > fMassCut &&
1149          TMath::Abs(deltaInvMassL) > fMassCut &&
1150          TMath::Abs(deltaInvMassAntiL) > fMassCut)
1151         continue;
1152       
1153       // Extract track information
1154       
1155       Short_t pcharge  = pTrack->Charge();
1156       Float_t ppt      = pTrack->Pt();
1157       Float_t pp       = pTrack->P(); 
1158       Float_t peta     = pTrack->Eta();
1159       Float_t pphi     = pTrack->Phi();
1160       Short_t pncl     = pTrack->GetTPCsignalN();
1161       Short_t pneff    = Short_t(pTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1162       Float_t pdedx    = pTrack->GetTPCsignal();
1163       
1164       Float_t ptpcchi  = 0;
1165       if(pTrack->GetTPCNcls() > 0)
1166         ptpcchi = pTrack->GetTPCchi2()/Float_t(pTrack->GetTPCNcls());
1167       
1168       Short_t ncharge  = nTrack->Charge();
1169       Float_t npt      = nTrack->Pt();
1170       Float_t np       = nTrack->P(); 
1171       Float_t neta     = nTrack->Eta();
1172       Float_t nphi     = nTrack->Phi();
1173       Short_t nncl     = nTrack->GetTPCsignalN();
1174       Short_t nneff    = Short_t(nTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1175       Float_t ndedx    = nTrack->GetTPCsignal();
1176       Float_t ntpcchi  = 0;
1177       if(nTrack->GetTPCNcls() > 0)
1178         ntpcchi = nTrack->GetTPCchi2()/Float_t(nTrack->GetTPCNcls());
1179       
1180       Float_t b[2];
1181       Float_t bCov[3];
1182       pTrack->GetImpactParameters(b,bCov);
1183       Float_t pdcaxy   = b[0];
1184       Float_t pdcaz    = b[1];
1185       nTrack->GetImpactParameters(b,bCov);
1186       Float_t ndcaxy   = b[0];
1187       Float_t ndcaz    = b[1];
1188       
1189       Int_t   primaryV0     = 0; // 0 means that the tracks are not both daughters of a primary particle (1 means they are)
1190       Int_t   pdgV0         = 0; // 0 means that they don't have same origin for MC (1 means they have the same original mother)
1191       Float_t p_ptMC        = 0;
1192       Short_t p_pidCode     = 0; // 0 = real data / no mc track!
1193       Short_t p_primaryFlag = 0; // 0 = real data / not primary mc track  
1194       Int_t   p_pdgMother   = 0;
1195       Float_t n_ptMC        = 0;
1196       Short_t n_pidCode     = 0; // 0 = real data / no mc track!
1197       Short_t n_primaryFlag = 0; // 0 = real data / not primary mc track  
1198       Int_t   n_pdgMother   = 0;
1199       if(fAnalysisMC) {
1200         
1201         Int_t p_mother_label = 0;
1202         Int_t p_mother_steps = 0;
1203         Int_t n_mother_label = 0;
1204         Int_t n_mother_steps = 0;
1205         
1206         // positive track
1207         const Int_t p_label = TMath::Abs(pTrack->GetLabel());
1208         TParticle* p_mcTrack = fMCStack->Particle(p_label);         
1209         if (p_mcTrack){
1210           
1211           if(fMCStack->IsPhysicalPrimary(p_label))
1212             p_primaryFlag = 1;
1213           
1214           Int_t p_pdgCode = p_mcTrack->GetPdgCode();
1215           p_pidCode = GetPidCode(p_pdgCode);
1216           
1217           p_ptMC      = p_mcTrack->Pt();
1218           
1219           p_mother_label = FindPrimaryMotherLabel(fMCStack, p_label, 
1220                                                   p_mother_steps);
1221           if(p_mother_label>0) {
1222             TParticle* p_mother = fMCStack->Particle(p_mother_label);
1223             p_pdgMother = p_mother->GetPdgCode();
1224           }
1225         }
1226         
1227         // negative track
1228         const Int_t n_label = TMath::Abs(nTrack->GetLabel());
1229         TParticle* n_mcTrack = fMCStack->Particle(n_label);         
1230         if (n_mcTrack){
1231           
1232           if(fMCStack->IsPhysicalPrimary(n_label))
1233             n_primaryFlag = 1;
1234           
1235           Int_t n_pdgCode = n_mcTrack->GetPdgCode();
1236           n_pidCode = GetPidCode(n_pdgCode);
1237           
1238           n_ptMC      = n_mcTrack->Pt();
1239           
1240           n_mother_label = FindPrimaryMotherLabel(fMCStack, n_label, 
1241                                                   n_mother_steps);
1242           if(n_mother_label>0) {
1243             TParticle* n_mother = fMCStack->Particle(n_mother_label);
1244             n_pdgMother = n_mother->GetPdgCode();
1245           }
1246         }
1247         
1248         // Check if V0 is primary = first and the same mother of both partciles
1249         if(p_mother_label>0 && n_mother_label>0 && p_mother_label == n_mother_label) {
1250           pdgV0 = p_pdgMother;
1251           if(p_mother_steps == 1 && n_mother_steps == 1) {
1252             primaryV0 = 1;
1253           }
1254         }
1255       }
1256
1257  
1258     
1259       if(fTreeOption) {
1260         
1261
1262
1263         DeDxV0* v0data = new((*fV0ArrayGlobalPar)[nadded]) DeDxV0();
1264         nadded++;
1265         
1266         // v0 data
1267         v0data->p       = esdV0->P();
1268         v0data->pt      = esdV0->Pt();
1269         v0data->eta     = esdV0->Eta();
1270         v0data->phi     = esdV0->Phi();
1271         v0data->pdca    = TMath::Sqrt(pdcaxy*pdcaxy + pdcaz*pdcaz);
1272         v0data->ndca    = TMath::Sqrt(ndcaxy*ndcaxy + ndcaz*ndcaz);
1273         v0data->dmassG  = deltaInvMassG;
1274         v0data->dmassK0 = deltaInvMassK0s;
1275         v0data->dmassL  = deltaInvMassL;
1276         v0data->dmassAL = deltaInvMassAntiL;
1277         v0data->alpha   = alpha;
1278         v0data->ptarm   = ptarm;
1279         v0data->decayr  = lV0Radius;
1280         v0data->decayl  = lV0DecayLength;
1281         
1282         // New parameters
1283         v0data->status  = esdV0->GetOnFlyStatus();
1284         v0data->chi2    = esdV0->GetChi2V0();
1285         v0data->cospt   = esdV0->GetV0CosineOfPointingAngle(); 
1286         // cospt: as I understand this means that the pointing to the vertex
1287         // is fine so I remove the dcaxy and dcaz for the V= class
1288         v0data->dcadaughters = esdV0->GetDcaV0Daughters();
1289         v0data->primary = primaryV0;
1290         v0data->pdg     = pdgV0;
1291         
1292         // positive track
1293         v0data->ptrack.p       = pp;
1294         v0data->ptrack.pt      = ppt;
1295         //        v0data->ptrack.ptcon   = ppt_con;
1296         //        v0data->ptrack.tpcchi  = ptpcchi;
1297         v0data->ptrack.eta     = peta;
1298         v0data->ptrack.phi     = pphi;
1299         v0data->ptrack.q       = pcharge;
1300         v0data->ptrack.ncl     = pncl;
1301         v0data->ptrack.neff    = pneff;
1302         v0data->ptrack.dedx    = pdedx;
1303         v0data->ptrack.dcaxy   = pdcaxy;
1304         v0data->ptrack.dcaz    = pdcaz;
1305         v0data->ptrack.pid     = p_pidCode;
1306         v0data->ptrack.primary = p_primaryFlag;
1307         v0data->ptrack.pttrue  = p_ptMC;
1308         v0data->ptrack.mother  = p_pdgMother;
1309         v0data->ptrack.filter  = filterFlag_p;
1310         v0data->ptrack.filterset1 = filterCut_Set1_p;
1311         v0data->ptrack.filterset2 = filterCut_Set2_p;
1312         v0data->ptrack.filterset3 = filterCut_Set3_p;
1313         
1314         // negative track
1315         v0data->ntrack.p       = np;
1316         v0data->ntrack.pt      = npt;
1317         //        v0data->ntrack.ptcon   = npt_con;
1318         //        v0data->ntrack.tpcchi  = ntpcchi;
1319         v0data->ntrack.eta     = neta;
1320         v0data->ntrack.phi     = nphi;
1321         v0data->ntrack.q       = ncharge;
1322         v0data->ntrack.ncl     = nncl;
1323         v0data->ntrack.neff    = nneff;
1324         v0data->ntrack.dedx    = ndedx;
1325         v0data->ntrack.dcaxy   = ndcaxy;
1326         v0data->ntrack.dcaz    = ndcaz;
1327         v0data->ntrack.pid     = n_pidCode;
1328         v0data->ntrack.primary = n_primaryFlag;
1329         v0data->ntrack.pttrue  = n_ptMC;
1330         v0data->ntrack.mother  = n_pdgMother;
1331         v0data->ntrack.filter  = filterFlag_n;
1332         v0data->ntrack.filterset1 = filterCut_Set1_n;
1333         v0data->ntrack.filterset2 = filterCut_Set2_n;
1334         v0data->ntrack.filterset3 = filterCut_Set3_n;
1335
1336
1337       }
1338       
1339       // clean up loop over v0
1340
1341       delete negPiKF;
1342       delete posPiKF;
1343       delete posPKF;
1344       delete negAPKF;
1345   
1346
1347
1348     }
1349   
1350     // clean up event
1351     //delete myPrimaryVertex;
1352   }
1353   else if( analysisMode==kTPCTrk ){  
1354     cout<<"&&&&&&&&&&&&&&&&&&&&&                Hello world"<<endl;
1355     const AliESDVertex *vtxSPD = ESDevent->GetPrimaryVertexSPD();
1356     if( vtxSPD->GetNContributors() < 1 || TMath::Abs(vtxSPD->GetZ()) > 10.0 ) return;
1357     
1358     
1359     // ################################
1360     // #### BEGINNING OF V0 CODE ######
1361     // ################################
1362
1363     
1364     for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
1365       
1366       // This is the begining of the V0 loop  
1367       AliESDv0 *esdV0 = ESDevent->GetV0(iV0);
1368       if (!esdV0) continue;
1369       
1370       // AliESDTrack (V0 Daughters)
1371       UInt_t lKeyPos = (UInt_t)TMath::Abs(esdV0->GetPindex());
1372       UInt_t lKeyNeg = (UInt_t)TMath::Abs(esdV0->GetNindex());
1373       
1374       AliESDtrack *pTrack = ESDevent->GetTrack(lKeyPos);
1375       AliESDtrack *nTrack = ESDevent->GetTrack(lKeyNeg);
1376       if (!pTrack || !nTrack) {
1377         Printf("ERROR: Could not retreive one of the daughter track");
1378         continue;
1379       }
1380
1381       // Remove like-sign
1382       if (pTrack->GetSign() == nTrack->GetSign()) {
1383         //cout<< "like sign, continue"<< endl;
1384         continue;
1385       } 
1386
1387       AliESDtrack *pTrackTPC = AliESDtrackCuts::GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(ESDevent),pTrack->GetID());
1388       AliESDtrack *nTrackTPC = AliESDtrackCuts::GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(ESDevent),nTrack->GetID());
1389
1390       if (!pTrackTPC || !nTrackTPC) {
1391         Printf("ERROR: Could not retreive one of the daughter TPC track");
1392         continue;
1393       }
1394
1395       //filter for positive track
1396       UInt_t selectDebug_p = 0;
1397       UShort_t filterFlag_p = 0;
1398       selectDebug_p = fTrackFilterTPC->IsSelected(pTrackTPC);
1399       //if(selectDebug_p==0) continue;
1400
1401       //filter for negative track
1402       UInt_t selectDebug_n = 0;
1403       UShort_t filterFlag_n = 0;
1404       selectDebug_n = fTrackFilterTPC->IsSelected(nTrackTPC);
1405       if(selectDebug_n==0 || selectDebug_p==0) continue;
1406
1407       if(selectDebug_p)
1408         filterFlag_p += 1;
1409
1410       if(pTrackTPC->Pt()>0.){
1411         // only constrain tracks above threshold
1412         AliExternalTrackParam exParamp;
1413         // take the B-field from the ESD, no 3D fieldMap available at this point
1414         Bool_t relate_p = false;
1415         relate_p = pTrackTPC->RelateToVertexTPC(vtxSPD,ESDevent->GetMagneticField(),
1416                                                 kVeryBig,&exParamp);
1417         if(!relate_p){
1418           delete pTrackTPC;
1419           continue;
1420         }
1421         pTrackTPC->Set(exParamp.GetX(),exParamp.GetAlpha(),exParamp.GetParameter(),
1422                        exParamp.GetCovariance());
1423       }
1424       else continue;     
1425       
1426
1427       //filter for negative track
1428  
1429       if(selectDebug_n)
1430         filterFlag_n += 1;
1431
1432   
1433       if(nTrackTPC->Pt()>0.){
1434         // only constrain tracks above threshold
1435         AliExternalTrackParam exParamn;
1436         // take the B-field from the ESD, no 3D fieldMap available at this point
1437         Bool_t relate_n = false;
1438         relate_n = nTrackTPC->RelateToVertexTPC(vtxSPD,ESDevent->GetMagneticField(),
1439                                                 kVeryBig,&exParamn);
1440         if(!relate_n){
1441           delete nTrackTPC;
1442           continue;
1443         }
1444         nTrackTPC->Set(exParamn.GetX(),exParamn.GetAlpha(),exParamn.GetParameter(),
1445                        exParamn.GetCovariance());
1446       }
1447       else continue;  
1448       
1449       
1450       // Eta cut on decay products
1451       if(TMath::Abs(pTrackTPC->Eta()) > fEtaCut || TMath::Abs(nTrackTPC->Eta()) > fEtaCut)
1452         continue;
1453       
1454       // Pt cut on decay products
1455       if (esdV0->Pt() < fMinPt)
1456         //      if (pTrack->Pt() < fMinPt && nTrack->Pt() < fMinPt)
1457         continue;
1458  
1459       
1460       // Check if switch does anything!
1461       Bool_t isSwitched = kFALSE;
1462       if (pTrackTPC->GetSign() < 0) { // switch
1463         
1464         isSwitched = kTRUE;
1465         AliESDtrack* helpTrack = nTrack;
1466         nTrackTPC = pTrackTPC;
1467         pTrackTPC = helpTrack;
1468       } 
1469       
1470
1471       // Extract track information
1472       
1473       Short_t pcharge  = pTrackTPC->Charge();
1474       Float_t ppt      = pTrackTPC->Pt();
1475       Float_t pp       = pTrackTPC->P(); 
1476       Float_t peta     = pTrackTPC->Eta();
1477       Float_t pphi     = pTrackTPC->Phi();
1478       Short_t pncl     = pTrackTPC->GetTPCsignalN();
1479       Short_t pneff    = Short_t(pTrackTPC->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1480       Float_t pdedx    = pTrackTPC->GetTPCsignal();
1481       
1482       Float_t ptpcchi  = 0;
1483       if(pTrackTPC->GetTPCNcls() > 0)
1484         ptpcchi = pTrackTPC->GetTPCchi2()/Float_t(pTrackTPC->GetTPCNcls());
1485       
1486       Short_t ncharge  = nTrackTPC->Charge();
1487       Float_t npt      = nTrackTPC->Pt();
1488       Float_t np       = nTrackTPC->P(); 
1489       Float_t neta     = nTrackTPC->Eta();
1490       Float_t nphi     = nTrackTPC->Phi();
1491       Short_t nncl     = nTrackTPC->GetTPCsignalN();
1492       Short_t nneff    = Short_t(nTrackTPC->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1493       Float_t ndedx    = nTrackTPC->GetTPCsignal();
1494       Float_t ntpcchi  = 0;
1495       if(nTrackTPC->GetTPCNcls() > 0)
1496         ntpcchi = nTrackTPC->GetTPCchi2()/Float_t(nTrackTPC->GetTPCNcls());
1497       
1498       Float_t bp[2]={0,0};
1499       Float_t bCovp[3]={0,0,0};
1500       pTrackTPC->GetImpactParameters(bp,bCovp);
1501       Float_t pdcaxy   = bp[0];
1502       Float_t pdcaz    = bp[1];
1503
1504       Float_t bn[2]={0,0};
1505       Float_t bCovn[3]={0,0,0};
1506       nTrackTPC->GetImpactParameters(bn,bCovn);
1507       Float_t ndcaxy   = bn[0];
1508       Float_t ndcaz    = bn[1];
1509
1510
1511       Float_t alpha = esdV0->AlphaV0();
1512       Float_t ptarm = esdV0->PtArmV0();
1513       // Double_t pVtxPos= v0->PrimaryVtxPosition();      
1514       
1515       Double_t  lV0Position[3];
1516       esdV0->GetXYZ(lV0Position[0], lV0Position[1], lV0Position[2]);
1517       
1518       Double_t lV0Radius      = TMath::Sqrt(lV0Position[0]*lV0Position[0]+lV0Position[1]*lV0Position[1]);
1519       Double_t lV0DecayLength = TMath::Sqrt(TMath::Power(lV0Position[0] - lPrimaryVtxPosition[0],2) +
1520                                             TMath::Power(lV0Position[1] - lPrimaryVtxPosition[1],2) +
1521                                             TMath::Power(lV0Position[2] - lPrimaryVtxPosition[2],2 ));
1522       AliKFVertex primaryVtxKF( *myPrimaryVertex );
1523       AliKFParticle::SetField(ESDevent->GetMagneticField());
1524       
1525       // Also implement switch here!!!!!!
1526       AliKFParticle* negEKF  = 0; // e-
1527       AliKFParticle* posEKF  = 0; // e+
1528       AliKFParticle* negPiKF = 0; // pi -
1529       AliKFParticle* posPiKF = 0; // pi +
1530       AliKFParticle* posPKF  = 0; // p
1531       AliKFParticle* negAPKF = 0; // p-bar
1532       
1533       
1534       if(!isSwitched) {
1535         /*      
1536                 negEKF  = new AliKFParticle( *(esdV0->GetParamN()) , 11);
1537                 posEKF  = new AliKFParticle( *(esdV0->GetParamP()) ,-11);
1538                 negPiKF = new AliKFParticle( *(esdV0->GetParamN()) ,-211);
1539                 posPiKF = new AliKFParticle( *(esdV0->GetParamP()) , 211);
1540                 posPKF  = new AliKFParticle( *(esdV0->GetParamP()) , 2212);
1541                 negAPKF = new AliKFParticle( *(esdV0->GetParamN()) ,-2212);
1542         */
1543
1544         negEKF  = new AliKFParticle( *(dynamic_cast<AliVTrack*>(nTrackTPC)) , 11);
1545         posEKF  = new AliKFParticle( *(dynamic_cast<AliVTrack*>(pTrackTPC)) ,-11);
1546         negPiKF = new AliKFParticle( *(dynamic_cast<AliVTrack*>(nTrackTPC)) ,-211);
1547         posPiKF = new AliKFParticle( *(dynamic_cast<AliVTrack*>(pTrackTPC)) , 211);
1548         posPKF  = new AliKFParticle( *(dynamic_cast<AliVTrack*>(pTrackTPC)) , 2212);
1549         negAPKF = new AliKFParticle( *(dynamic_cast<AliVTrack*>(nTrackTPC)) ,-2212);
1550
1551
1552
1553       } else { // switch + and - 
1554         /*
1555         negEKF  = new AliKFParticle( *(esdV0->GetParamP()) , 11);
1556         posEKF  = new AliKFParticle( *(esdV0->GetParamN()) ,-11);
1557         negPiKF = new AliKFParticle( *(esdV0->GetParamP()) ,-211);
1558         posPiKF = new AliKFParticle( *(esdV0->GetParamN()) , 211);
1559         posPKF  = new AliKFParticle( *(esdV0->GetParamN()) , 2212);
1560         negAPKF = new AliKFParticle( *(esdV0->GetParamP()) ,-2212);
1561         */
1562
1563         negEKF  = new AliKFParticle(  *(dynamic_cast<AliVTrack*>(pTrackTPC)), 11);
1564         posEKF  = new AliKFParticle(  *(dynamic_cast<AliVTrack*>(nTrackTPC)),-11);
1565         negPiKF = new AliKFParticle(  *(dynamic_cast<AliVTrack*>(pTrackTPC)),-211);
1566         posPiKF = new AliKFParticle(  *(dynamic_cast<AliVTrack*>(nTrackTPC)), 211);
1567         posPKF  = new AliKFParticle(  *(dynamic_cast<AliVTrack*>(nTrackTPC)), 2212);
1568         negAPKF = new AliKFParticle(  *(dynamic_cast<AliVTrack*>(pTrackTPC)),-2212);
1569
1570
1571       }
1572    
1573
1574
1575
1576       
1577       AliKFParticle v0GKF;  // Gamma e.g. from pi0
1578       v0GKF+=(*negEKF);
1579       v0GKF+=(*posEKF);
1580       v0GKF.SetProductionVertex(primaryVtxKF);
1581       
1582       AliKFParticle v0K0sKF; // K0 short
1583       v0K0sKF+=(*negPiKF);
1584       v0K0sKF+=(*posPiKF);
1585       v0K0sKF.SetProductionVertex(primaryVtxKF);
1586       
1587       AliKFParticle v0LambdaKF; // Lambda
1588       v0LambdaKF+=(*negPiKF);
1589       v0LambdaKF+=(*posPKF);    
1590       v0LambdaKF.SetProductionVertex(primaryVtxKF);
1591       
1592       AliKFParticle v0AntiLambdaKF; // Lambda-bar
1593       v0AntiLambdaKF+=(*posPiKF);
1594       v0AntiLambdaKF+=(*negAPKF);
1595       v0AntiLambdaKF.SetProductionVertex(primaryVtxKF);
1596       
1597       Double_t deltaInvMassG     = v0GKF.GetMass();
1598       Double_t deltaInvMassK0s   = v0K0sKF.GetMass()-0.498;
1599       Double_t deltaInvMassL     = v0LambdaKF.GetMass()-1.116;
1600       Double_t deltaInvMassAntiL = v0AntiLambdaKF.GetMass()-1.116;
1601       
1602       if(TMath::Abs(deltaInvMassK0s) > fMassCut &&
1603          TMath::Abs(deltaInvMassL) > fMassCut &&
1604          TMath::Abs(deltaInvMassAntiL) > fMassCut)
1605         continue;
1606       
1607       // TODO: Whe should these be different? Different mass hypothesis = energy loss
1608       // This is not important for us as we focus on the decay products!
1609       // Double_t ptK0s        = v0K0sKF.GetPt(); 
1610       // Double_t ptL          = v0LambdaKF.GetPt();
1611       // Double_t ptAntiL      = v0AntiLambdaKF.GetPt();     
1612       
1613       Int_t   primaryV0     = 0; // 0 means that the tracks are not both daughters of a primary particle (1 means they are)
1614       Int_t   pdgV0         = 0; // 0 means that they don't have same origin for MC (1 means they have the same original mother)
1615       Float_t p_ptMC        = 0;
1616       Short_t p_pidCode     = 0; // 0 = real data / no mc track!
1617       Short_t p_primaryFlag = 0; // 0 = real data / not primary mc track  
1618       Int_t   p_pdgMother   = 0;
1619       Float_t n_ptMC        = 0;
1620       Short_t n_pidCode     = 0; // 0 = real data / no mc track!
1621       Short_t n_primaryFlag = 0; // 0 = real data / not primary mc track  
1622       Int_t   n_pdgMother   = 0;
1623       if(fAnalysisMC) {
1624         
1625         Int_t p_mother_label = 0;
1626         Int_t p_mother_steps = 0;
1627         Int_t n_mother_label = 0;
1628         Int_t n_mother_steps = 0;
1629         
1630         // positive track
1631         const Int_t p_label = TMath::Abs(pTrackTPC->GetLabel());
1632         TParticle* p_mcTrack = fMCStack->Particle(p_label);         
1633         if (p_mcTrack){
1634           
1635           if(fMCStack->IsPhysicalPrimary(p_label))
1636             p_primaryFlag = 1;
1637           
1638           Int_t p_pdgCode = p_mcTrack->GetPdgCode();
1639           p_pidCode = GetPidCode(p_pdgCode);
1640           
1641           p_ptMC      = p_mcTrack->Pt();
1642           
1643           p_mother_label = FindPrimaryMotherLabel(fMCStack, p_label, 
1644                                                   p_mother_steps);
1645           if(p_mother_label>0) {
1646             TParticle* p_mother = fMCStack->Particle(p_mother_label);
1647             p_pdgMother = p_mother->GetPdgCode();
1648           }
1649         }
1650         
1651         // negative track
1652         const Int_t n_label = TMath::Abs(nTrackTPC->GetLabel());
1653         TParticle* n_mcTrack = fMCStack->Particle(n_label);         
1654         if (n_mcTrack){
1655           
1656           if(fMCStack->IsPhysicalPrimary(n_label))
1657             n_primaryFlag = 1;
1658           
1659           Int_t n_pdgCode = n_mcTrack->GetPdgCode();
1660           n_pidCode = GetPidCode(n_pdgCode);
1661           
1662           n_ptMC      = n_mcTrack->Pt();
1663           
1664           n_mother_label = FindPrimaryMotherLabel(fMCStack, n_label, 
1665                                                   n_mother_steps);
1666           if(n_mother_label>0) {
1667             TParticle* n_mother = fMCStack->Particle(n_mother_label);
1668             n_pdgMother = n_mother->GetPdgCode();
1669           }
1670         }
1671         
1672         // Check if V0 is primary = first and the same mother of both partciles
1673         if(p_mother_label>0 && n_mother_label>0 && p_mother_label == n_mother_label) {
1674           pdgV0 = p_pdgMother;
1675           if(p_mother_steps == 1 && n_mother_steps == 1) {
1676             primaryV0 = 1;
1677           }
1678         }
1679       }
1680       
1681
1682       if(fTreeOption) {
1683         
1684         DeDxV0* v0datatpc = new((*fV0ArrayTPCPar)[nadded]) DeDxV0();
1685         nadded++;
1686         
1687         // v0 data
1688         v0datatpc->p       = esdV0->P();
1689         v0datatpc->pt      = esdV0->Pt();
1690         v0datatpc->eta     = esdV0->Eta();
1691         v0datatpc->phi     = esdV0->Phi();
1692         v0datatpc->pdca    = TMath::Sqrt(pdcaxy*pdcaxy + pdcaz*pdcaz);
1693         v0datatpc->ndca    = TMath::Sqrt(ndcaxy*ndcaxy + ndcaz*ndcaz);
1694         v0datatpc->dmassG  = deltaInvMassG;
1695         v0datatpc->dmassK0 = deltaInvMassK0s;
1696         v0datatpc->dmassL  = deltaInvMassL;
1697         v0datatpc->dmassAL = deltaInvMassAntiL;
1698         v0datatpc->alpha   = alpha;
1699         v0datatpc->ptarm   = ptarm;
1700         v0datatpc->decayr  = lV0Radius;
1701         v0datatpc->decayl  = lV0DecayLength;
1702         
1703         // New parameters
1704         v0datatpc->status  = esdV0->GetOnFlyStatus();
1705         v0datatpc->chi2    = esdV0->GetChi2V0();
1706         v0datatpc->cospt   = esdV0->GetV0CosineOfPointingAngle(); 
1707         // cospt: as I understand this means that the pointing to the vertex
1708         // is fine so I remove the dcaxy and dcaz for the V= class
1709         v0datatpc->dcadaughters = esdV0->GetDcaV0Daughters();
1710         v0datatpc->primary = primaryV0;
1711         v0datatpc->pdg     = pdgV0;
1712         
1713         // positive track
1714         v0datatpc->ptrack.p       = pp;
1715         v0datatpc->ptrack.pt      = ppt;
1716         //        v0data->ptrack.ptcon   = ppt_con;
1717         //        v0data->ptrack.tpcchi  = ptpcchi;
1718         v0datatpc->ptrack.eta     = peta;
1719         v0datatpc->ptrack.phi     = pphi;
1720         v0datatpc->ptrack.q       = pcharge;
1721         v0datatpc->ptrack.ncl     = pncl;
1722         v0datatpc->ptrack.neff    = pneff;
1723         v0datatpc->ptrack.dedx    = pdedx;
1724         v0datatpc->ptrack.dcaxy   = pdcaxy;
1725         v0datatpc->ptrack.dcaz    = pdcaz;
1726         v0datatpc->ptrack.pid     = p_pidCode;
1727         v0datatpc->ptrack.primary = p_primaryFlag;
1728         v0datatpc->ptrack.pttrue  = p_ptMC;
1729         v0datatpc->ptrack.mother  = p_pdgMother;
1730         v0datatpc->ptrack.filter  = filterFlag_p;
1731         v0datatpc->ptrack.filterset1 = 0;
1732         v0datatpc->ptrack.filterset2 = 0;
1733         v0datatpc->ptrack.filterset3 = 0;
1734         
1735         // negative track
1736         v0datatpc->ntrack.p       = np;
1737         v0datatpc->ntrack.pt      = npt;
1738         //        v0data->ntrack.ptcon   = npt_con;
1739         //        v0data->ntrack.tpcchi  = ntpcchi;
1740         v0datatpc->ntrack.eta     = neta;
1741         v0datatpc->ntrack.phi     = nphi;
1742         v0datatpc->ntrack.q       = ncharge;
1743         v0datatpc->ntrack.ncl     = nncl;
1744         v0datatpc->ntrack.neff    = nneff;
1745         v0datatpc->ntrack.dedx    = ndedx;
1746         v0datatpc->ntrack.dcaxy   = ndcaxy;
1747         v0datatpc->ntrack.dcaz    = ndcaz;
1748         v0datatpc->ntrack.pid     = n_pidCode;
1749         v0datatpc->ntrack.primary = n_primaryFlag;
1750         v0datatpc->ntrack.pttrue  = n_ptMC;
1751         v0datatpc->ntrack.mother  = n_pdgMother;
1752         v0datatpc->ntrack.filter  = filterFlag_n;
1753         v0datatpc->ntrack.filterset1 = 0;
1754         v0datatpc->ntrack.filterset2 = 0;
1755         v0datatpc->ntrack.filterset3 = 0;
1756       }
1757       
1758       // clean up loop over v0
1759       
1760       delete negPiKF;
1761       delete posPiKF;
1762       delete posPKF;
1763       delete negAPKF;
1764  
1765
1766
1767     }
1768
1769  
1770  
1771
1772
1773   }
1774   delete myPrimaryVertex;
1775
1776   if(fTreeOption) {
1777     
1778     if( analysisMode==kGlobalTrk ){
1779
1780       
1781       fEvent->trackmult = trackmult;
1782       fEvent->n         = nadded;
1783
1784
1785     }
1786  
1787     
1788   }
1789
1790
1791   
1792 }
1793 //_______________________________________________________________________________________________________________
1794 void AliAnalysisTaskHighPtDeDxV0::ProduceArrayTrksAOD( AliAODEvent *AODevent, AnalysisMode analysisMode ){
1795   Int_t nv0s = AODevent->GetNumberOfV0s();
1796   if(nv0s<1)return;
1797   Int_t     trackmult = 0; // no pt cuts
1798   Int_t     nadded    = 0;
1799   
1800   AliAODVertex *myBestPrimaryVertex = AODevent->GetPrimaryVertex();
1801   if (!myBestPrimaryVertex) return;
1802   
1803   
1804   if( analysisMode == kGlobalTrk ){
1805     if(fV0ArrayGlobalPar)
1806       fV0ArrayGlobalPar->Clear();
1807   } else if( analysisMode == kTPCTrk ){
1808     if(fV0ArrayTPCPar)
1809       fV0ArrayTPCPar->Clear();
1810   }
1811   
1812   if( analysisMode==kGlobalTrk ){  
1813     
1814     
1815     // ################################
1816     // #### BEGINNING OF V0 CODE ######
1817     // ################################
1818     // This is the begining of the V0 loop  
1819     for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
1820       AliAODv0 *aodV0 = AODevent->GetV0(iV0);
1821       if (!aodV0) continue;
1822       
1823       // common part
1824       
1825       // AliAODTrack (V0 Daughters)
1826       AliAODVertex* vertex = aodV0->GetSecondaryVtx();
1827       if (!vertex) {
1828         Printf("ERROR: Could not retrieve vertex");
1829         continue;
1830       }
1831       
1832       AliAODTrack *pTrack = (AliAODTrack*)vertex->GetDaughter(0);
1833       AliAODTrack *nTrack = (AliAODTrack*)vertex->GetDaughter(1);
1834       if (!pTrack || !nTrack) {
1835         Printf("ERROR: Could not retrieve one of the daughter track");
1836         continue;
1837       }
1838       
1839       // Remove like-sign
1840       if (pTrack->Charge() == nTrack->Charge()) {
1841         //cout<< "like sign, continue"<< endl;
1842         continue;
1843       } 
1844       
1845       // Make sure charge ordering is ok
1846       if (pTrack->Charge() < 0) {
1847         AliAODTrack* helpTrack = pTrack;
1848         pTrack = nTrack;
1849         nTrack = helpTrack;
1850       } 
1851       
1852       // Eta cut on decay products
1853       if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)
1854         continue;
1855       
1856       // Pt cut on decay products
1857       if (aodV0->Pt() < fMinPt)
1858         //      if (pTrack->Pt() < fMinPt && nTrack->Pt() < fMinPt)
1859         continue;
1860       
1861       //check positive tracks
1862       UShort_t filterFlag_p = 0;
1863       Bool_t filterCut_Set1_p = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
1864       Bool_t filterCut_Set2_p = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
1865       Bool_t filterCut_Set3_p = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
1866       
1867       if (fTrackFilterGolden) {  
1868         // ITSTPC2010 cuts is bit 32 according to above macro, new versions of aliroot includes the golden cuts
1869         
1870         if(pTrack->TestFilterBit(32)) {
1871           filterFlag_p +=1;
1872           filterCut_Set3_p = kTRUE;
1873         }
1874       }
1875       
1876       
1877       if (fTrackFilterTPC) {
1878         // TPC only cuts is bit 1 according to above macro
1879         // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX, 
1880         if(pTrack->TestFilterBit(1)){
1881           filterFlag_p +=2;
1882           filterCut_Set1_p = kTRUE;
1883           
1884         }
1885       }
1886       
1887       if(filterFlag_p==0)
1888         continue;
1889       
1890       //check negative tracks
1891       UShort_t filterFlag_n = 0;
1892       Bool_t filterCut_Set1_n = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
1893       Bool_t filterCut_Set2_n = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
1894       Bool_t filterCut_Set3_n = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
1895       
1896       
1897       if (fTrackFilterGolden) {
1898         
1899         // ITSTPC2010 cuts is bit 32 according to above macro, new versions of aliroot includes the golden cuts
1900         if(nTrack->TestFilterBit(32)) {
1901           filterFlag_n +=1;
1902           filterCut_Set3_n = kTRUE;
1903         }
1904       }
1905       
1906       
1907       if (fTrackFilterTPC) {
1908         // TPC only cuts is bit 1 according to above macro
1909         // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX, 
1910         if(nTrack->TestFilterBit(1)){
1911           filterFlag_n +=2;
1912           filterCut_Set1_n = kTRUE;
1913           
1914         }
1915       }
1916       
1917       if(filterFlag_n==0)
1918         continue;
1919       
1920       
1921       Float_t alpha = aodV0->AlphaV0();
1922       Float_t ptarm = aodV0->PtArmV0();
1923       // Double_t pVtxPos= v0->PrimaryVtxPosition();      
1924       
1925       Double_t lV0Radius      = aodV0->RadiusV0();
1926       Double_t lV0DecayLength = aodV0->DecayLength(myBestPrimaryVertex);
1927       
1928       Double_t deltaInvMassG     = aodV0->InvMass2Prongs(0,1,11,11);
1929       Double_t deltaInvMassK0s   = aodV0->MassK0Short()-0.498;
1930       Double_t deltaInvMassL     = aodV0->MassLambda()-1.116;
1931       Double_t deltaInvMassAntiL = aodV0->MassAntiLambda()-1.116;
1932       
1933       if(TMath::Abs(deltaInvMassK0s) > fMassCut &&
1934          TMath::Abs(deltaInvMassL) > fMassCut &&
1935          TMath::Abs(deltaInvMassAntiL) > fMassCut)
1936         continue;
1937       
1938       // TODO: Why should these be different? Different mass hypothesis = energy loss
1939       // This is not important for us as we focus on the decay products!
1940       // Double_t ptK0s        = v0K0sKF.GetPt(); 
1941       // Double_t ptL          = v0LambdaKF.GetPt();
1942       // Double_t ptAntiL      = v0AntiLambdaKF.GetPt();     
1943       
1944       // Extract track information
1945       
1946       Double_t b[2], cov[3];
1947       if(!pTrack->PropagateToDCA(myBestPrimaryVertex, AODevent->GetMagneticField(), kVeryBig, b, cov))
1948         filterFlag_p += 32; // propagation failed!!!!!
1949       
1950       Float_t pdcaxy   = b[0];
1951       Float_t pdcaz    = b[1];
1952       if(!nTrack->PropagateToDCA(myBestPrimaryVertex, AODevent->GetMagneticField(), kVeryBig, b, cov))
1953         filterFlag_n += 32; // propagation failed!!!!!
1954       Float_t ndcaxy   = b[0];
1955       Float_t ndcaz    = b[1];
1956       
1957       Short_t pcharge  = pTrack->Charge();
1958       Float_t ppt      = pTrack->Pt();
1959       Float_t pp       = pTrack->P(); 
1960       Float_t peta     = pTrack->Eta();
1961       Float_t pphi     = pTrack->Phi();
1962       //        Float_t ptpcchi  = pTrack->Chi2perNDF();
1963       
1964       AliAODPid* pPid = pTrack->GetDetPid();
1965       Short_t pncl     = -10;
1966       Short_t pneff    = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1967       Float_t pdedx    = -10;
1968       Float_t pbeta = -99;
1969       if(pPid) {
1970         pncl     = pPid->GetTPCsignalN();
1971         pdedx    = pPid->GetTPCsignal();
1972         //TOF
1973         if (pTrack->GetStatus()&AliESDtrack::kTOFpid){
1974           Double_t tof[5];
1975           pPid->GetIntegratedTimes(tof);
1976           pbeta = tof[0]/pPid->GetTOFsignal();
1977         }
1978       }
1979       
1980       Short_t ncharge  = nTrack->Charge();
1981       Float_t npt      = nTrack->Pt();
1982       Float_t np       = nTrack->P(); 
1983       Float_t neta     = nTrack->Eta();
1984       Float_t nphi     = nTrack->Phi();
1985       //        Float_t ntpcchi  = nTrack->Chi2perNDF();
1986       
1987       AliAODPid* nPid = nTrack->GetDetPid();
1988       Short_t nncl     = -10;
1989       Short_t nneff    = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1990       Float_t ndedx    = -10;
1991       Float_t nbeta = -99;
1992       if(pPid) {
1993         nncl     = nPid->GetTPCsignalN();
1994         ndedx    = nPid->GetTPCsignal();
1995         //TOF
1996         if (nTrack->GetStatus()&AliESDtrack::kTOFpid){
1997           Double_t tof[5];
1998           nPid->GetIntegratedTimes(tof);
1999           nbeta = tof[0]/nPid->GetTOFsignal();
2000         }
2001       }
2002       
2003       Int_t   primaryV0     = 0; // 0 means that the tracks are not both daughters of a primary particle (1 means they are)
2004       Int_t   pdgV0         = 0; // 0 means that they don't have same origin for MC (1 means they have the same original mother)
2005       Float_t p_ptMC        = 0;
2006       Short_t p_pidCode     = 0; // 0 = real data / no mc track!
2007       Short_t p_primaryFlag = 0; // 0 = real data / not primary mc track  
2008       Int_t   p_pdgMother   = 0;
2009       Float_t n_ptMC        = 0;
2010       Short_t n_pidCode     = 0; // 0 = real data / no mc track!
2011       Short_t n_primaryFlag = 0; // 0 = real data / not primary mc track  
2012       Int_t   n_pdgMother   = 0;
2013       if(fAnalysisMC) {
2014         
2015         AliAODMCParticle* p_mother = 0;
2016         Int_t p_mother_steps = 0;
2017         AliAODMCParticle* n_mother = 0;
2018         Int_t n_mother_steps = 0;
2019         
2020         // positive track
2021         const Int_t p_label = TMath::Abs(pTrack->GetLabel());
2022         
2023         AliAODMCParticle* p_mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(p_label));
2024         if (p_mcTrack){
2025           
2026           if(p_mcTrack->IsPhysicalPrimary())
2027             p_primaryFlag = 1;
2028           
2029           Int_t p_pdgCode = p_mcTrack->GetPdgCode();
2030           p_pidCode = GetPidCode(p_pdgCode);
2031           
2032           p_ptMC      = p_mcTrack->Pt();
2033           
2034           p_mother = FindPrimaryMotherAOD(p_mcTrack, p_mother_steps);
2035           if(p_mother)
2036             p_pdgMother = p_mother->GetPdgCode();
2037         }
2038         
2039         // negative track
2040         const Int_t n_label = TMath::Abs(pTrack->GetLabel());
2041         
2042         AliAODMCParticle* n_mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(n_label));
2043         if (n_mcTrack){
2044           
2045           if(n_mcTrack->IsPhysicalPrimary())
2046             n_primaryFlag = 1;
2047           
2048           Int_t n_pdgCode = n_mcTrack->GetPdgCode();
2049           n_pidCode = GetPidCode(n_pdgCode);
2050           
2051           n_ptMC      = n_mcTrack->Pt();
2052           
2053           n_mother = FindPrimaryMotherAOD(n_mcTrack, n_mother_steps);
2054           if(n_mother)
2055             n_pdgMother = n_mother->GetPdgCode();
2056         }
2057         
2058         // Check if V0 is primary = first and the same mother of both partciles
2059         if(p_mother && n_mother && p_mother == n_mother) {
2060           pdgV0 = p_pdgMother;
2061           if(p_mother_steps == 1 && n_mother_steps == 1) {
2062             primaryV0 = 1;
2063           }
2064         }
2065       }
2066       
2067       if(fTreeOption) {
2068         
2069         DeDxV0* v0data = new((*fV0ArrayGlobalPar)[nadded]) DeDxV0();
2070         nadded++;
2071         
2072         // v0 data
2073         v0data->p       = aodV0->P();
2074         v0data->pt      = aodV0->Pt();
2075         v0data->eta     = aodV0->Eta();
2076         v0data->phi     = aodV0->Phi();
2077         v0data->pdca    = aodV0->DcaPosToPrimVertex();
2078         v0data->ndca    = aodV0->DcaNegToPrimVertex();
2079         v0data->dmassG  = deltaInvMassG;
2080         v0data->dmassK0 = deltaInvMassK0s;
2081         v0data->dmassL  = deltaInvMassL;
2082         v0data->dmassAL = deltaInvMassAntiL;
2083         v0data->alpha   = alpha;
2084         v0data->ptarm   = ptarm;
2085         v0data->decayr  = lV0Radius;
2086         v0data->decayl  = lV0DecayLength;
2087         // v0data->pdca    = TMath::Sqrt(pdcaxy*pdcaxy + pdcaz*pdcaz);
2088         // v0data->ndca    = TMath::Sqrt(ndcaxy*ndcaxy + ndcaz*ndcaz);
2089         
2090         // New parameters
2091         v0data->status  = aodV0->GetOnFlyStatus();
2092         v0data->chi2    = aodV0->Chi2V0();
2093         v0data->cospt   = aodV0->CosPointingAngle(myBestPrimaryVertex);
2094         // cospt: as I understand this means that the pointing to the vertex
2095         // is fine so I remove the dcaxy and dcaz for the V= class
2096         v0data->dcav0   = aodV0->DcaV0ToPrimVertex();
2097         v0data->dcadaughters = aodV0->DcaV0Daughters();
2098         v0data->primary = primaryV0;
2099         v0data->pdg     = pdgV0;
2100         
2101         // positive track
2102         v0data->ptrack.p       = pp;
2103         v0data->ptrack.pt      = ppt;
2104         //        v0data->ptrack.ptcon   = ppt_con;
2105         //        v0data->ptrack.tpcchi  = ptpcchi;
2106         v0data->ptrack.eta     = peta;
2107         v0data->ptrack.phi     = pphi;
2108         v0data->ptrack.q       = pcharge;
2109         v0data->ptrack.ncl     = pncl;
2110         v0data->ptrack.neff    = pneff;
2111         v0data->ptrack.dedx    = pdedx;
2112         v0data->ptrack.dcaxy   = pdcaxy;
2113         v0data->ptrack.dcaz    = pdcaz;
2114         v0data->ptrack.pid     = p_pidCode;
2115         v0data->ptrack.primary = p_primaryFlag;
2116         v0data->ptrack.pttrue  = p_ptMC;
2117         v0data->ptrack.mother  = p_pdgMother;
2118         v0data->ptrack.filter  = filterFlag_p;
2119         v0data->ptrack.filterset1 = filterCut_Set1_p;
2120         v0data->ptrack.filterset2 = filterCut_Set2_p;
2121         v0data->ptrack.filterset3 = filterCut_Set3_p;
2122         
2123         
2124         // negative track
2125         v0data->ntrack.p       = np;
2126         v0data->ntrack.pt      = npt;
2127         //        v0data->ntrack.ptcon   = npt_con;
2128         //        v0data->ntrack.tpcchi  = ntpcchi;
2129         v0data->ntrack.eta     = neta;
2130         v0data->ntrack.phi     = nphi;
2131         v0data->ntrack.q       = ncharge;
2132         v0data->ntrack.ncl     = nncl;
2133         v0data->ntrack.neff    = nneff;
2134         v0data->ntrack.dedx    = ndedx;
2135         v0data->ntrack.dcaxy   = ndcaxy;
2136         v0data->ntrack.dcaz    = ndcaz;
2137         v0data->ntrack.pid     = n_pidCode;
2138         v0data->ntrack.primary = n_primaryFlag;
2139         v0data->ntrack.pttrue  = n_ptMC;
2140         v0data->ntrack.mother  = n_pdgMother;
2141         v0data->ntrack.filter  = filterFlag_n;
2142         v0data->ntrack.filterset1 = filterCut_Set1_n;
2143         v0data->ntrack.filterset2 = filterCut_Set2_n;
2144         v0data->ntrack.filterset3 = filterCut_Set3_n;
2145         
2146         
2147       }
2148     }//end loop over v0's
2149     
2150     
2151   }else if( analysisMode==kTPCTrk ){  
2152     
2153     const AliAODVertex* vertexSPD= (AliAODVertex*)AODevent->GetPrimaryVertexSPD();//GetPrimaryVertex()
2154     cout<<"&&&&&&&&&&&&&&&&&&&&&                Hello world  0"<<endl;
2155     if( vertexSPD->GetNContributors() < 1 || TMath::Abs(vertexSPD->GetZ()) > 10.0 ) return;
2156     
2157     cout<<"&&&&&&&&&&&&&&&&&&&&&                Hello world"<<endl;
2158     // ################################
2159     // #### BEGINNING OF V0 CODE ######
2160     // ################################
2161     // This is the begining of the V0 loop  
2162
2163
2164
2165      for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
2166        AliAODv0 *aodV0 = AODevent->GetV0(iV0);
2167        if (!aodV0) continue;
2168        
2169        // common part
2170        
2171        // AliAODTrack (V0 Daughters)
2172        AliAODVertex* vertex = aodV0->GetSecondaryVtx();
2173        if (!vertex) {
2174          Printf("ERROR: Could not retrieve vertex");
2175          continue;
2176        }
2177        
2178        AliAODTrack *pTrack = (AliAODTrack*)vertex->GetDaughter(0);
2179        AliAODTrack *nTrack = (AliAODTrack*)vertex->GetDaughter(1);
2180        if (!pTrack || !nTrack) {
2181          Printf("ERROR: Could not retrieve one of the daughter track");
2182          continue;
2183        }
2184        
2185        
2186        // Remove like-sign
2187        if (pTrack->Charge() == nTrack->Charge()) {
2188          //cout<< "like sign, continue"<< endl;
2189          continue;
2190        } 
2191        
2192        // Make sure charge ordering is ok
2193        if (pTrack->Charge() < 0) {
2194          AliAODTrack* helpTrack = pTrack;
2195          pTrack = nTrack;
2196          nTrack = helpTrack;
2197        } 
2198        
2199        // Eta cut on decay products
2200        if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)
2201          continue;
2202
2203        cout<<"Eta positive track:"<<pTrack->Eta()<<endl;
2204
2205        
2206        // Pt cut on decay products
2207        if (aodV0->Pt() < fMinPt)
2208          //     if (pTrack->Pt() < fMinPt && nTrack->Pt() < fMinPt)
2209          continue;
2210        
2211        //check positive tracks
2212        UShort_t filterFlag_p = 0;
2213        Bool_t filterCut_Set1_p = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
2214        Bool_t filterCut_Set2_p = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
2215        Bool_t filterCut_Set3_p = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
2216        
2217        // TPC only cuts is bit 1 according to above macro
2218        // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX, 
2219        if(pTrack->TestFilterBit(128)) {
2220          cout<<"este track paso el corte bit 128"<<endl;
2221          filterFlag_p +=1;
2222        }
2223        cout<<"filterFlag_p="<<filterFlag_p<<endl;      
2224
2225
2226
2227
2228        if(filterFlag_p==0)
2229          continue;
2230  
2231
2232
2233        //check negative tracks
2234        UShort_t filterFlag_n = 0;
2235        Bool_t filterCut_Set1_n = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
2236        Bool_t filterCut_Set2_n = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
2237        Bool_t filterCut_Set3_n = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
2238        
2239        // TPC only cuts is bit 1 according to above macro
2240        // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX, 
2241        if(nTrack->TestFilterBit(128)) {
2242          filterFlag_n +=1;
2243        }
2244        
2245        if(filterFlag_n==0)
2246          continue;
2247        
2248        
2249        Float_t alpha = aodV0->AlphaV0();
2250        Float_t ptarm = aodV0->PtArmV0();
2251        // Double_t pVtxPos= v0->PrimaryVtxPosition();      
2252        
2253        Double_t lV0Radius      = aodV0->RadiusV0();
2254        Double_t lV0DecayLength = aodV0->DecayLength(myBestPrimaryVertex);
2255        
2256        Double_t deltaInvMassG     = aodV0->InvMass2Prongs(0,1,11,11);
2257        Double_t deltaInvMassK0s   = aodV0->MassK0Short()-0.498;
2258        Double_t deltaInvMassL     = aodV0->MassLambda()-1.116;
2259        Double_t deltaInvMassAntiL = aodV0->MassAntiLambda()-1.116;
2260        
2261        if(TMath::Abs(deltaInvMassK0s) > fMassCut &&
2262           TMath::Abs(deltaInvMassL) > fMassCut &&
2263           TMath::Abs(deltaInvMassAntiL) > fMassCut)
2264          continue;
2265        
2266        // TODO: Why should these be different? Different mass hypothesis = energy loss
2267        // This is not important for us as we focus on the decay products!
2268        // Double_t ptK0s        = v0K0sKF.GetPt(); 
2269        // Double_t ptL          = v0LambdaKF.GetPt();
2270        // Double_t ptAntiL      = v0AntiLambdaKF.GetPt();     
2271        
2272        // Extract track information
2273        
2274        Double_t b[2], cov[3];
2275        if(!pTrack->PropagateToDCA(vertexSPD, AODevent->GetMagneticField(), kVeryBig, b, cov))
2276          filterFlag_p += 32; // propagation failed!!!!!
2277        
2278        Float_t pdcaxy   = b[0];
2279        Float_t pdcaz    = b[1];
2280        if(!nTrack->PropagateToDCA(vertexSPD, AODevent->GetMagneticField(), kVeryBig, b, cov))
2281          filterFlag_n += 32; // propagation failed!!!!!
2282        Float_t ndcaxy   = b[0];
2283        Float_t ndcaz    = b[1];
2284        
2285        Short_t pcharge  = pTrack->Charge();
2286        Float_t ppt      = pTrack->Pt();
2287        Float_t pp       = pTrack->P(); 
2288        Float_t peta     = pTrack->Eta();
2289        Float_t pphi     = pTrack->Phi();
2290        //       Float_t ptpcchi  = pTrack->Chi2perNDF();
2291        
2292        AliAODPid* pPid = pTrack->GetDetPid();
2293        Short_t pncl     = -10;
2294        Short_t pneff    = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
2295        Float_t pdedx    = -10;
2296        Float_t pbeta = -99;
2297        if(pPid) {
2298          pncl     = pPid->GetTPCsignalN();
2299          pdedx    = pPid->GetTPCsignal();
2300          //TOF
2301          if (pTrack->GetStatus()&AliESDtrack::kTOFpid){
2302            Double_t tof[5];
2303            pPid->GetIntegratedTimes(tof);
2304            pbeta = tof[0]/pPid->GetTOFsignal();
2305          }
2306        }
2307        
2308        Short_t ncharge  = nTrack->Charge();
2309        Float_t npt      = nTrack->Pt();
2310        Float_t np       = nTrack->P(); 
2311        Float_t neta     = nTrack->Eta();
2312        Float_t nphi     = nTrack->Phi();
2313        //       Float_t ntpcchi  = nTrack->Chi2perNDF();
2314        
2315        AliAODPid* nPid = nTrack->GetDetPid();
2316        Short_t nncl     = -10;
2317        Short_t nneff    = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
2318        Float_t ndedx    = -10;
2319        Float_t nbeta = -99;
2320        if(pPid) {
2321          nncl     = nPid->GetTPCsignalN();
2322          ndedx    = nPid->GetTPCsignal();
2323          //TOF
2324          if (nTrack->GetStatus()&AliESDtrack::kTOFpid){
2325            Double_t tof[5];
2326            nPid->GetIntegratedTimes(tof);
2327            nbeta = tof[0]/nPid->GetTOFsignal();
2328          }
2329        }
2330        
2331        Int_t   primaryV0     = 0; // 0 means that the tracks are not both daughters of a primary particle (1 means they are)
2332        Int_t   pdgV0         = 0; // 0 means that they don't have same origin for MC (1 means they have the same original mother)
2333        Float_t p_ptMC        = 0;
2334        Short_t p_pidCode     = 0; // 0 = real data / no mc track!
2335        Short_t p_primaryFlag = 0; // 0 = real data / not primary mc track  
2336        Int_t   p_pdgMother   = 0;
2337        Float_t n_ptMC        = 0;
2338        Short_t n_pidCode     = 0; // 0 = real data / no mc track!
2339        Short_t n_primaryFlag = 0; // 0 = real data / not primary mc track  
2340        Int_t   n_pdgMother   = 0;
2341        if(fAnalysisMC) {
2342          
2343          AliAODMCParticle* p_mother = 0;
2344          Int_t p_mother_steps = 0;
2345          AliAODMCParticle* n_mother = 0;
2346          Int_t n_mother_steps = 0;
2347          
2348          // positive track
2349          const Int_t p_label = TMath::Abs(pTrack->GetLabel());
2350          
2351          AliAODMCParticle* p_mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(p_label));
2352          if (p_mcTrack){
2353            
2354            if(p_mcTrack->IsPhysicalPrimary())
2355              p_primaryFlag = 1;
2356            
2357            Int_t p_pdgCode = p_mcTrack->GetPdgCode();
2358            p_pidCode = GetPidCode(p_pdgCode);
2359            
2360            p_ptMC      = p_mcTrack->Pt();
2361            
2362            p_mother = FindPrimaryMotherAOD(p_mcTrack, p_mother_steps);
2363            if(p_mother)
2364              p_pdgMother = p_mother->GetPdgCode();
2365          }
2366          
2367          // negative track
2368          const Int_t n_label = TMath::Abs(pTrack->GetLabel());
2369          
2370          AliAODMCParticle* n_mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(n_label));
2371          if (n_mcTrack){
2372            
2373            if(n_mcTrack->IsPhysicalPrimary())
2374              n_primaryFlag = 1;
2375            
2376            Int_t n_pdgCode = n_mcTrack->GetPdgCode();
2377            n_pidCode = GetPidCode(n_pdgCode);
2378            
2379            n_ptMC      = n_mcTrack->Pt();
2380            
2381            n_mother = FindPrimaryMotherAOD(n_mcTrack, n_mother_steps);
2382            if(n_mother)
2383              n_pdgMother = n_mother->GetPdgCode();
2384          }
2385          
2386          // Check if V0 is primary = first and the same mother of both partciles
2387          if(p_mother && n_mother && p_mother == n_mother) {
2388            pdgV0 = p_pdgMother;
2389            if(p_mother_steps == 1 && n_mother_steps == 1) {
2390              primaryV0 = 1;
2391            }
2392          }
2393        }
2394        
2395        if(fTreeOption) {
2396          
2397          DeDxV0* v0datatpc = new((*fV0ArrayTPCPar)[nadded]) DeDxV0();
2398          nadded++;
2399          
2400          // v0 data
2401          v0datatpc->p       = aodV0->P();
2402          v0datatpc->pt      = aodV0->Pt();
2403          v0datatpc->eta     = aodV0->Eta();
2404          v0datatpc->phi     = aodV0->Phi();
2405          v0datatpc->pdca    = aodV0->DcaPosToPrimVertex();
2406          v0datatpc->ndca    = aodV0->DcaNegToPrimVertex();
2407          v0datatpc->dmassG  = deltaInvMassG;
2408          v0datatpc->dmassK0 = deltaInvMassK0s;
2409          v0datatpc->dmassL  = deltaInvMassL;
2410          v0datatpc->dmassAL = deltaInvMassAntiL;
2411          v0datatpc->alpha   = alpha;
2412          v0datatpc->ptarm   = ptarm;
2413          v0datatpc->decayr  = lV0Radius;
2414          v0datatpc->decayl  = lV0DecayLength;
2415          // v0data->pdca    = TMath::Sqrt(pdcaxy*pdcaxy + pdcaz*pdcaz);
2416          // v0data->ndca    = TMath::Sqrt(ndcaxy*ndcaxy + ndcaz*ndcaz);
2417          
2418          // New parameters
2419          v0datatpc->status  = aodV0->GetOnFlyStatus();
2420          v0datatpc->chi2    = aodV0->Chi2V0();
2421          v0datatpc->cospt   = aodV0->CosPointingAngle(myBestPrimaryVertex);
2422          // cospt: as I understand this means that the pointing to the vertex
2423          // is fine so I remove the dcaxy and dcaz for the V= class
2424          v0datatpc->dcav0   = aodV0->DcaV0ToPrimVertex();
2425          v0datatpc->dcadaughters = aodV0->DcaV0Daughters();
2426          v0datatpc->primary = primaryV0;
2427          v0datatpc->pdg     = pdgV0;
2428          
2429          // positive track
2430          v0datatpc->ptrack.p       = pp;
2431          v0datatpc->ptrack.pt      = ppt;
2432          //       v0data->ptrack.ptcon   = ppt_con;
2433          //       v0data->ptrack.tpcchi  = ptpcchi;
2434          v0datatpc->ptrack.eta     = peta;
2435          v0datatpc->ptrack.phi     = pphi;
2436          v0datatpc->ptrack.q       = pcharge;
2437          v0datatpc->ptrack.ncl     = pncl;
2438          v0datatpc->ptrack.neff    = pneff;
2439          v0datatpc->ptrack.dedx    = pdedx;
2440          v0datatpc->ptrack.dcaxy   = pdcaxy;
2441          v0datatpc->ptrack.dcaz    = pdcaz;
2442          v0datatpc->ptrack.pid     = p_pidCode;
2443          v0datatpc->ptrack.primary = p_primaryFlag;
2444          v0datatpc->ptrack.pttrue  = p_ptMC;
2445          v0datatpc->ptrack.mother  = p_pdgMother;
2446          v0datatpc->ptrack.filter  = filterFlag_p;
2447          v0datatpc->ptrack.filterset1 = 0;
2448          v0datatpc->ptrack.filterset2 = 0;
2449          v0datatpc->ptrack.filterset3 = 0;
2450          
2451          
2452          // negative track
2453          v0datatpc->ntrack.p       = np;
2454          v0datatpc->ntrack.pt      = npt;
2455          //       v0data->ntrack.ptcon   = npt_con;
2456          //       v0data->ntrack.tpcchi  = ntpcchi;
2457          v0datatpc->ntrack.eta     = neta;
2458          v0datatpc->ntrack.phi     = nphi;
2459          v0datatpc->ntrack.q       = ncharge;
2460          v0datatpc->ntrack.ncl     = nncl;
2461          v0datatpc->ntrack.neff    = nneff;
2462          v0datatpc->ntrack.dedx    = ndedx;
2463          v0datatpc->ntrack.dcaxy   = ndcaxy;
2464          v0datatpc->ntrack.dcaz    = ndcaz;
2465          v0datatpc->ntrack.pid     = n_pidCode;
2466          v0datatpc->ntrack.primary = n_primaryFlag;
2467          v0datatpc->ntrack.pttrue  = n_ptMC;
2468          v0datatpc->ntrack.mother  = n_pdgMother;
2469          v0datatpc->ntrack.filter  = filterFlag_n;
2470          v0datatpc->ntrack.filterset1 = 0;
2471          v0datatpc->ntrack.filterset2 = 0;
2472          v0datatpc->ntrack.filterset3 = 0;
2473          
2474          
2475        }
2476      }//end v0's loop
2477     
2478   }
2479   
2480   if(fTreeOption) {
2481     
2482     if( analysisMode==kGlobalTrk ){
2483
2484       
2485       fEvent->trackmult = trackmult;
2486       fEvent->n         = nadded;
2487       
2488       
2489     }
2490     
2491     
2492   }
2493   
2494   
2495   
2496 }
2497
2498
2499