Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / grid / AliAnalysisTaskHighPtDeDx.cxx
1 #include "AliAnalysisTaskHighPtDeDx.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 #include <AliESDtrackCuts.h>
20 #include <AliESDVZERO.h>
21 #include <AliAODVZERO.h>
22
23 #include <AliMCEventHandler.h>
24 #include <AliMCEvent.h>
25 #include <AliStack.h>
26
27 #include <TTreeStream.h>
28
29 #include <AliHeader.h>
30 #include <AliGenPythiaEventHeader.h>
31 #include <AliGenDPMjetEventHeader.h>
32
33 #include "AliCentrality.h" 
34
35 #include <AliAODTrack.h> 
36 #include <AliAODPid.h> 
37 #include <AliAODMCHeader.h> 
38
39
40 // STL includes
41 #include <iostream>
42 using namespace std;
43
44
45 //
46 // Responsible:
47 // Alexandru Dobrin (Wayne State) 
48 // Peter Christiansen (Lund)
49 //
50
51 /* 
52 To do:
53
54 Separate the code into two
55
56 */
57
58
59
60
61 ClassImp(AliAnalysisTaskHighPtDeDx)
62
63 const Double_t AliAnalysisTaskHighPtDeDx::fgkClight = 2.99792458e-2;
64
65 //_____________________________________________________________________________
66 AliAnalysisTaskHighPtDeDx::AliAnalysisTaskHighPtDeDx():
67   AliAnalysisTaskSE(),
68   fESD(0x0),
69   fAOD(0x0),
70   fMC(0x0),
71   fMCStack(0x0),
72   fMCArray(0x0),
73   fTrackFilter(0x0),
74   fTrackFilterGolden(0x0),
75   fTrackFilterTPC(0x0),
76   fAnalysisType("ESD"),
77   fAnalysisMC(kFALSE),
78   fAnalysisPbPb(kFALSE),
79   fTPCBranch(kFALSE),
80   ftrigBit1(0x0),
81   ftrigBit2(0x0),
82   fRandom(0x0),
83   fEvent(0x0),
84   fTrackArrayGlobalPar(0x0),
85   fTrackArrayTPCPar(0x0),
86   fTrackArrayMC(0x0),
87   fVZEROArray(0x0),
88   fVtxCut(10.0),  
89   fEtaCut(0.9),  
90   fMinPt(0.1),
91   fMinCent(0.0),
92   fMaxCent(100.0),
93   fLowPtFraction(0.01),
94   fTreeOption(0),
95   fMcProcessType(-999),
96   fTriggeredEventMB(-999),
97   fVtxStatus(-999),
98   fZvtx(-999),
99   fZvtxMC(-999),
100   fRun(-999),
101   fEventId(-999),
102   fListOfObjects(0), 
103   fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),
104   fTree(0x0),
105   fn1(0),
106   fn2(0),
107   fcent(0)
108
109
110 {
111   // Default constructor (should not be used)
112
113   //  fRandom = new TRandom(0); // 0 means random seed
114 }
115
116 //______________________________________________________________________________
117 AliAnalysisTaskHighPtDeDx::AliAnalysisTaskHighPtDeDx(const char *name):
118   AliAnalysisTaskSE(name),
119   fESD(0x0),
120   fAOD(0x0),
121   fMC(0x0),
122   fMCStack(0x0),
123   fMCArray(0x0),
124   fTrackFilter(0x0),
125   fTrackFilterGolden(0x0),
126   fTrackFilterTPC(0x0),
127   fAnalysisType("ESD"),
128   fAnalysisMC(kFALSE),
129   fAnalysisPbPb(kFALSE),
130   fTPCBranch(kFALSE),
131   ftrigBit1(0x0),
132   ftrigBit2(0x0),
133   fRandom(0x0),
134   fEvent(0x0),
135   fTrackArrayGlobalPar(0x0),
136   fTrackArrayMC(0x0),
137   fVtxCut(10.0),  
138   fEtaCut(0.9),  
139   fMinPt(0.1),
140   fLowPtFraction(0.01),
141   fTreeOption(0),
142   fMcProcessType(-999),
143   fTriggeredEventMB(-999),
144   fVtxStatus(-999),
145   fZvtx(-999),
146   fZvtxMC(-999),
147   fRun(-999),
148   fEventId(-999),
149   fListOfObjects(0), 
150   fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),
151   fTree(0x0),
152   fn1(0),
153   fn2(0),
154   fcent(0)
155
156 {
157
158   if(fTree)fTree=0;
159   // Output slot #1 writes into a TList
160   DefineOutput(1, TList::Class());
161 }
162
163 //_____________________________________________________________________________
164 AliAnalysisTaskHighPtDeDx::~AliAnalysisTaskHighPtDeDx()
165 {
166   // Destructor
167   // histograms are in the output list and deleted when the output
168   // list is deleted by the TSelector dtor
169   if (fListOfObjects && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
170     delete fListOfObjects;
171     fListOfObjects = 0;
172   }
173   if (fRandom) delete fRandom;
174   fRandom=0;
175   
176   // //for proof running; I cannot create tree do to memory limitations -> work with THnSparse 
177   // if (fListOfObjects  && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutputList;
178   
179   
180   
181 }
182
183 //______________________________________________________________________________
184 void AliAnalysisTaskHighPtDeDx::UserCreateOutputObjects()
185
186   // This method is called once per worker node
187   // Here we define the output: histograms and debug tree if requested 
188   // We also create the random generator here so it might get different seeds...
189   fRandom = new TRandom(0); // 0 means random seed
190
191   OpenFile(1);
192   fListOfObjects = new TList();
193   fListOfObjects->SetOwner();
194   
195   //
196   // Histograms
197   //  
198   fEvents = new TH1I("fEvents","Number of analyzed events; Events; Counts", 3, 0, 3);
199   fListOfObjects->Add(fEvents);
200
201   fn1=new TH1F("fn1","fn1",5001,-1,5000);
202   fListOfObjects->Add(fn1);
203
204   fn2=new TH1F("fn2","fn2",5001,-1,5000);
205   fListOfObjects->Add(fn2);
206
207   fcent=new TH1F("fcent","fcent",104,-2,102);
208   fListOfObjects->Add(fcent);
209
210   fVtx = new TH1I("fVtx","Vtx info (0=no, 1=yes); Vtx; Counts", 2, -0.5, 1.5);
211   fListOfObjects->Add(fVtx);
212
213   fVtxBeforeCuts = new TH1F("fVtxBeforeCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
214   fListOfObjects->Add(fVtxBeforeCuts);
215   
216   fVtxAfterCuts = new TH1F("fVtxAfterCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
217   fListOfObjects->Add(fVtxAfterCuts);
218
219   if (fAnalysisMC) {    
220     fVtxMC = new TH1I("fVtxMC","Vtx info - no trigger cut (0=no, 1=yes); Vtx; Counts", 2, -0.5, 1.5);
221     fListOfObjects->Add(fVtxMC);
222   }
223
224   if (fTreeOption) {
225
226     fTree = new TTree("tree","Event data");
227     fEvent = new DeDxEvent();
228     fTree->Branch("event", &fEvent);
229
230     fTrackArrayGlobalPar = new TClonesArray("DeDxTrack", 1000);
231     fTree->Bronch("trackGlobalPar", "TClonesArray", &fTrackArrayGlobalPar);
232     if(fTPCBranch){
233       fTrackArrayTPCPar = new TClonesArray("DeDxTrack", 1000);
234       fTree->Bronch("trackTPCPar", "TClonesArray", &fTrackArrayTPCPar);
235     }
236
237     fVZEROArray = new TClonesArray("VZEROCell", 1000);
238     fTree->Bronch("cellVZERO", "TClonesArray", &fVZEROArray);
239
240     if (fAnalysisMC) {    
241       fTrackArrayMC = new TClonesArray("DeDxTrackMC", 1000);
242       fTree->Bronch("trackMC", "TClonesArray", &fTrackArrayMC);
243     }
244
245     fTree->SetDirectory(0);
246
247     fListOfObjects->Add(fTree);
248
249   }
250
251   // Post output data.
252   PostData(1, fListOfObjects);
253 }
254
255 //______________________________________________________________________________
256 void AliAnalysisTaskHighPtDeDx::UserExec(Option_t *) 
257 {
258   // Main loop
259
260   //
261   // First we make sure that we have valid input(s)!
262   //
263   AliVEvent *event = InputEvent();
264   if (!event) {
265      Error("UserExec", "Could not retrieve event");
266      return;
267   }
268
269
270   if (fAnalysisType == "ESD"){
271     fESD = dynamic_cast<AliESDEvent*>(event);
272     if(!fESD){
273       Printf("%s:%d ESDEvent not found in Input Manager",(char*)__FILE__,__LINE__);
274       this->Dump();
275       return;
276     }    
277   } else {
278     fAOD = dynamic_cast<AliAODEvent*>(event);
279     if(!fAOD){
280       Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
281       this->Dump();
282       return;
283     }    
284   }
285
286
287
288   if (fAnalysisMC) {
289
290     if (fAnalysisType == "ESD"){
291       fMC = dynamic_cast<AliMCEvent*>(MCEvent());
292       if(!fMC){
293         Printf("%s:%d MCEvent not found in Input Manager",(char*)__FILE__,__LINE__);
294         this->Dump();
295         return;
296       }    
297
298       fMCStack = fMC->Stack();
299       
300       if(!fMCStack){
301         Printf("%s:%d MCStack not found in Input Manager",(char*)__FILE__,__LINE__);
302         this->Dump();
303         return;
304       }    
305     } else { // AOD
306
307       fMC = dynamic_cast<AliMCEvent*>(MCEvent());
308       if(fMC)
309         fMC->Dump();
310
311       fMCArray = (TClonesArray*)fAOD->FindListObject("mcparticles");
312       if(!fMCArray){
313         Printf("%s:%d AOD MC array not found in Input Manager",(char*)__FILE__,__LINE__);
314         this->Dump();
315         return;
316       }    
317     }
318   }
319
320   
321   // Get trigger decision
322   fTriggeredEventMB = 0; //init
323  
324
325   if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
326      ->IsEventSelected() & ftrigBit1 ){
327     fn1->Fill(1);
328     fTriggeredEventMB = 1;  //event triggered as minimum bias
329   }
330   if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
331      ->IsEventSelected() & ftrigBit2 ){
332     // From AliVEvent:
333     //    kINT7         = BIT(1), // V0AND trigger, offline V0 selection
334     fTriggeredEventMB += 2;  
335     fn2->Fill(1);
336   }
337
338
339   // Get process type for MC
340   fMcProcessType = 0; // -1=invalid, 0=data, 1=ND, 2=SD, 3=DD
341
342   // real data that are not triggered we skip
343   if(!fAnalysisMC && !fTriggeredEventMB)
344     return; 
345   
346   if (fAnalysisMC) {
347     
348     if (fAnalysisType == "ESD"){
349
350       AliHeader* headerMC = fMC->Header();
351       if (headerMC) {
352         
353         AliGenEventHeader* genHeader = headerMC->GenEventHeader();
354         TArrayF vtxMC(3); // primary vertex  MC 
355         vtxMC[0]=9999; vtxMC[1]=9999;  vtxMC[2]=9999; //initialize with dummy
356         if (genHeader) {
357           genHeader->PrimaryVertex(vtxMC);
358         }
359         fZvtxMC = vtxMC[2];
360         
361         // PYTHIA:
362         AliGenPythiaEventHeader* pythiaGenHeader =
363           dynamic_cast<AliGenPythiaEventHeader*>(headerMC->GenEventHeader());
364         if (pythiaGenHeader) {  //works only for pythia
365           fMcProcessType =  GetPythiaEventProcessType(pythiaGenHeader->ProcessType());
366         }
367         // PHOJET:
368         AliGenDPMjetEventHeader* dpmJetGenHeader =
369           dynamic_cast<AliGenDPMjetEventHeader*>(headerMC->GenEventHeader());
370         if (dpmJetGenHeader) {
371           fMcProcessType = GetDPMjetEventProcessType(dpmJetGenHeader->ProcessType());
372         }
373       }
374     } else { // AOD
375       
376       AliAODMCHeader* mcHeader = dynamic_cast<AliAODMCHeader*>(fAOD->FindListObject("mcHeader")); 
377       if(mcHeader) {
378         fZvtxMC = mcHeader->GetVtxZ();
379         
380         if(strstr(mcHeader->GetGeneratorName(), "Pythia")) {
381           fMcProcessType =  GetPythiaEventProcessType(mcHeader->GetEventType());
382         } else {
383           fMcProcessType =  GetDPMjetEventProcessType(mcHeader->GetEventType());
384         }
385       }
386     }
387   }
388   
389   // There are 3 cases
390   // Vertex: NO  - status -1
391   // Vertex: YES : outside cut - status 0
392   //             : inside cut  - status 1
393   // We have to be careful how we normalize because we probably want to
394   // normalize to:
395   // Nevents=(No vertex + outside + inside)/(out + in)*in
396   
397   
398   if (fAnalysisType == "ESD")
399     fZvtx = GetVertex(fESD);
400   else // AOD
401     fZvtx = GetVertex(fAOD);
402     
403   fVtxStatus = -999;
404   
405   if(fZvtx<-990) {
406     
407     fVtxStatus = -1;
408     if(fTriggeredEventMB)
409       fVtx->Fill(0);
410     if(fAnalysisMC)
411       fVtxMC->Fill(0);
412   } else {
413     
414     if(fTriggeredEventMB)
415       fVtx->Fill(1);
416     if(fAnalysisMC)
417       fVtxMC->Fill(1);
418     fVtxBeforeCuts->Fill(fZvtx);
419     fVtxStatus = 0;
420     if (TMath::Abs(fZvtx) < fVtxCut) {  
421       fVtxAfterCuts->Fill(fZvtx);
422       fVtxStatus = 1;
423     }
424   }
425   
426   // store MC event data no matter what
427   if(fAnalysisMC) {
428
429     if (fAnalysisType == "ESD") {
430       ProcessMCTruthESD();
431     } else { // AOD
432       ProcessMCTruthAOD();
433     }
434   }      
435   
436
437
438   Float_t centrality = -10;
439   // only analyze triggered events
440   if(fTriggeredEventMB) {
441     
442     if (fAnalysisType == "ESD"){
443       if(fAnalysisPbPb){
444         AliCentrality *centObject = fESD->GetCentrality();
445         centrality = centObject->GetCentralityPercentile("V0M"); 
446         if((centrality>fMaxCent)||(centrality<fMinCent))return;
447       }
448       fcent->Fill(centrality);
449       AnalyzeESD(fESD);
450     } else { // AOD
451       if(fAnalysisPbPb){
452         AliCentrality *centObject = fAOD->GetCentrality();
453         if(centObject)
454           centrality = centObject->GetCentralityPercentile("V0M"); 
455         if((centrality>fMaxCent)||(centrality<fMinCent))return;
456       }
457       fcent->Fill(centrality);
458       AnalyzeAOD(fAOD);
459     }
460   }
461
462   if( fTreeOption) {
463     
464     fEvent->process = fMcProcessType;
465     fEvent->trig    = fTriggeredEventMB;
466     fEvent->zvtxMC  = fZvtxMC;
467     fEvent->cent      = centrality;
468
469     fTree->Fill();
470     fTrackArrayGlobalPar->Clear();
471     if(fTPCBranch)
472       fTrackArrayTPCPar->Clear();
473     fVZEROArray->Clear();
474
475     if (fAnalysisMC)    
476       fTrackArrayMC->Clear();
477   }
478   
479   // Post output data.
480   PostData(1, fListOfObjects);
481 }
482
483 //________________________________________________________________________
484 void AliAnalysisTaskHighPtDeDx::AnalyzeESD(AliESDEvent* esdEvent)
485 {
486   fRun  = esdEvent->GetRunNumber();
487   fEventId = 0;
488   if(esdEvent->GetHeader())
489     fEventId = GetEventIdAsLong(esdEvent->GetHeader());
490   
491   Short_t isPileup = esdEvent->IsPileupFromSPD();
492   
493   //  Int_t     event     = esdEvent->GetEventNumberInFile();
494   UInt_t    time      = esdEvent->GetTimeStamp();
495   //  ULong64_t trigger   = esdEvent->GetTriggerMask();
496   Float_t   magf      = esdEvent->GetMagneticField();
497
498
499
500
501
502   if(fTriggeredEventMB) {// Only MC case can we have not triggered events
503     
504     // accepted event
505     fEvents->Fill(0);
506     
507     
508     if(fVtxStatus!=1) return; // accepted vertex
509     Int_t nESDTracks = esdEvent->GetNumberOfTracks();
510     
511     if(nESDTracks<1)return;
512     cout<<"Multiplicity="<<nESDTracks<<endl;
513     
514     ProduceArrayTrksESD( esdEvent, kGlobalTrk );//produce array with global track parameters
515     if(fTPCBranch)
516       ProduceArrayTrksESD( esdEvent, kTPCTrk );//array with tpc track parametes
517
518     fEvents->Fill(1);
519     AliESDVZERO *esdV0 = esdEvent->GetVZEROData();// loop sobre canales del V0 para obtener las multiplicidad
520
521     for (Int_t iCh=0; iCh<64; ++iCh) { 
522       Float_t multv=esdV0->GetMultiplicity(iCh); 
523       Int_t intexv=iCh;
524       VZEROCell* cellv0 = new((*fVZEROArray)[iCh]) VZEROCell();
525       cellv0->cellindex=intexv;
526       cellv0->cellmult= multv;
527     }   
528
529
530
531   } // end if triggered
532   
533   if(fTreeOption) {
534
535     fEvent->run       = fRun;
536     fEvent->eventid   = fEventId;
537     fEvent->time      = time;
538     //fEvent->cent      = centrality;
539     fEvent->mag       = magf;
540     fEvent->zvtx      = fZvtx;
541     fEvent->vtxstatus = fVtxStatus;
542     fEvent->pileup    = isPileup;
543
544   }
545
546
547
548
549 }
550
551 //________________________________________________________________________
552 void AliAnalysisTaskHighPtDeDx::AnalyzeAOD(AliAODEvent* aodEvent)
553 {
554   fRun  = aodEvent->GetRunNumber();
555   fEventId = 0;
556   if(aodEvent->GetHeader())
557     fEventId = GetEventIdAsLong(aodEvent->GetHeader());
558    
559   UInt_t    time      = 0; // Missing AOD info? aodEvent->GetTimeStamp();
560   Float_t   magf      = aodEvent->GetMagneticField();
561
562   //Int_t     trackmult = 0; // no pt cuts
563   //Int_t     nadded    = 0;
564
565   Short_t isPileup = aodEvent->IsPileupFromSPD();
566
567
568
569
570   if(fTriggeredEventMB) {// Only MC case can we have not triggered events
571     
572     // accepted event
573     fEvents->Fill(0);
574     
575     if(fVtxStatus!=1) return; // accepted vertex
576     Int_t nAODTracks = aodEvent->GetNumberOfTracks();
577     if(nAODTracks<1)return;      
578
579     ProduceArrayTrksAOD( aodEvent, kGlobalTrk );
580     if(fTPCBranch)
581       ProduceArrayTrksAOD( aodEvent, kTPCTrk );
582     fEvents->Fill(1);
583
584     AliAODVZERO *esdV0 = aodEvent->GetVZEROData();// loop sobre canales del V0 para obtener las multiplicidad
585
586     for (Int_t iCh=0; iCh<64; ++iCh) { 
587       Float_t multv=esdV0->GetMultiplicity(iCh); 
588       Int_t intexv=iCh;
589       VZEROCell* cellv0 = new((*fVZEROArray)[iCh]) VZEROCell();
590       cellv0->cellindex=intexv;
591       cellv0->cellmult= multv;
592     }   
593
594
595
596
597   } // end if triggered
598   
599   if(fTreeOption) {
600
601     //Sort(fTrackArrayGlobalPar, kFALSE);
602
603     fEvent->run       = fRun;
604     fEvent->eventid   = fEventId;
605     fEvent->time      = time;
606     //fEvent->cent      = centrality;
607     fEvent->mag       = magf;
608     fEvent->zvtx      = fZvtx;
609     fEvent->vtxstatus = fVtxStatus;
610     //fEvent->trackmult = trackmult;
611     //fEvent->n         = nadded;
612     fEvent->pileup    = isPileup;
613   }
614 }
615
616 //_____________________________________________________________________________
617 Float_t AliAnalysisTaskHighPtDeDx::GetVertex(const AliVEvent* event) const
618 {
619   Float_t zvtx = -999;
620   
621   const AliVVertex* primaryVertex = event->GetPrimaryVertex(); 
622   
623   if(primaryVertex->GetNContributors()>0)
624     zvtx = primaryVertex->GetZ();
625
626   return zvtx;
627 }
628
629 //_____________________________________________________________________________
630 Short_t AliAnalysisTaskHighPtDeDx::GetPidCode(Int_t pdgCode) const 
631 {
632   // return our internal code for pions, kaons, and protons
633   
634   Short_t pidCode = 6;
635   
636   switch (TMath::Abs(pdgCode)) {
637   case 211:
638     pidCode = 1; // pion
639     break;
640   case 321:
641     pidCode = 2; // kaon
642     break;
643   case 2212:
644     pidCode = 3; // proton
645     break;
646   case 11:
647     pidCode = 4; // electron
648     break;
649   case 13:
650     pidCode = 5; // muon
651     break;
652   default:
653     pidCode = 6;  // something else?
654   };
655   
656   return pidCode;
657 }
658
659
660 //_____________________________________________________________________________
661 void AliAnalysisTaskHighPtDeDx::ProcessMCTruthESD() 
662 {
663   // Fill the special MC histogram with the MC truth info
664   
665   Short_t trackmult = 0;
666   Short_t nadded    = 0;
667   const Int_t nTracksMC = fMCStack->GetNtrack();
668   
669   for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
670     
671     //Cuts
672     if(!(fMCStack->IsPhysicalPrimary(iTracks)))
673       continue;
674     
675     TParticle* trackMC = fMCStack->Particle(iTracks);
676     
677     Double_t chargeMC = trackMC->GetPDG()->Charge();
678     if (chargeMC == 0)
679       continue;
680     
681     if (TMath::Abs(trackMC->Eta()) > fEtaCut )
682       continue;
683     
684     trackmult++;
685     
686     //   "charge:pt:p:eta:phi:pidCode"
687     Float_t ptMC      = trackMC->Pt();
688     Float_t pMC       = trackMC->P();
689     Float_t etaMC     = trackMC->Eta();
690     Float_t phiMC     = trackMC->Phi();
691     
692     Int_t pdgCode = trackMC->GetPdgCode();
693     Short_t pidCodeMC = 0;
694     pidCodeMC = GetPidCode(pdgCode);
695     
696     // Here we want to add some of the MC histograms!
697     
698     // And therefore we first cut here!
699     if (trackMC->Pt() < fMinPt) {
700       
701       // Keep small fraction of low pT tracks
702       if(fRandom->Rndm() > fLowPtFraction)
703         continue;
704     } // else {
705     // Here we want to add the high pt part of the MC histograms!
706     //    }
707     
708     if(fTreeOption) {
709       
710       DeDxTrackMC* track = new((*fTrackArrayMC)[nadded]) DeDxTrackMC();
711       nadded++;
712       
713       track->pMC   = pMC;
714       track->ptMC  = ptMC;
715       track->etaMC = etaMC;
716       track->phiMC = phiMC;
717       track->qMC   = Short_t(chargeMC);
718       track->pidMC = pidCodeMC;
719       track->pdgMC = pdgCode;
720     }
721     
722   }//MC track loop
723   
724   if(fTreeOption) {
725     
726     Sort(fTrackArrayMC, kTRUE);
727
728     fEvent->trackmultMC = trackmult;
729     fEvent->nMC         = nadded;
730   }
731   
732 }
733
734 //_____________________________________________________________________________
735 void AliAnalysisTaskHighPtDeDx::ProcessMCTruthAOD() 
736 {
737   // Fill the special MC histogram with the MC truth info
738
739   Short_t trackmult = 0;
740   Short_t nadded    = 0;
741   const Int_t nTracksMC = fMCArray->GetEntriesFast();
742   
743   for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
744     
745     AliAODMCParticle* trackMC = dynamic_cast<AliAODMCParticle*>(fMCArray->At(iTracks));
746     
747     //Cuts
748     if(!(trackMC->IsPhysicalPrimary()))
749       continue;
750     
751     
752     Double_t chargeMC = trackMC->Charge();
753     if (chargeMC == 0)
754       continue;
755     
756     if (TMath::Abs(trackMC->Eta()) > fEtaCut )
757       continue;
758     
759     trackmult++;
760     
761     //   "charge:pt:p:eta:phi:pidCode"
762     Float_t ptMC      = trackMC->Pt();
763     Float_t pMC       = trackMC->P();
764     Float_t etaMC     = trackMC->Eta();
765     Float_t phiMC     = trackMC->Phi();
766     
767     Int_t pdgCode = trackMC->PdgCode();
768     Short_t pidCodeMC = 0;
769     pidCodeMC = GetPidCode(pdgCode);
770     
771     // Here we want to add some of the MC histograms!
772     
773     // And therefore we first cut here!
774     if (trackMC->Pt() < fMinPt) {
775       
776       // Keep small fraction of low pT tracks
777       if(fRandom->Rndm() > fLowPtFraction)
778         continue;
779     } // else {
780     // Here we want to add the high pt part of the MC histograms!
781     //    }
782     
783     if(fTreeOption) {
784       
785       DeDxTrackMC* track = new((*fTrackArrayMC)[nadded]) DeDxTrackMC();
786       nadded++;
787       
788       track->pMC   = pMC;
789       track->ptMC  = ptMC;
790       track->etaMC = etaMC;
791       track->phiMC = phiMC;
792       track->qMC   = Short_t(chargeMC);
793       track->pidMC = pidCodeMC;
794       track->pdgMC = pdgCode;
795     }
796     
797   }//MC track loop
798   
799   if(fTreeOption) {
800     
801     Sort(fTrackArrayMC, kTRUE);
802
803     fEvent->trackmultMC = trackmult;
804     fEvent->nMC         = nadded;
805   }
806   
807 }
808
809 //_____________________________________________________________________________
810 Short_t AliAnalysisTaskHighPtDeDx::GetPythiaEventProcessType(Int_t pythiaType) {
811   //
812   // Get the process type of the event.  PYTHIA
813   //
814   // source PWG0   dNdpt 
815
816   Short_t globalType = -1; //init
817       
818   if(pythiaType==92||pythiaType==93){
819     globalType = 2; //single diffractive
820   }
821   else if(pythiaType==94){
822     globalType = 3; //double diffractive
823   }
824   //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB??
825   else {
826     globalType = 1;  //non diffractive
827   }
828   return globalType;
829 }
830
831 //_____________________________________________________________________________
832 Short_t AliAnalysisTaskHighPtDeDx::GetDPMjetEventProcessType(Int_t dpmJetType) {
833   //
834   // get the process type of the event.  PHOJET
835   //
836   //source PWG0   dNdpt 
837   // can only read pythia headers, either directly or from cocktalil header
838   Short_t globalType = -1;
839   
840   if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction
841     globalType = 1;
842   }
843   else if (dpmJetType==5 || dpmJetType==6) {
844     globalType = 2;
845   }
846   else if (dpmJetType==7) {
847     globalType = 3;
848   }
849   return globalType;
850 }
851
852 //_____________________________________________________________________________
853 ULong64_t AliAnalysisTaskHighPtDeDx::GetEventIdAsLong(AliVHeader* header) const
854 {
855   // To have a unique id for each event in a run!
856   // Modified from AliRawReader.h
857   return ((ULong64_t)header->GetBunchCrossNumber()+
858           (ULong64_t)header->GetOrbitNumber()*3564+
859           (ULong64_t)header->GetPeriodNumber()*16777215*3564);
860 }
861
862
863 //____________________________________________________________________
864 TParticle* AliAnalysisTaskHighPtDeDx::FindPrimaryMother(AliStack* stack, Int_t label)
865 {
866   //
867   // Finds the first mother among the primary particles of the particle identified by <label>,
868   // i.e. the primary that "caused" this particle
869   //
870   // Taken from AliPWG0Helper class
871   //
872
873   Int_t motherLabel = FindPrimaryMotherLabel(stack, label);
874   if (motherLabel < 0)
875     return 0;
876
877   return stack->Particle(motherLabel);
878 }
879
880 //____________________________________________________________________
881 Int_t AliAnalysisTaskHighPtDeDx::FindPrimaryMotherLabel(AliStack* stack, Int_t label)
882 {
883   //
884   // Finds the first mother among the primary particles of the particle identified by <label>,
885   // i.e. the primary that "caused" this particle
886   //
887   // returns its label
888   //
889   // Taken from AliPWG0Helper class
890   //
891   const Int_t nPrim  = stack->GetNprimary();
892   
893   while (label >= nPrim) {
894
895     //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
896
897     TParticle* particle = stack->Particle(label);
898     if (!particle) {
899       
900       AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label));
901       return -1;
902     }
903     
904     // find mother
905     if (particle->GetMother(0) < 0) {
906
907       AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
908       return -1;
909     }
910     
911     label = particle->GetMother(0);
912   }
913   
914   return label;
915 }
916
917 //____________________________________________________________________
918 AliAODMCParticle* AliAnalysisTaskHighPtDeDx::FindPrimaryMotherAOD(AliAODMCParticle* startParticle)
919 {
920   //
921   // Finds the first mother among the primary particles of the particle identified by <label>,
922   // i.e. the primary that "caused" this particle
923   //
924   // Taken from AliPWG0Helper class
925   //
926
927   AliAODMCParticle* mcPart = startParticle;
928
929   while (mcPart)
930     {
931       
932       if(mcPart->IsPrimary())
933         return mcPart;
934       
935       Int_t mother = mcPart->GetMother();
936
937       mcPart = dynamic_cast<AliAODMCParticle*>(fMCArray->At(mother));
938     }
939
940   return 0;
941 }
942  
943 //_____________________________________________________________________________
944 void AliAnalysisTaskHighPtDeDx::Sort(TClonesArray* array, Bool_t isMC) 
945 {
946   const Int_t n = array->GetEntriesFast(); 
947   if(n==0) {
948
949     if(isMC) 
950       fEvent->ptmaxMC = 0;
951     else
952       fEvent->ptmax   = 0;
953       
954     return;
955   }
956
957   Float_t ptArray[n];
958   Int_t   indexArray[n];
959   
960   for(Int_t i = 0; i < n; i++) {
961
962     Float_t pt = 0;
963     if(isMC) {
964       DeDxTrackMC* track = (DeDxTrackMC*)array->At(i);
965       pt = track->ptMC;
966     } else {
967       DeDxTrack* track = (DeDxTrack*)array->At(i);
968       pt = track->pt;
969     }
970     ptArray[i]    = pt;
971     indexArray[i] = i;
972   }
973
974   TMath::Sort(n, ptArray, indexArray);
975   
976   // set max pt
977   if(isMC) {
978     DeDxTrackMC* track = (DeDxTrackMC*)array->At(indexArray[0]);
979     fEvent->ptmaxMC = track->ptMC;
980   } else {
981     DeDxTrack* track = (DeDxTrack*)array->At(indexArray[0]);
982     fEvent->ptmax   = track->pt;
983   }
984   
985   // set order of each track
986   for(Int_t i = 0; i < n; i++) {
987     
988     if(isMC) {
989       DeDxTrackMC* track = (DeDxTrackMC*)array->At(indexArray[i]);
990       track->orderMC = i;
991     } else {
992       DeDxTrack* track = (DeDxTrack*)array->At(indexArray[i]);
993       track->order = i;
994     }
995   }
996 }
997 //__________________________________________________________________
998 void AliAnalysisTaskHighPtDeDx::ProduceArrayTrksESD( AliESDEvent *ESDevent, AnalysisMode analysisMode ){
999   
1000   const Int_t nESDTracks = ESDevent->GetNumberOfTracks();
1001   Int_t trackmult=0;
1002   Int_t nadded=0;
1003   if( analysisMode == kGlobalTrk ){
1004     if(fTrackArrayGlobalPar)
1005       fTrackArrayGlobalPar->Clear();
1006   } else if( analysisMode == kTPCTrk ){
1007     if(fTrackArrayTPCPar)
1008       fTrackArrayTPCPar->Clear();
1009   }
1010   
1011   if( analysisMode==kGlobalTrk ){  
1012     
1013     for(Int_t iT = 0; iT < nESDTracks; iT++) {
1014       
1015       AliESDtrack* esdTrack = ESDevent->GetTrack(iT);
1016       
1017       
1018       if(TMath::Abs(esdTrack->Eta()) > fEtaCut)
1019         continue;
1020       
1021       UShort_t filterFlag = 0;
1022       Bool_t filterCut_Set1 = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
1023       Bool_t filterCut_Set2 = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
1024       Bool_t filterCut_Set3 = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
1025       
1026       UInt_t selectDebug = 0;
1027       if (fTrackFilterGolden) {
1028         selectDebug = fTrackFilterGolden->IsSelected(esdTrack);
1029         if (selectDebug) {
1030           filterFlag +=1;
1031           filterCut_Set3=kTRUE;
1032         }
1033       }
1034       
1035       if (fTrackFilterTPC) {
1036         
1037         selectDebug = fTrackFilterTPC->IsSelected(esdTrack);
1038         if (selectDebug){//only tracks which pass the TPC-only track cuts
1039           trackmult++;
1040           filterFlag +=2;
1041           filterCut_Set1=kTRUE;
1042           
1043         }
1044         
1045       }
1046       
1047       if (fTrackFilter) {
1048         selectDebug = fTrackFilter->IsSelected(esdTrack);
1049         if (selectDebug) {
1050           filterCut_Set2=kTRUE;
1051         }
1052       }
1053       
1054      
1055       if(filterFlag==0)
1056         continue;
1057       
1058       
1059       //
1060       // Track was accepted
1061       //      
1062       
1063       // Here we want to add histograms!
1064       
1065       if (esdTrack->Pt() < fMinPt) {
1066         
1067         // Keep small fraction of low pT tracks
1068         if(fRandom->Rndm() > fLowPtFraction)
1069           continue;
1070       } // else {
1071       // Here we want to add the high pt part of the histograms!
1072       //    }
1073     
1074       Short_t charge  = esdTrack->Charge();
1075       Float_t pt      = esdTrack->Pt();
1076       Float_t p       = esdTrack->P(); 
1077       Float_t eta     = esdTrack->Eta();
1078       Float_t phi     = esdTrack->Phi();
1079       Short_t ncl     = esdTrack->GetTPCsignalN();
1080       Short_t neff    = Short_t(esdTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1081       //          Short_t nclf    = esdTrack->GetTPCNclsF();
1082       Float_t dedx    = esdTrack->GetTPCsignal();
1083       Float_t tpcchi  = 0;
1084       if(esdTrack->GetTPCNcls() > 0)
1085         tpcchi = esdTrack->GetTPCchi2()/Float_t(esdTrack->GetTPCNcls());
1086       Float_t b[2];
1087       Float_t bCov[3];
1088       esdTrack->GetImpactParameters(b,bCov);
1089       Float_t dcaxy   = b[0];
1090       Float_t dcaz    = b[1];
1091       Double_t p_con[3] = {0, 0, 0};
1092       esdTrack->GetConstrainedPxPyPz(p_con);
1093       
1094       
1095       //        Float_t pt_con = (Float_t)TMath::Sqrt(p_con[0]*p_con[0] + p_con[1]*p_con[1]);
1096       // const AliExternalTrackParam* tpcParam = esdTrack->GetTPCInnerParam();
1097       // Float_t pttpc   = tpcParam->Pt();
1098       // Float_t ptpc    = tpcParam->P();
1099       
1100       Float_t ptMC        = 0;
1101       Short_t pidCode     = 0; // 0 = real data / no mc track!
1102       Short_t primaryFlag = 0; // 0 = real data / not primary mc track  
1103       Int_t   pdgMother   = 0;
1104       
1105       
1106       //with Globaltrack parameters????
1107       if(fAnalysisMC) {
1108         
1109         const Int_t label = TMath::Abs(esdTrack->GetLabel());
1110         TParticle* mcTrack = fMCStack->Particle(label);     
1111         if (mcTrack){
1112           
1113           if(fMCStack->IsPhysicalPrimary(label))
1114             primaryFlag = 1;
1115           
1116           Int_t pdgCode = mcTrack->GetPdgCode();
1117           pidCode = GetPidCode(pdgCode);
1118           
1119           ptMC      = mcTrack->Pt();
1120           
1121           TParticle* mother = FindPrimaryMother(fMCStack, label);
1122           pdgMother = mother->GetPdgCode();
1123         }
1124       }
1125       
1126     
1127       //TOF
1128       Float_t beta = -99;
1129       if (esdTrack->GetStatus()&AliESDtrack::kTOFpid){
1130         if ((esdTrack->GetIntegratedLength() != 0) && (esdTrack->GetTOFsignal() != 0))
1131           beta = esdTrack->GetIntegratedLength()/esdTrack->GetTOFsignal()/fgkClight;
1132       }
1133       
1134       if(fTreeOption) {
1135         
1136         DeDxTrack* track = new((*fTrackArrayGlobalPar)[nadded]) DeDxTrack();
1137         nadded++;
1138         
1139         track->p          = p;
1140         track->pt         = pt;
1141         //        track->ptcon   = pt_con;
1142         //        track->tpcchi  = tpcchi;
1143         track->eta        = eta;
1144         track->phi        = phi;
1145         track->q          = charge;
1146         track->filter     = filterFlag;
1147         track->ncl        = ncl;
1148         track->neff       = neff;
1149         track->dedx       = dedx;
1150         track->beta       = beta;
1151         track->dcaxy      = dcaxy;
1152         track->dcaz       = dcaz;
1153         track->pid        = pidCode;
1154         track->primary    = primaryFlag;
1155         track->pttrue     = ptMC;
1156         track->mother     = pdgMother;
1157         track->filterset1 = filterCut_Set1;
1158         track->filterset2 = filterCut_Set2;
1159         track->filterset3 = filterCut_Set3;
1160         
1161         
1162       }
1163     }//end of track loop
1164
1165   }//end global
1166   
1167   else if( analysisMode==kTPCTrk ){  
1168     const AliESDVertex *vtxSPD = ESDevent->GetPrimaryVertexSPD();
1169     if( vtxSPD->GetNContributors() < 1 || TMath::Abs(vtxSPD->GetZ()) > 10.0 ) return;
1170  
1171
1172     for(Int_t iT = 0; iT < nESDTracks; iT++) {
1173       
1174       AliESDtrack* esdTrack = ESDevent->GetTrack(iT);
1175
1176       AliESDtrack *trackTPC = AliESDtrackCuts::GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(ESDevent),esdTrack->GetID());
1177       
1178       if(!trackTPC) continue;
1179     
1180       UInt_t selectDebug = 0;
1181       selectDebug = fTrackFilterTPC->IsSelected(trackTPC);
1182       if(selectDebug==0) continue;
1183     
1184       if(trackTPC->Pt()>0.){
1185         // only constrain tracks above threshold
1186         AliExternalTrackParam exParam;
1187         // take the B-field from the ESD, no 3D fieldMap available at this point
1188         Bool_t relate = false;
1189         relate = trackTPC->RelateToVertexTPC(vtxSPD,ESDevent->GetMagneticField(),
1190                                            kVeryBig,&exParam);
1191         if(!relate){
1192           delete trackTPC;
1193           continue;
1194         }
1195         trackTPC->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),
1196                       exParam.GetCovariance());
1197       }
1198       else continue;     
1199       
1200  
1201
1202
1203       if(TMath::Abs(trackTPC->Eta()) > fEtaCut)
1204         continue;
1205       
1206       //
1207       // Track was accepted
1208       //      
1209       
1210       // Here we want to add histograms!
1211       
1212       if (trackTPC->Pt() < fMinPt) {
1213         
1214         // Keep small fraction of low pT tracks
1215         if(fRandom->Rndm() > fLowPtFraction)
1216           continue;
1217       } // else {
1218       // Here we want to add the high pt part of the histograms!
1219       //    }
1220       
1221       Short_t charge  = trackTPC->Charge();
1222       Float_t pt      = trackTPC->Pt();
1223       Float_t p       = trackTPC->P(); 
1224       Float_t eta     = trackTPC->Eta();
1225       Float_t phi     = trackTPC->Phi();
1226       Short_t ncl     = trackTPC->GetTPCsignalN();
1227       Short_t neff    = Short_t(trackTPC->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1228       //          Short_t nclf    = esdTrack->GetTPCNclsF();
1229       Float_t dedx    = trackTPC->GetTPCsignal();
1230       Float_t tpcchi  = 0;
1231       if(trackTPC->GetTPCNcls() > 0)
1232         tpcchi = trackTPC->GetTPCchi2()/Float_t(trackTPC->GetTPCNcls());
1233       //Float_t b[2];
1234       //Float_t bCov[3];
1235       //trackTPC->GetImpactParameters(b,bCov);
1236       //Float_t dcaxy   = b[0];
1237       //Float_t dcaz    = b[1];
1238
1239       Float_t dcaxy = 0.;
1240       Float_t dcaz = 0.;
1241       trackTPC->GetImpactParameters(dcaxy,dcaz);
1242
1243
1244       Double_t p_con[3] = {0, 0, 0};
1245       trackTPC->GetConstrainedPxPyPz(p_con);
1246       
1247     
1248       // Float_t pt_con = (Float_t)TMath::Sqrt(p_con[0]*p_con[0] + p_con[1]*p_con[1]);
1249       // const AliExternalTrackParam* tpcParam = esdTrack->GetTPCInnerParam();
1250       // Float_t pttpc   = tpcParam->Pt();
1251       // Float_t ptpc    = tpcParam->P();
1252       
1253       Float_t ptMC        = 0;
1254       Short_t pidCode     = 0; // 0 = real data / no mc track!
1255       Short_t primaryFlag = 0; // 0 = real data / not primary mc track  
1256       Int_t   pdgMother   = 0;
1257       
1258       
1259       //with Globaltrack parameters????
1260       if(fAnalysisMC) {
1261         
1262         const Int_t label = TMath::Abs(esdTrack->GetLabel());
1263         TParticle* mcTrack = fMCStack->Particle(label);     
1264         if (mcTrack){
1265           
1266           if(fMCStack->IsPhysicalPrimary(label))
1267             primaryFlag = 1;
1268           
1269           Int_t pdgCode = mcTrack->GetPdgCode();
1270           pidCode = GetPidCode(pdgCode);
1271           
1272           ptMC      = mcTrack->Pt();
1273           
1274           TParticle* mother = FindPrimaryMother(fMCStack, label);
1275           pdgMother = mother->GetPdgCode();
1276         }
1277       }
1278     
1279     
1280       //TOF
1281       Float_t beta = -99;
1282       if (esdTrack->GetStatus()&AliESDtrack::kTOFpid){
1283         if ((esdTrack->GetIntegratedLength() != 0) && (esdTrack->GetTOFsignal() != 0))
1284           beta = esdTrack->GetIntegratedLength()/esdTrack->GetTOFsignal()/fgkClight;
1285       }
1286       
1287       if(fTreeOption) {
1288         
1289         DeDxTrack* tracktpc = new((*fTrackArrayTPCPar)[nadded]) DeDxTrack();
1290         nadded++;
1291         
1292         tracktpc->p          = p;
1293         tracktpc->pt         = pt;
1294         //        track->ptcon   = pt_con;
1295         //        track->tpcchi  = tpcchi;
1296         tracktpc->eta        = eta;
1297         tracktpc->phi        = phi;
1298         tracktpc->q          = charge;
1299         tracktpc->filter     = 1;
1300         tracktpc->ncl        = ncl;
1301         tracktpc->neff       = neff;
1302         tracktpc->dedx       = dedx;
1303         tracktpc->beta       = beta;//computed with Global tracks
1304         tracktpc->dcaxy      = dcaxy;
1305         tracktpc->dcaz       = dcaz;
1306         tracktpc->pid        = pidCode;
1307         tracktpc->primary    = primaryFlag;
1308         tracktpc->pttrue     = ptMC;
1309         tracktpc->mother     = pdgMother;
1310         tracktpc->filterset1 = 0;
1311         tracktpc->filterset2 = 0;
1312         tracktpc->filterset3 = 0;
1313         
1314       }
1315
1316
1317       if(trackTPC)
1318         delete trackTPC;
1319
1320
1321     }//end of track loop
1322  
1323   }//end of: if isglobal==kFALSE
1324
1325
1326   if(fTreeOption) {
1327
1328     if( analysisMode==kGlobalTrk ){
1329       Sort(fTrackArrayGlobalPar, kFALSE);
1330       
1331       fEvent->trackmult = trackmult;
1332       fEvent->n         = nadded;
1333     }
1334     else if( analysisMode==kTPCTrk ){
1335       Sort(fTrackArrayTPCPar, kFALSE);
1336     }
1337
1338   }
1339
1340
1341 }
1342 //__________________________________________________________________
1343 void AliAnalysisTaskHighPtDeDx::ProduceArrayTrksAOD( AliAODEvent *AODevent, AnalysisMode analysisMode ){
1344
1345
1346   Int_t     trackmult = 0; // no pt cuts
1347   Int_t     nadded    = 0;
1348   Int_t nAODTracks = AODevent->GetNumberOfTracks();
1349   if( analysisMode == kGlobalTrk ){
1350     
1351     if(fTrackArrayGlobalPar)
1352       fTrackArrayGlobalPar->Clear();
1353     
1354   } 
1355   if( analysisMode == kTPCTrk ){
1356     if(fTrackArrayTPCPar)
1357       fTrackArrayTPCPar->Clear();
1358   }
1359
1360
1361
1362
1363   //const AliAODVertex* vertexSPD= (AliAODVertex*)AODevent->GetPrimaryVertexSPD();//GetPrimaryVertex()
1364   //if( vertexSPD->GetNContributors() < 1 || TMath::Abs( vertexSPD->GetZ()) > 10.0 ) return;
1365
1366
1367   
1368   if( analysisMode==kGlobalTrk ){  
1369  
1370
1371      const AliAODVertex* vertex = AODevent->GetPrimaryVertex();
1372     for(Int_t iT = 0; iT < nAODTracks; iT++) {
1373       
1374       AliAODTrack* aodTrack = AODevent->GetTrack(iT);
1375       
1376  
1377       
1378       UShort_t filterFlag = 0;
1379       Bool_t filterCut_Set1 = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
1380       Bool_t filterCut_Set2 = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
1381       Bool_t filterCut_Set3 = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
1382   
1383       
1384       
1385       if (fTrackFilterGolden) {
1386         
1387         // ITSTPC2010 cuts is bit 32 according to above macro, new versions of aliroot includes the golden cuts
1388         if(aodTrack->TestFilterBit(32)) {
1389           filterFlag +=1;
1390           filterCut_Set3 = kTRUE;
1391         }
1392       }
1393       
1394       
1395       if (fTrackFilterTPC) {
1396         // TPC only cuts is bit 1 according to above macro
1397         // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX, 
1398         if(aodTrack->TestFilterBit(1)){
1399           filterFlag +=2;
1400           filterCut_Set1 = kTRUE;
1401           trackmult++;
1402
1403         }
1404       }
1405       
1406       
1407       
1408       if(filterFlag==0)
1409         continue;
1410       
1411       
1412       Double_t b[2], cov[3];
1413       if (!(aodTrack->PropagateToDCA(vertex, AODevent->GetMagneticField(), kVeryBig, b, cov)))
1414         filterFlag = 32; // propagation failed!!!!!;
1415       Float_t dcaxy   = b[0];
1416       Float_t dcaz    = b[1];
1417       
1418
1419
1420
1421
1422       // As I understand this routine recalculates the momentum so it should be called first!
1423       //Double_t b[2], cov[3];
1424   
1425       
1426       //if(!aodTrack->PropagateToDCA(vertex, AODevent->GetMagneticField(), kVeryBig, b, cov))
1427       //        filterFlag = 32; // propagation failed!!!!!
1428       
1429       if(TMath::Abs(aodTrack->Eta()) > fEtaCut)
1430         continue;
1431       if (aodTrack->Pt() < fMinPt) {
1432         
1433         // Keep small fraction of low pT tracks
1434         if(fRandom->Rndm() > fLowPtFraction)
1435           continue;
1436       } // else {
1437       // Here we want to add the high pt part of the histograms!
1438       //    }
1439      
1440      
1441
1442       
1443       Short_t charge  = aodTrack->Charge();
1444       Float_t pt      = aodTrack->Pt();
1445       Float_t p       = aodTrack->P(); 
1446       Float_t eta     = aodTrack->Eta();
1447       Float_t phi     = aodTrack->Phi();
1448       AliAODPid* aodPid = aodTrack->GetDetPid();
1449       Short_t ncl     = -10;
1450       Short_t neff    = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1451       //          Short_t nclf    = aodTrack->GetTPCNclsF();
1452       Float_t dedx    = -10;
1453       Float_t beta = -99;
1454       if(aodPid) {
1455         ncl     = aodPid->GetTPCsignalN();
1456         dedx    = aodPid->GetTPCsignal();
1457         //TOF
1458         if (aodTrack->GetStatus()&AliESDtrack::kTOFpid){
1459           Double_t tof[5];
1460           aodPid->GetIntegratedTimes(tof);
1461           beta = tof[0]/aodPid->GetTOFsignal();
1462         }
1463       }
1464       //        Float_t tpcchi = aodTrack->Chi2perNDF();
1465       
1466       // Double_t p_con[3] = {0, 0, 0};
1467       // aodTrack->GetConstrainedPxPyPz(p_con);
1468       //        Float_t pt_con = 0; // This is not there! (Float_t)TMath::Sqrt(p_con[0]*p_con[0] + p_con[1]*p_con[1]);
1469       // const AliExternalTrackParam* tpcParam = aodTrack->GetTPCInnerParam();
1470       // Float_t pttpc   = tpcParam->Pt();
1471       // Float_t ptpc    = tpcParam->P();
1472       
1473       Float_t ptMC        = 0;
1474       Short_t pidCode     = 0; // 0 = real data / no mc track!
1475       Short_t primaryFlag = 0; // 0 = real data / not primary mc track  
1476       Int_t   pdgMother   = 0;
1477       
1478       if(fAnalysisMC) {
1479         
1480         const Int_t label = TMath::Abs(aodTrack->GetLabel());
1481         
1482         AliAODMCParticle* mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(label));
1483         if (mcTrack){
1484           
1485           if(mcTrack->IsPhysicalPrimary())
1486             primaryFlag = 1;
1487           
1488           Int_t pdgCode = mcTrack->GetPdgCode();
1489           pidCode = GetPidCode(pdgCode);
1490           
1491           ptMC      = mcTrack->Pt();
1492           
1493           AliAODMCParticle* mother = FindPrimaryMotherAOD(mcTrack);
1494           pdgMother = mother->GetPdgCode();         
1495         }
1496       }
1497     
1498     
1499       if(fTreeOption) {
1500         
1501         DeDxTrack* track = new((*fTrackArrayGlobalPar)[nadded]) DeDxTrack();
1502         nadded++;
1503         
1504         track->p          = p;
1505         track->pt         = pt;
1506         //        track->ptcon   = pt_con;
1507         //        track->tpcchi  = tpcchi;
1508         track->eta        = eta;
1509         track->phi        = phi;
1510         track->q          = charge;
1511         track->filter     = filterFlag;
1512         track->ncl        = ncl;
1513         track->neff       = neff;
1514         track->dedx       = dedx;
1515         track->beta       = beta;
1516         track->dcaxy      = dcaxy;
1517         track->dcaz       = dcaz;
1518         track->pid        = pidCode;
1519         track->primary    = primaryFlag;
1520         track->pttrue     = ptMC;
1521         track->mother     = pdgMother;
1522         track->filterset1 = filterCut_Set1;
1523         track->filterset2 = filterCut_Set2;
1524         track->filterset3 = filterCut_Set3;
1525       }
1526     }//end of track loop
1527
1528  
1529
1530   }//end of global track analysis
1531   
1532
1533
1534   else if( analysisMode==kTPCTrk ){  
1535  
1536     const AliAODVertex* vertexSPD= (AliAODVertex*)AODevent->GetPrimaryVertexSPD();//GetPrimaryVertex()
1537     if( vertexSPD->GetNContributors() < 1 || TMath::Abs(vertexSPD->GetZ()) > 10.0 ) return;
1538
1539
1540     //const AliAODVertex* vertex = AODevent->GetPrimaryVertex();
1541     for(Int_t iT = 0; iT < nAODTracks; iT++) {
1542       
1543       AliAODTrack* aodTrack = AODevent->GetTrack(iT);
1544       
1545  
1546       
1547       UShort_t filterFlag = 0;
1548       Bool_t filterCut_Set1 = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
1549       Bool_t filterCut_Set2 = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
1550       Bool_t filterCut_Set3 = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
1551   
1552       
1553      
1554       // TPC only cuts is bit 1 according to above macro
1555       // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX, 
1556       if(aodTrack->TestFilterBit(128)){
1557         filterFlag +=2;
1558         trackmult++;
1559       }
1560       
1561       
1562       if(filterFlag==0)
1563         continue;
1564       
1565       
1566       Double_t b[2], cov[3];
1567       //AliAODVertex* vertex = AODevent->GetPrimaryVertex();
1568       //Double_t b[2] = {-99., -99.};
1569       //Double_t bCov[3] = {-99., -99., -99.};
1570       //aodTrack->PropagateToDCA(aod->GetPrimaryVertex(), aod->GetMagneticField(), 100., b, bCov);
1571       if (!(aodTrack->PropagateToDCA(vertexSPD, AODevent->GetMagneticField(), kVeryBig, b, cov)))
1572         filterFlag = 32; // propagation failed!!!!!;
1573       Float_t dcaxy   = b[0];
1574       Float_t dcaz    = b[1];
1575
1576
1577
1578
1579       // As I understand this routine recalculates the momentum so it should be called first!
1580       //Double_t b[2], cov[3];
1581   
1582       
1583       //if(!aodTrack->PropagateToDCA(vertex, AODevent->GetMagneticField(), kVeryBig, b, cov))
1584       //        filterFlag = 32; // propagation failed!!!!!
1585       
1586       if(TMath::Abs(aodTrack->Eta()) > fEtaCut)
1587         continue;
1588       if (aodTrack->Pt() < fMinPt) {
1589         
1590         // Keep small fraction of low pT tracks
1591         if(fRandom->Rndm() > fLowPtFraction)
1592           continue;
1593       } // else {
1594       // Here we want to add the high pt part of the histograms!
1595       //    }
1596      
1597      
1598
1599       
1600       Short_t charge  = aodTrack->Charge();
1601       Float_t pt      = aodTrack->Pt();
1602       Float_t p       = aodTrack->P(); 
1603       Float_t eta     = aodTrack->Eta();
1604       Float_t phi     = aodTrack->Phi();
1605       AliAODPid* aodPid = aodTrack->GetDetPid();
1606       Short_t ncl     = -10;
1607       Short_t neff    = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
1608       //          Short_t nclf    = aodTrack->GetTPCNclsF();
1609       Float_t dedx    = -10;
1610       Float_t beta = -99;
1611       if(aodPid) {
1612         ncl     = aodPid->GetTPCsignalN();
1613         dedx    = aodPid->GetTPCsignal();
1614         //TOF
1615         if (aodTrack->GetStatus()&AliESDtrack::kTOFpid){
1616           Double_t tof[5];
1617           aodPid->GetIntegratedTimes(tof);
1618           beta = tof[0]/aodPid->GetTOFsignal();
1619         }
1620       }
1621       //        Float_t tpcchi = aodTrack->Chi2perNDF();
1622       
1623       // Double_t p_con[3] = {0, 0, 0};
1624       // aodTrack->GetConstrainedPxPyPz(p_con);
1625       //        Float_t pt_con = 0; // This is not there! (Float_t)TMath::Sqrt(p_con[0]*p_con[0] + p_con[1]*p_con[1]);
1626       // const AliExternalTrackParam* tpcParam = aodTrack->GetTPCInnerParam();
1627       // Float_t pttpc   = tpcParam->Pt();
1628       // Float_t ptpc    = tpcParam->P();
1629       
1630       Float_t ptMC        = 0;
1631       Short_t pidCode     = 0; // 0 = real data / no mc track!
1632       Short_t primaryFlag = 0; // 0 = real data / not primary mc track  
1633       Int_t   pdgMother   = 0;
1634       
1635       if(fAnalysisMC) {
1636         
1637         const Int_t label = TMath::Abs(aodTrack->GetLabel());
1638         
1639         AliAODMCParticle* mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(label));
1640         if (mcTrack){
1641           
1642           if(mcTrack->IsPhysicalPrimary())
1643             primaryFlag = 1;
1644           
1645           Int_t pdgCode = mcTrack->GetPdgCode();
1646           pidCode = GetPidCode(pdgCode);
1647           
1648           ptMC      = mcTrack->Pt();
1649           
1650           AliAODMCParticle* mother = FindPrimaryMotherAOD(mcTrack);
1651           pdgMother = mother->GetPdgCode();         
1652         }
1653       }
1654     
1655     
1656       if(fTreeOption) {
1657         
1658         DeDxTrack* tracktpc = new((*fTrackArrayTPCPar)[nadded]) DeDxTrack();
1659         nadded++;
1660         
1661         tracktpc->p          = p;
1662         tracktpc->pt         = pt;
1663         //        track->ptcon   = pt_con;
1664         //        track->tpcchi  = tpcchi;
1665         tracktpc->eta        = eta;
1666         tracktpc->phi        = phi;
1667         tracktpc->q          = charge;
1668         tracktpc->filter     = filterFlag;
1669         tracktpc->ncl        = ncl;
1670         tracktpc->neff       = neff;
1671         tracktpc->dedx       = dedx;
1672         tracktpc->beta       = beta;
1673         tracktpc->dcaxy      = dcaxy;
1674         tracktpc->dcaz       = dcaz;
1675         tracktpc->pid        = pidCode;
1676         tracktpc->primary    = primaryFlag;
1677         tracktpc->pttrue     = ptMC;
1678         tracktpc->mother     = pdgMother;
1679         tracktpc->filterset1 = filterCut_Set1;
1680         tracktpc->filterset2 = filterCut_Set2;
1681         tracktpc->filterset3 = filterCut_Set3;
1682       }
1683     }//end of track loop
1684
1685
1686   }//end of global track analysis
1687
1688
1689   if(fTreeOption) {
1690     
1691     if( analysisMode==kGlobalTrk ){
1692       Sort(fTrackArrayGlobalPar, kFALSE);
1693       
1694       fEvent->trackmult = trackmult;
1695       fEvent->n         = nadded;
1696
1697
1698     }
1699     if( analysisMode==kTPCTrk ){
1700       Sort(fTrackArrayTPCPar, kFALSE);
1701     }
1702     
1703   }
1704   
1705 }