Adapt add macro and particle/cluster containers for the analysis on jets
[u/mrichter/AliRoot.git] / PWGJE / AliPWG4HighPtSpectra.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 //-----------------------------------------------------------------------
17 // Efficiency between different steps of the procedure.
18 // The ouptut of the task is a AliCFContainer from which the efficiencies
19 // can be calculated
20 //-----------------------------------------------------------------------
21 // Author : Marta Verweij - UU
22 //-----------------------------------------------------------------------
23
24
25 #ifndef ALIPWG4HIGHPTSPECTRA_CXX
26 #define ALIPWG4HIGHPTSPECTRA_CXX
27
28 #include "AliPWG4HighPtSpectra.h"
29
30 #include "TVector3.h"
31 #include <iostream>
32 #include "TH1.h"
33 #include "TH2.h"
34 #include "TH3.h"
35 #include "TProfile.h"
36 #include "TList.h"
37 #include "TChain.h"
38 #include "TKey.h"
39 #include "TSystem.h"
40 #include "TFile.h"
41
42 #include "AliAnalysisManager.h"
43 #include "AliESDInputHandler.h"
44 #include "AliAODInputHandler.h"
45 #include "AliESDtrack.h"
46 #include "AliAODTrack.h"
47 #include "AliAODMCParticle.h"
48 #include "AliESDtrackCuts.h"
49 #include "AliExternalTrackParam.h"
50 #include "AliCentrality.h"
51
52 #include "AliLog.h"
53
54 #include "AliStack.h"
55 #include "TParticle.h"
56 #include "AliMCEvent.h"
57 #include "AliMCEventHandler.h"
58 #include "AliCFContainer.h"
59 #include "AliGenPythiaEventHeader.h"
60 #include "AliGenHijingEventHeader.h"
61 #include "AliGenCocktailEventHeader.h"
62 #include "AliAODMCHeader.h"
63 #include "AliESDVertex.h"
64 //#include "$ALICE_ROOT/PWG4/JetTasks/AliAnalysisHelperJetTasks.h"
65
66 //#include <iostream>
67 using namespace std; //required for resolving the 'cout' symbol
68
69 ClassImp(AliPWG4HighPtSpectra)
70
71 //__________________________________________________________________________
72 AliPWG4HighPtSpectra::AliPWG4HighPtSpectra() : AliAnalysisTask("AliPWG4HighPtSpectra", ""), 
73   fReadAODData(0),
74   fNoPythiaInfo(0),
75   fCFManagerPos(0x0),
76   fCFManagerNeg(0x0),
77   fESD(0x0),
78   fAOD(0x0),
79   fMC(0x0),
80   fStack(0x0),
81   fArrayMCAOD(0x0),
82   fVtx(0x0),
83   fTriggerMask(AliVEvent::kMB),
84   fIsPbPb(0),
85   fCentClass(10),
86   fTrackType(0),
87   fTrackCuts(0x0),
88   fTrackCutsReject(0x0),
89   fFilterMask(0),
90   fbSelectHIJING(kFALSE),
91   fSigmaConstrainedMax(100.),
92   fAvgTrials(1),
93   fHistList(0),
94   fNEventAll(0),
95   fNEventSel(0),
96   fNEventReject(0),
97   fh1Centrality(0x0),
98   fh1Xsec(0),
99   fh1Trials(0),
100   fh1PtHard(0),
101   fh1PtHardTrials(0),
102   fPtRelUncertainty1PtPrim(0x0),
103   fPtRelUncertainty1PtSec(0x0)
104 {
105   //
106   //Default ctor
107   //
108 }
109
110 //___________________________________________________________________________
111 AliPWG4HighPtSpectra::AliPWG4HighPtSpectra(const Char_t* name) :
112   AliAnalysisTask(name,""),
113   fReadAODData(0),
114   fNoPythiaInfo(0),
115   fCFManagerPos(0x0),
116   fCFManagerNeg(0x0),
117   fESD(0x0),
118   fAOD(0x0),
119   fMC(0x0),
120   fStack(0x0),
121   fArrayMCAOD(0x0),
122   fVtx(0x0),
123   fTriggerMask(AliVEvent::kMB),
124   fIsPbPb(0),
125   fCentClass(10),
126   fTrackType(0),
127   fTrackCuts(0x0),
128   fTrackCutsReject(0x0),
129   fFilterMask(0),
130   fbSelectHIJING(kFALSE),
131   fSigmaConstrainedMax(100.),
132   fAvgTrials(1),
133   fHistList(0),
134   fNEventAll(0),
135   fNEventSel(0),
136   fNEventReject(0),
137   fh1Centrality(0x0),
138   fh1Xsec(0),
139   fh1Trials(0),
140   fh1PtHard(0),
141   fh1PtHardTrials(0),
142   fPtRelUncertainty1PtPrim(0x0),
143   fPtRelUncertainty1PtSec(0x0)
144 {
145   //
146   // Constructor. Initialization of Inputs and Outputs
147   //
148   AliDebug(2,Form("AliPWG4HighPtSpectra Calling Constructor"));
149   // Input slot #0 works with a TChain ESD
150   DefineInput(0, TChain::Class());
151   // Output slot #0 writes into a TList
152   DefineOutput(0,TList::Class());
153   // Output slot #1, #2 writes into a AliCFContainer
154   DefineOutput(1,AliCFContainer::Class());
155   DefineOutput(2,AliCFContainer::Class());
156   // Output slot #3 writes into a AliESDtrackCuts
157   DefineOutput(3, AliESDtrackCuts::Class());
158   DefineOutput(4, AliESDtrackCuts::Class());
159 }
160
161 //________________________________________________________________________
162 void AliPWG4HighPtSpectra::LocalInit()
163 {
164   //
165   // Only called once at beginning
166   //
167   PostData(3,fTrackCuts);
168 }
169
170 //________________________________________________________________________
171 void AliPWG4HighPtSpectra::ConnectInputData(Option_t *) 
172 {
173   // Connect ESD here
174   // Called once
175   AliDebug(2,Form(">> AliPWG4HighPtSpectra::ConnectInputData \n"));
176
177   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
178   if (!tree) {
179     AliDebug(2,Form( "ERROR: Could not read chain from input slot 0 \n"));
180     return;
181   }
182
183   if(fReadAODData) {
184     AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
185   
186     if (!aodH) {
187       AliDebug(2,Form("ERROR: Could not get AODInputHandler"));
188       return;
189     } else
190       fAOD = aodH->GetEvent();
191   } else {
192     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
193   
194     if (!esdH) {
195       AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
196       return;
197     } else
198       fESD = esdH->GetEvent();
199   }
200
201   AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
202   if (!eventHandler) {
203     AliDebug(2,Form( "ERROR: Could not retrieve MC event handler \n"));
204   }
205   else
206     fMC = eventHandler->MCEvent();
207
208 }
209
210 //________________________________________________________________________
211 Bool_t AliPWG4HighPtSpectra::SelectEvent()
212 {
213   //
214   // Decide if event should be selected for analysis
215   //
216
217   // Checks following requirements:
218   // - fESD available
219   // - trigger info from AliPhysicsSelection
220   // - number of reconstructed tracks > 1
221   // - primary vertex reconstructed
222   // - z-vertex < 10 cm
223
224   Bool_t selectEvent = kTRUE;
225
226   //fESD or fAOD object available?
227   if (!fESD && !fAOD) {
228     AliDebug(2,Form("ERROR: fInputEvent not available\n"));
229     fNEventReject->Fill("noESD/AOD",1);
230     selectEvent = kFALSE;
231     return selectEvent;
232   }
233
234   //Trigger
235   UInt_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
236   if(fTriggerMask != AliVEvent::kAny && !(isSelected&fTriggerMask)) { //Select collison candidates
237     AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
238     fNEventReject->Fill("Trigger",1);
239     selectEvent = kFALSE;
240     return selectEvent;
241   }
242
243   //Check if number of reconstructed tracks is larger than 1
244   if(!fReadAODData) {
245     if(!fESD->GetNumberOfTracks() || fESD->GetNumberOfTracks()<2)  {
246       fNEventReject->Fill("NTracks<2",1);
247       selectEvent = kFALSE;
248       return selectEvent;
249     }
250   }
251
252   //Check if vertex is reconstructed
253   if(fESD)
254     fVtx = fESD->GetPrimaryVertexSPD();
255   if(fAOD)
256     fVtx = fAOD->GetPrimaryVertex();
257
258   if(!fVtx) {
259     fNEventReject->Fill("noVTX",1);
260     selectEvent = kFALSE;
261     return selectEvent;
262   }
263
264   if(!fReadAODData){
265     const AliESDVertex* esdVtx = fESD->GetPrimaryVertexSPD();
266     if(!esdVtx->GetStatus()) {
267       fNEventReject->Fill("VtxStatus",1);
268       selectEvent = kFALSE;
269       return selectEvent;
270     }
271   }
272
273   // Need vertex cut
274   //  TString vtxName(fVtx->GetName());
275   if(fVtx->GetNContributors()<2) {
276     fNEventReject->Fill("NCont<2",1);
277     selectEvent = kFALSE;
278     return selectEvent;
279   }
280
281   //Check if z-vertex < 10 cm
282   double primVtx[3];
283   fVtx->GetXYZ(primVtx);
284   if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
285     fNEventReject->Fill("ZVTX>10",1);
286     selectEvent = kFALSE;
287     return selectEvent;
288   }
289
290   //Centrality selection should only be done in case of PbPb
291   if(IsPbPb()) {
292     Float_t cent = 0.;
293     Int_t calccent = 0;
294     if(fReadAODData)
295       calccent = CalculateCentrality(fAOD);
296     else
297       calccent = CalculateCentrality(fESD);
298     if(fCentClass!=calccent && fCentClass!=10) {
299       fNEventReject->Fill("cent",1);
300       selectEvent = kFALSE;
301       return selectEvent;
302     }
303     else {
304       if(fReadAODData) {
305         if(dynamic_cast<AliAODEvent*>(fAOD)->GetCentrality()) {
306           cent = dynamic_cast<AliAODEvent*>(fAOD)->GetCentrality()->GetCentralityPercentile("V0M");
307         }
308       }
309       else {
310         if(dynamic_cast<AliESDEvent*>(fESD)->GetCentrality()) {
311           cent = dynamic_cast<AliESDEvent*>(fESD)->GetCentrality()->GetCentralityPercentile("V0M");
312         }
313       }
314
315       if(cent>90.) {
316         fNEventReject->Fill("cent>90",1);
317         selectEvent = kFALSE;
318         return selectEvent;
319       }
320       fh1Centrality->Fill(cent);
321     }
322   }
323
324   return selectEvent;
325
326 }
327
328 //________________________________________________________________________
329 Int_t AliPWG4HighPtSpectra::CalculateCentrality(AliVEvent *event)
330 {
331   Float_t cent = 999;
332
333   if(event){
334     if(event->GetCentrality()){
335       cent = event->GetCentrality()->GetCentralityPercentile("V0M");
336     }
337   }
338
339   if(cent<0) return 5;
340   if(cent>80)return 4;
341   if(cent>50)return 3;
342   if(cent>30)return 2;
343   if(cent>10)return 1;
344   return 0;
345 }
346
347 //_________________________________________________
348 void AliPWG4HighPtSpectra::Exec(Option_t *)
349 {
350   //
351   // Main loop function
352   //
353   AliDebug(2,Form(">> AliPWG4HighPtSpectra::Exec \n"));  
354
355   AliVEvent* event;
356   if(fReadAODData)
357     event = fAOD;
358   else {
359     event = fESD;
360     if(!fMC) {
361       AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
362       if (!eventHandler) {
363         AliDebug(2,Form( "ERROR: Could not retrieve MC event handler \n"));
364       }
365       else
366         fMC = eventHandler->MCEvent();
367     }
368   }
369
370   // All events without selection
371   fNEventAll->Fill(0.);
372
373   if(!SelectEvent()) {
374     // Post output data
375     PostData(0,fHistList);
376     PostData(1,fCFManagerPos->GetParticleContainer());
377     PostData(2,fCFManagerNeg->GetParticleContainer());
378     return;
379   }
380
381   if(fReadAODData){
382     //Open the MC particles in the case of AODanalysis
383     fArrayMCAOD = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
384   }
385   else {
386     //In case of ESD event: MCEvent available? if yes: get MC particles stack
387     if(fMC) {
388       AliDebug(2,Form("MC particles: %d", fMC->GetNumberOfTracks()));
389       fStack = fMC->Stack();                //Particles Stack
390       AliDebug(2,Form("MC particles stack: %d", fStack->GetNtrack()));
391     }
392   }
393
394   Int_t nTracks = event->GetNumberOfTracks();
395   AliDebug(2,Form("nTracks %d", nTracks));
396
397   if((!fTrackCuts && !fReadAODData) || (fFilterMask==0 && fReadAODData)) { //if the TrackCuts are missing in ESD analysis or the FilterBit is missing in AOD analysis
398     fNEventReject->Fill("noTrackCuts",1);
399     // Post output data
400     PostData(0,fHistList);
401     PostData(1,fCFManagerPos->GetParticleContainer());
402     PostData(2,fCFManagerNeg->GetParticleContainer());
403     return;
404   }
405
406   // Selected events for analysis
407   fNEventSel->Fill(0.);
408   
409   const Int_t nvar = 4;
410   
411   Double_t containerInputRec[nvar]       = {0.,0.,0.,0.};
412   Double_t containerInputMC[nvar]        = {0.,0.,0.,0.};
413   Double_t containerInputRecMC[nvar]     = {0.,0.,0.,0.}; //reconstructed yield as function of MC variable
414
415   //Now go to rec level
416   if(fReadAODData) {//AOD analysis
417     for (Int_t iTrack = 0; iTrack<nTracks; iTrack++) 
418     {
419       //Get track for analysis
420       AliVTrack *track = 0x0;
421
422       AliAODTrack *aodtrack = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTrack));
423       if(!aodtrack) AliFatal("Not a standard AOD");
424       if(!aodtrack)
425         continue;
426       if((!aodtrack->TestFilterBit(fFilterMask)) || 
427          ((!fCFManagerPos->CheckParticleCuts(kStepReconstructed,aodtrack)) && (aodtrack->Charge()>0.)) ||
428          ((!fCFManagerNeg->CheckParticleCuts(kStepReconstructed,aodtrack)) && (aodtrack->Charge()<0))    )
429         continue;
430       else {
431         track = static_cast<AliVTrack*>(aodtrack);
432       }
433       //fill the container
434       containerInputRec[0] = track->Pt();
435       containerInputRec[1] = track->Phi();
436       containerInputRec[2] = track->Eta();
437       containerInputRec[3] = track->GetTPCNcls();
438
439       if(track->Charge()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
440       if(track->Charge()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
441
442       if(fArrayMCAOD) {
443         Int_t label = TMath::Abs(track->GetLabel());
444         if(label>fArrayMCAOD->GetEntries()) {
445           continue;
446         }
447         AliAODMCParticle *particle = (AliAODMCParticle*) fArrayMCAOD->At(label);
448         if(!particle) {
449           continue;
450         }
451         //Only select particles generated by HIJING if requested
452         if(fbSelectHIJING) {
453           if(!IsHIJINGParticle(label)) {
454             continue;
455           }
456         }
457         containerInputRecMC[0] = particle->Pt();      
458         containerInputRecMC[1] = particle->Phi();      
459         containerInputRecMC[2] = particle->Eta();  
460         containerInputRecMC[3] = track->GetTPCNcls();
461
462         //Container with primaries
463         if(particle->IsPhysicalPrimary()) {
464           if(particle->Charge()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
465           if(particle->Charge()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
466           //Fill pT resolution plots for primaries
467           //fPtRelUncertainty1PtPrim->Fill(containerInputRec[0],containerInputRec[0]*TMath::Sqrt(track->GetSigma1Pt2())); //This has not been implemented in AOD analysis, since they are also produced by the AddTaskPWG4HighPtTrackQA.C macro
468         }
469
470         //Container with secondaries
471         if (!particle->IsPhysicalPrimary() ) {
472           if(particle->Charge()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRecMC,kStepSecondaries);
473           if(particle->Charge()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRecMC,kStepSecondaries);
474           //Fill pT resolution plots for primaries
475           //fPtRelUncertainty1PtSec->Fill(containerInputRec[0],containerInputRec[0]*TMath::Sqrt(track->GetSigma1Pt2())); //This has not been implemented in AOD analysis, since they are also produced by the AddTaskPWG4HighPtTrackQA.C macro
476         }
477       }
478     }//track loop
479
480     //Fill MC containers if particles are findable
481     if(fArrayMCAOD) {
482       int noPart = fArrayMCAOD->GetEntriesFast();
483
484       for(int iPart = 1; iPart<noPart; iPart++) {
485         AliAODMCParticle *mcPart = (AliAODMCParticle*) fArrayMCAOD->At(iPart);
486         if(!mcPart) continue;
487
488         //Only select particles generated by HIJING if requested
489         if(fbSelectHIJING) {
490           if(!IsHIJINGParticle(iPart))
491             continue;
492         }
493
494         Int_t pdg = mcPart->PdgCode();
495
496         // select charged pions, protons, kaons , electrons, muons
497         if(TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) ==
498            321 || TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13){
499   
500           //fill the container
501           containerInputMC[0] = mcPart->Pt();
502           containerInputMC[1] = mcPart->Phi();      
503           containerInputMC[2] = mcPart->Eta();  
504           //      AliESDtrack *esdtrack = fESD->GetTrack(mcPart->GetLabel());
505           containerInputMC[3] = 159.;
506           
507           if(mcPart->IsPhysicalPrimary()) {
508             if(mcPart->Charge()>0. && fCFManagerPos->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
509             if(mcPart->Charge()<0. && fCFManagerNeg->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
510           }
511         }
512       }
513     }
514
515   }//end of AOD analysis
516   else {//ESD analysis
517     for (Int_t iTrack = 0; iTrack<nTracks; iTrack++) 
518     {
519       //Get track for analysis
520       AliESDtrack *track = 0x0;
521       AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
522       if(!esdtrack)
523         continue;
524       
525
526       if(fTrackType==4) {
527         if (!(fTrackCuts->AcceptTrack(esdtrack)))
528           continue;
529       }
530
531       if(fTrackType==1)
532         track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
533       else if(fTrackType==2 || fTrackType==4) {
534         track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(fESD),esdtrack->GetID());
535         if(!track)
536           continue;
537         
538         AliExternalTrackParam exParam;
539         Bool_t relate = track->RelateToVertexTPC(fESD->GetPrimaryVertexSPD(),fESD->GetMagneticField(),kVeryBig,&exParam);
540         if( !relate ) {
541           if(track) delete track;
542           continue;
543         }
544         track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
545       }
546       else if(fTrackType==5 || fTrackType==6) {
547         if(fTrackCuts->AcceptTrack(esdtrack)) {
548           continue;
549         }
550         else {
551           if( !(fTrackCutsReject->AcceptTrack(esdtrack)) && fTrackCuts->AcceptTrack(esdtrack) ) {
552
553             if(fTrackType==5) {
554               //use TPConly constrained track
555               track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
556               if(!track)
557                 continue;
558               
559               AliExternalTrackParam exParam;
560               Bool_t relate = track->RelateToVertexTPC(fESD->GetPrimaryVertexSPD(),fESD->GetMagneticField(),kVeryBig,&exParam);
561               if( !relate ) {
562                 if(track) delete track;
563                 continue;
564               }
565               track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
566             }
567             else if(fTrackType==6) {
568               //use global constrained track
569               track = new AliESDtrack(*esdtrack);
570               track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
571             }
572           }
573         }
574       }
575       else if(fTrackType==7) {
576         //use global constrained track
577         track = new AliESDtrack(*esdtrack);
578       }
579       else
580         track = esdtrack;
581     
582       if(!track)
583         continue;
584       
585
586       if(fTrackType==2 || fTrackType==4 || fTrackType==5) {
587         //Cut on chi2 of constrained fit
588         if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax && fSigmaConstrainedMax>0.) {
589           if(track) delete track;
590           continue;
591         }
592       }
593       if (!(fTrackCuts->AcceptTrack(track)) && fTrackType!=4 && fTrackType!=5 && fTrackType!=6) {
594         if(fTrackType==1 || fTrackType==2 || fTrackType==7) {
595           if(track) delete track;
596         }
597         continue;
598       }
599
600       if(fTrackType==7) {
601         if(fTrackCutsReject ) {
602           if(fTrackCutsReject->AcceptTrack(track) ) {
603             if(track) delete track;
604             continue;
605           }
606         }
607         if(esdtrack->GetConstrainedParam()) 
608           track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
609       }
610
611       if(!track) {
612         if(fTrackType==1 || fTrackType==2 || fTrackType==4 || fTrackType==5 || fTrackType==6 || fTrackType==7) {
613           if(track) delete track;
614         }
615         continue;
616       }
617
618       //fill the container
619       containerInputRec[0] = track->Pt();
620       containerInputRec[1] = track->Phi();
621       containerInputRec[2] = track->Eta();
622       containerInputRec[3] = track->GetTPCNclsIter1();
623
624       if(track->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
625       if(track->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
626
627         
628       //Only fill the MC containers if MC information is available
629       if(fMC) {
630         Int_t label = TMath::Abs(track->GetLabel());
631         if(label>fStack->GetNtrack()) {
632           if(fTrackType==1 || fTrackType==2 || fTrackType==7)
633             delete track;
634           continue;
635         }
636         //Only select particles generated by HIJING if requested
637         if(fbSelectHIJING) {
638           if(!IsHIJINGParticle(label)) {
639             if(fTrackType==1 || fTrackType==2 || fTrackType==7)
640               delete track;
641             continue;
642           }
643         }
644         TParticle *particle = fStack->Particle(label) ;
645         if(!particle) {
646           if(fTrackType==1 || fTrackType==2 || fTrackType==7)
647             delete track;
648           continue;
649         }
650         containerInputRecMC[0] = particle->Pt();      
651         containerInputRecMC[1] = particle->Phi();      
652         containerInputRecMC[2] = particle->Eta();  
653         containerInputRecMC[3] = track->GetTPCNclsIter1();
654
655         //Container with primaries
656         if(fStack->IsPhysicalPrimary(label)) {
657           if(particle->GetPDG()->Charge()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
658           if(particle->GetPDG()->Charge()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
659           //Fill pT resolution plots for primaries
660           fPtRelUncertainty1PtPrim->Fill(containerInputRec[0],containerInputRec[0]*TMath::Sqrt(track->GetSigma1Pt2()));
661         }
662
663         //Container with secondaries
664         if (!fStack->IsPhysicalPrimary(label) ) {
665           if(particle->GetPDG()->Charge()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepSecondaries);
666           if(particle->GetPDG()->Charge()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepSecondaries);
667           //Fill pT resolution plots for primaries
668           fPtRelUncertainty1PtSec->Fill(containerInputRec[0],containerInputRec[0]*TMath::Sqrt(track->GetSigma1Pt2()));
669         }
670       }
671
672       if(fTrackType==1  || fTrackType==2 || fTrackType==4 || fTrackType==5 || fTrackType==6 || fTrackType==7) {
673         if(track) delete track;
674       }
675     }//track loop
676
677     //Fill MC containers if particles are findable
678     if(fMC) {
679       for(int iPart = 1; iPart<(fMC->GetNumberOfPrimaries()); iPart++) {
680         AliMCParticle *mcPart  = (AliMCParticle*)fMC->GetTrack(iPart);
681         if(!mcPart) continue;
682   
683         //Only select particles generated by HIJING if requested
684         if(fbSelectHIJING) {
685           if(!IsHIJINGParticle(iPart))
686             continue;
687         }
688   
689         Int_t pdg = mcPart->PdgCode();
690         
691         // select charged pions, protons, kaons , electrons, muons
692         if(TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) ==
693            321 || TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13){
694   
695           //fill the container
696           containerInputMC[0] = mcPart->Pt();
697           containerInputMC[1] = mcPart->Phi();      
698           containerInputMC[2] = mcPart->Eta();  
699           //      AliESDtrack *esdtrack = fESD->GetTrack(mcPart->GetLabel());
700           containerInputMC[3] = 159.;
701           
702           if(fStack->IsPhysicalPrimary(iPart)) {
703             if(mcPart->Charge()>0. && fCFManagerPos->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
704             if(mcPart->Charge()<0. && fCFManagerNeg->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
705           }
706         }
707       }
708     }
709   }//end of ESD analysis
710
711   PostData(0,fHistList);
712   PostData(1,fCFManagerPos->GetParticleContainer());
713   PostData(2,fCFManagerNeg->GetParticleContainer());
714   
715 }
716
717 //________________________________________________________________________
718 Bool_t AliPWG4HighPtSpectra::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials)
719 {
720   //
721   // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
722   // This is to called in Notify and should provide the path to the AOD/ESD file
723   // Copied from AliAnalysisTaskJetSpectrum2
724   //
725
726   TString file(currFile);  
727   fXsec = 0;
728   fTrials = 1;
729
730   if(file.Contains("root_archive.zip#")){
731     Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
732     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
733     file.Replace(pos+1,20,"");
734   }
735   else {
736     // not an archive take the basename....
737     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
738   }
739   //  Printf("%s",file.Data());
740    
741   TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root
742   if(!fxsec){
743     // next trial fetch the histgram file
744     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
745     if(!fxsec){
746       // not a severe condition but indicates that we have no information
747       return kFALSE;
748     }
749     else{
750       // find the tlist we want to be independtent of the name so use the Tkey
751       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0); 
752       if(!key){
753         fxsec->Close();
754         return kFALSE;
755       }
756       TList *list = dynamic_cast<TList*>(key->ReadObj());
757       if(!list){
758         fxsec->Close();
759         return kFALSE;
760       }
761       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
762       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
763       fxsec->Close();
764     }
765   } // no tree pyxsec.root
766   else {
767     TTree *xtree = (TTree*)fxsec->Get("Xsection");
768     if(!xtree){
769       fxsec->Close();
770       return kFALSE;
771     }
772     UInt_t   ntrials  = 0;
773     Double_t  xsection  = 0;
774     xtree->SetBranchAddress("xsection",&xsection);
775     xtree->SetBranchAddress("ntrials",&ntrials);
776     xtree->GetEntry(0);
777     fTrials = ntrials;
778     fXsec = xsection;
779     fxsec->Close();
780   }
781   return kTRUE;
782 }
783
784 //________________________________________________________________________
785 Bool_t AliPWG4HighPtSpectra::Notify()
786 {
787   //
788   // Implemented Notify() to read the cross sections
789   // and number of trials from pyxsec.root
790   // Copied from AliAnalysisTaskJetSpectrum2
791   // 
792
793   if(fNoPythiaInfo)
794     return kTRUE;
795
796   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
797   Float_t xsection = 0;
798   Float_t ftrials  = 1;
799
800   if(tree){
801     TFile *curfile = tree->GetCurrentFile();
802     if (!curfile) {
803       Error("Notify","No current file");
804       return kFALSE;
805     }
806     if(!fh1Xsec||!fh1Trials){
807       //      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
808       return kFALSE;
809     }
810     PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
811     fh1Xsec->Fill("<#sigma>",xsection);
812     fh1Trials->Fill("#sum{ntrials}",ftrials);
813   } 
814   return kTRUE;
815 }
816
817 //________________________________________________________________________
818 AliGenPythiaEventHeader*  AliPWG4HighPtSpectra::GetPythiaEventHeader(AliMCEvent *mcEvent)
819 {
820   
821   if(!mcEvent)return 0;
822   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
823   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
824   if(!pythiaGenHeader){
825     // cocktail ??
826     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
827     
828     if (!genCocktailHeader) {
829       //      AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
830       //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
831       return 0;
832     }
833     TList* headerList = genCocktailHeader->GetHeaders();
834     for (Int_t i=0; i<headerList->GetEntries(); i++) {
835       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
836       if (pythiaGenHeader)
837         break;
838     }
839     if(!pythiaGenHeader){
840       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
841       return 0;
842     }
843   }
844   return pythiaGenHeader;
845
846 }
847
848 //________________________________________________________________________
849 AliGenHijingEventHeader*  AliPWG4HighPtSpectra::GetHijingEventHeader(AliMCEvent *mcEvent)
850 {
851   
852   if(!mcEvent)return 0;
853   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
854   AliGenHijingEventHeader* hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
855   if(!hijingGenHeader){
856     // cocktail ??
857     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
858     
859     if (!genCocktailHeader) {
860       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
861       //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
862       return 0;
863     }
864     TList* headerList = genCocktailHeader->GetHeaders();
865     for (Int_t i=0; i<headerList->GetEntries(); i++) {
866       hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(headerList->At(i));
867       if (hijingGenHeader)
868         break;
869     }
870     if(!hijingGenHeader){
871       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Hijing event header not found");
872       return 0;
873     }
874   }
875   return hijingGenHeader;
876
877 }
878
879 //________________________________________________________________________
880 AliGenHijingEventHeader*  AliPWG4HighPtSpectra::GetHijingEventHeaderAOD(AliAODEvent *aodEvent)
881 {
882   AliGenHijingEventHeader* hijingGenHeader = 0x0;
883   if(aodEvent) {
884     AliAODMCHeader* mcHeader = (AliAODMCHeader*) aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
885     TList* headerList = mcHeader->GetCocktailHeaders();
886     for (Int_t i=0; i<headerList->GetEntries(); i++) {
887       hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(headerList->At(i));
888       if (hijingGenHeader)
889         break;
890     }
891   }
892   return hijingGenHeader;
893 }
894
895 //___________________________________________________________________________
896 void AliPWG4HighPtSpectra::Terminate(Option_t*)
897 {
898   // The Terminate() function is the last function to be called during
899   // a query. It always runs on the client, it can be used to present
900   // the results graphically or save the results to file.
901
902 }
903
904 //___________________________________________________________________________
905 void AliPWG4HighPtSpectra::CreateOutputObjects() 
906 {
907   //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
908   //TO BE SET BEFORE THE EXECUTION OF THE TASK
909   //
910   AliDebug(2,Form("CreateOutputObjects CreateOutputObjects of task %s", GetName()));
911
912   Bool_t oldStatus = TH1::AddDirectoryStatus();
913   TH1::AddDirectory(kFALSE); 
914
915   //slot #1
916   OpenFile(0);
917   fHistList = new TList();
918   fHistList->SetOwner(kTRUE);
919   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
920   fHistList->Add(fNEventAll);
921   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
922   fHistList->Add(fNEventSel);
923
924   fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
925   //Set labels
926   fNEventReject->Fill("noESD/AOD",0);
927   fNEventReject->Fill("Trigger",0);
928   fNEventReject->Fill("NTracks<2",0);
929   fNEventReject->Fill("noVTX",0);
930   fNEventReject->Fill("VtxStatus",0);
931   fNEventReject->Fill("NCont<2",0);
932   fNEventReject->Fill("ZVTX>10",0);
933   fNEventReject->Fill("cent",0);
934   fNEventReject->Fill("cent>90",0);
935   fHistList->Add(fNEventReject);
936
937   fh1Centrality = new TH1F("fh1Centrality","fh1Centrality; Centrality %",100,0,100);
938   fHistList->Add(fh1Centrality);
939
940   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
941   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
942   fHistList->Add(fh1Xsec);
943
944   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
945   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
946   fHistList->Add(fh1Trials);
947
948   fh1PtHard       = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
949   fHistList->Add(fh1PtHard);
950   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
951   fHistList->Add(fh1PtHardTrials);
952
953   Int_t fgkNPtBins = 100;
954   Float_t kMinPt   = 0.;
955   Float_t kMaxPt   = 100.;
956   Double_t *binsPt = new Double_t[fgkNPtBins+1];
957   for(Int_t i=0; i<=fgkNPtBins; i++) binsPt[i]=(Double_t)kMinPt + (kMaxPt-kMinPt)/fgkNPtBins*(Double_t)i ;
958
959   Int_t fgkNRel1PtUncertaintyBins=50;
960   Float_t fgkRel1PtUncertaintyMin = 0.;
961   Float_t fgkRel1PtUncertaintyMax = 1.;
962
963   Double_t *binsRel1PtUncertainty=new Double_t[fgkNRel1PtUncertaintyBins+1];
964   for(Int_t i=0; i<=fgkNRel1PtUncertaintyBins; i++) binsRel1PtUncertainty[i]=(Double_t)fgkRel1PtUncertaintyMin + (fgkRel1PtUncertaintyMax-fgkRel1PtUncertaintyMin)/fgkNRel1PtUncertaintyBins*(Double_t)i ;
965
966   fPtRelUncertainty1PtPrim = new TH2F("fPtRelUncertainty1PtPrim","fPtRelUncertainty1PtPrim",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
967   fHistList->Add(fPtRelUncertainty1PtPrim);
968
969   fPtRelUncertainty1PtSec = new TH2F("fPtRelUncertainty1PtSec","fPtRelUncertainty1PtSec",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
970   fHistList->Add(fPtRelUncertainty1PtSec);
971
972   TH1::AddDirectory(oldStatus);   
973
974   OpenFile(1);
975   OpenFile(2);
976
977   PostData(0,fHistList);
978   PostData(1,fCFManagerPos->GetParticleContainer());
979   PostData(2,fCFManagerNeg->GetParticleContainer());
980
981   if(binsPt)                delete [] binsPt;
982   if(binsRel1PtUncertainty) delete [] binsRel1PtUncertainty;
983
984 }
985
986 //________________________________________________________________________
987 Bool_t AliPWG4HighPtSpectra::IsHIJINGParticle(Int_t label)
988 {
989   //
990   // Return kTRUE in case particle is from HIJING event
991   //
992
993   AliGenHijingEventHeader* hijingHeader;
994   if(fReadAODData)
995     hijingHeader = GetHijingEventHeaderAOD(fAOD);
996   else
997     hijingHeader = GetHijingEventHeader(fMC);
998
999   Int_t nproduced = hijingHeader->NProduced();
1000
1001   if(fReadAODData) {
1002     AliAODMCParticle* mom = (AliAODMCParticle*) fArrayMCAOD->At(label);
1003     Int_t iMom = label;
1004     Int_t iParent = mom->GetMother();
1005     while(iParent!=-1){
1006       iMom = iParent;
1007       mom = (AliAODMCParticle*) fArrayMCAOD->At(iMom);
1008       iParent = mom->GetMother();
1009     }
1010
1011     if(iMom<nproduced) {
1012       return kTRUE;
1013     } else {
1014       return kFALSE;
1015     }
1016   }
1017   else { //ESD analysis
1018     TParticle * mom = fStack->Particle(label);
1019     Int_t iMom = label;
1020     Int_t iParent = mom->GetFirstMother();
1021     while(iParent!=-1){
1022       iMom = iParent;
1023       mom = fStack->Particle(iMom);
1024       iParent = mom->GetFirstMother();
1025     }
1026   
1027     if(iMom<nproduced) {
1028       return kTRUE;
1029     } else {
1030       return kFALSE;
1031     }
1032   }
1033 }
1034
1035 #endif