Add sigma2 jet shape and fill constituent info. for subtracted 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 = fAOD->GetTrack(iTrack);
423       if(!aodtrack)
424         continue;
425       if((!aodtrack->TestFilterBit(fFilterMask)) || 
426          ((!fCFManagerPos->CheckParticleCuts(kStepReconstructed,aodtrack)) && (aodtrack->Charge()>0.)) ||
427          ((!fCFManagerNeg->CheckParticleCuts(kStepReconstructed,aodtrack)) && (aodtrack->Charge()<0))    )
428         continue;
429       else {
430         track = static_cast<AliVTrack*>(aodtrack);
431       }
432       //fill the container
433       containerInputRec[0] = track->Pt();
434       containerInputRec[1] = track->Phi();
435       containerInputRec[2] = track->Eta();
436       containerInputRec[3] = track->GetTPCNcls();
437
438       if(track->Charge()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
439       if(track->Charge()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
440
441       if(fArrayMCAOD) {
442         Int_t label = TMath::Abs(track->GetLabel());
443         if(label>fArrayMCAOD->GetEntries()) {
444           continue;
445         }
446         AliAODMCParticle *particle = (AliAODMCParticle*) fArrayMCAOD->At(label);
447         if(!particle) {
448           continue;
449         }
450         //Only select particles generated by HIJING if requested
451         if(fbSelectHIJING) {
452           if(!IsHIJINGParticle(label)) {
453             continue;
454           }
455         }
456         containerInputRecMC[0] = particle->Pt();      
457         containerInputRecMC[1] = particle->Phi();      
458         containerInputRecMC[2] = particle->Eta();  
459         containerInputRecMC[3] = track->GetTPCNcls();
460
461         //Container with primaries
462         if(particle->IsPhysicalPrimary()) {
463           if(particle->Charge()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
464           if(particle->Charge()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
465           //Fill pT resolution plots for primaries
466           //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
467         }
468
469         //Container with secondaries
470         if (!particle->IsPhysicalPrimary() ) {
471           if(particle->Charge()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRecMC,kStepSecondaries);
472           if(particle->Charge()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRecMC,kStepSecondaries);
473           //Fill pT resolution plots for primaries
474           //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
475         }
476       }
477     }//track loop
478
479     //Fill MC containers if particles are findable
480     if(fArrayMCAOD) {
481       int noPart = fArrayMCAOD->GetEntriesFast();
482
483       for(int iPart = 1; iPart<noPart; iPart++) {
484         AliAODMCParticle *mcPart = (AliAODMCParticle*) fArrayMCAOD->At(iPart);
485         if(!mcPart) continue;
486
487         //Only select particles generated by HIJING if requested
488         if(fbSelectHIJING) {
489           if(!IsHIJINGParticle(iPart))
490             continue;
491         }
492
493         Int_t pdg = mcPart->PdgCode();
494
495         // select charged pions, protons, kaons , electrons, muons
496         if(TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) ==
497            321 || TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13){
498   
499           //fill the container
500           containerInputMC[0] = mcPart->Pt();
501           containerInputMC[1] = mcPart->Phi();      
502           containerInputMC[2] = mcPart->Eta();  
503           //      AliESDtrack *esdtrack = fESD->GetTrack(mcPart->GetLabel());
504           containerInputMC[3] = 159.;
505           
506           if(mcPart->IsPhysicalPrimary()) {
507             if(mcPart->Charge()>0. && fCFManagerPos->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
508             if(mcPart->Charge()<0. && fCFManagerNeg->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
509           }
510         }
511       }
512     }
513
514   }//end of AOD analysis
515   else {//ESD analysis
516     for (Int_t iTrack = 0; iTrack<nTracks; iTrack++) 
517     {
518       //Get track for analysis
519       AliESDtrack *track = 0x0;
520       AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
521       if(!esdtrack)
522         continue;
523       
524
525       if(fTrackType==4) {
526         if (!(fTrackCuts->AcceptTrack(esdtrack)))
527           continue;
528       }
529
530       if(fTrackType==1)
531         track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
532       else if(fTrackType==2 || fTrackType==4) {
533         track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(fESD),esdtrack->GetID());
534         if(!track)
535           continue;
536         
537         AliExternalTrackParam exParam;
538         Bool_t relate = track->RelateToVertexTPC(fESD->GetPrimaryVertexSPD(),fESD->GetMagneticField(),kVeryBig,&exParam);
539         if( !relate ) {
540           if(track) delete track;
541           continue;
542         }
543         track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
544       }
545       else if(fTrackType==5 || fTrackType==6) {
546         if(fTrackCuts->AcceptTrack(esdtrack)) {
547           continue;
548         }
549         else {
550           if( !(fTrackCutsReject->AcceptTrack(esdtrack)) && fTrackCuts->AcceptTrack(esdtrack) ) {
551
552             if(fTrackType==5) {
553               //use TPConly constrained track
554               track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
555               if(!track)
556                 continue;
557               
558               AliExternalTrackParam exParam;
559               Bool_t relate = track->RelateToVertexTPC(fESD->GetPrimaryVertexSPD(),fESD->GetMagneticField(),kVeryBig,&exParam);
560               if( !relate ) {
561                 if(track) delete track;
562                 continue;
563               }
564               track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
565             }
566             else if(fTrackType==6) {
567               //use global constrained track
568               track = new AliESDtrack(*esdtrack);
569               track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
570             }
571           }
572         }
573       }
574       else if(fTrackType==7) {
575         //use global constrained track
576         track = new AliESDtrack(*esdtrack);
577       }
578       else
579         track = esdtrack;
580     
581       if(!track)
582         continue;
583       
584
585       if(fTrackType==2 || fTrackType==4 || fTrackType==5) {
586         //Cut on chi2 of constrained fit
587         if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax && fSigmaConstrainedMax>0.) {
588           if(track) delete track;
589           continue;
590         }
591       }
592       if (!(fTrackCuts->AcceptTrack(track)) && fTrackType!=4 && fTrackType!=5 && fTrackType!=6) {
593         if(fTrackType==1 || fTrackType==2 || fTrackType==7) {
594           if(track) delete track;
595         }
596         continue;
597       }
598
599       if(fTrackType==7) {
600         if(fTrackCutsReject ) {
601           if(fTrackCutsReject->AcceptTrack(track) ) {
602             if(track) delete track;
603             continue;
604           }
605         }
606         if(esdtrack->GetConstrainedParam()) 
607           track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
608       }
609
610       if(!track) {
611         if(fTrackType==1 || fTrackType==2 || fTrackType==4 || fTrackType==5 || fTrackType==6 || fTrackType==7) {
612           if(track) delete track;
613         }
614         continue;
615       }
616
617       //fill the container
618       containerInputRec[0] = track->Pt();
619       containerInputRec[1] = track->Phi();
620       containerInputRec[2] = track->Eta();
621       containerInputRec[3] = track->GetTPCNclsIter1();
622
623       if(track->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
624       if(track->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
625
626         
627       //Only fill the MC containers if MC information is available
628       if(fMC) {
629         Int_t label = TMath::Abs(track->GetLabel());
630         if(label>fStack->GetNtrack()) {
631           if(fTrackType==1 || fTrackType==2 || fTrackType==7)
632             delete track;
633           continue;
634         }
635         //Only select particles generated by HIJING if requested
636         if(fbSelectHIJING) {
637           if(!IsHIJINGParticle(label)) {
638             if(fTrackType==1 || fTrackType==2 || fTrackType==7)
639               delete track;
640             continue;
641           }
642         }
643         TParticle *particle = fStack->Particle(label) ;
644         if(!particle) {
645           if(fTrackType==1 || fTrackType==2 || fTrackType==7)
646             delete track;
647           continue;
648         }
649         containerInputRecMC[0] = particle->Pt();      
650         containerInputRecMC[1] = particle->Phi();      
651         containerInputRecMC[2] = particle->Eta();  
652         containerInputRecMC[3] = track->GetTPCNclsIter1();
653
654         //Container with primaries
655         if(fStack->IsPhysicalPrimary(label)) {
656           if(particle->GetPDG()->Charge()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
657           if(particle->GetPDG()->Charge()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
658           //Fill pT resolution plots for primaries
659           fPtRelUncertainty1PtPrim->Fill(containerInputRec[0],containerInputRec[0]*TMath::Sqrt(track->GetSigma1Pt2()));
660         }
661
662         //Container with secondaries
663         if (!fStack->IsPhysicalPrimary(label) ) {
664           if(particle->GetPDG()->Charge()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepSecondaries);
665           if(particle->GetPDG()->Charge()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepSecondaries);
666           //Fill pT resolution plots for primaries
667           fPtRelUncertainty1PtSec->Fill(containerInputRec[0],containerInputRec[0]*TMath::Sqrt(track->GetSigma1Pt2()));
668         }
669       }
670
671       if(fTrackType==1  || fTrackType==2 || fTrackType==4 || fTrackType==5 || fTrackType==6 || fTrackType==7) {
672         if(track) delete track;
673       }
674     }//track loop
675
676     //Fill MC containers if particles are findable
677     if(fMC) {
678       for(int iPart = 1; iPart<(fMC->GetNumberOfPrimaries()); iPart++) {
679         AliMCParticle *mcPart  = (AliMCParticle*)fMC->GetTrack(iPart);
680         if(!mcPart) continue;
681   
682         //Only select particles generated by HIJING if requested
683         if(fbSelectHIJING) {
684           if(!IsHIJINGParticle(iPart))
685             continue;
686         }
687   
688         Int_t pdg = mcPart->PdgCode();
689         
690         // select charged pions, protons, kaons , electrons, muons
691         if(TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) ==
692            321 || TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13){
693   
694           //fill the container
695           containerInputMC[0] = mcPart->Pt();
696           containerInputMC[1] = mcPart->Phi();      
697           containerInputMC[2] = mcPart->Eta();  
698           //      AliESDtrack *esdtrack = fESD->GetTrack(mcPart->GetLabel());
699           containerInputMC[3] = 159.;
700           
701           if(fStack->IsPhysicalPrimary(iPart)) {
702             if(mcPart->Charge()>0. && fCFManagerPos->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
703             if(mcPart->Charge()<0. && fCFManagerNeg->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
704           }
705         }
706       }
707     }
708   }//end of ESD analysis
709
710   PostData(0,fHistList);
711   PostData(1,fCFManagerPos->GetParticleContainer());
712   PostData(2,fCFManagerNeg->GetParticleContainer());
713   
714 }
715
716 //________________________________________________________________________
717 Bool_t AliPWG4HighPtSpectra::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials)
718 {
719   //
720   // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
721   // This is to called in Notify and should provide the path to the AOD/ESD file
722   // Copied from AliAnalysisTaskJetSpectrum2
723   //
724
725   TString file(currFile);  
726   fXsec = 0;
727   fTrials = 1;
728
729   if(file.Contains("root_archive.zip#")){
730     Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
731     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
732     file.Replace(pos+1,20,"");
733   }
734   else {
735     // not an archive take the basename....
736     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
737   }
738   //  Printf("%s",file.Data());
739    
740   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
741   if(!fxsec){
742     // next trial fetch the histgram file
743     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
744     if(!fxsec){
745       // not a severe condition but indicates that we have no information
746       return kFALSE;
747     }
748     else{
749       // find the tlist we want to be independtent of the name so use the Tkey
750       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0); 
751       if(!key){
752         fxsec->Close();
753         return kFALSE;
754       }
755       TList *list = dynamic_cast<TList*>(key->ReadObj());
756       if(!list){
757         fxsec->Close();
758         return kFALSE;
759       }
760       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
761       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
762       fxsec->Close();
763     }
764   } // no tree pyxsec.root
765   else {
766     TTree *xtree = (TTree*)fxsec->Get("Xsection");
767     if(!xtree){
768       fxsec->Close();
769       return kFALSE;
770     }
771     UInt_t   ntrials  = 0;
772     Double_t  xsection  = 0;
773     xtree->SetBranchAddress("xsection",&xsection);
774     xtree->SetBranchAddress("ntrials",&ntrials);
775     xtree->GetEntry(0);
776     fTrials = ntrials;
777     fXsec = xsection;
778     fxsec->Close();
779   }
780   return kTRUE;
781 }
782
783 //________________________________________________________________________
784 Bool_t AliPWG4HighPtSpectra::Notify()
785 {
786   //
787   // Implemented Notify() to read the cross sections
788   // and number of trials from pyxsec.root
789   // Copied from AliAnalysisTaskJetSpectrum2
790   // 
791
792   if(fNoPythiaInfo)
793     return kTRUE;
794
795   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
796   Float_t xsection = 0;
797   Float_t ftrials  = 1;
798
799   if(tree){
800     TFile *curfile = tree->GetCurrentFile();
801     if (!curfile) {
802       Error("Notify","No current file");
803       return kFALSE;
804     }
805     if(!fh1Xsec||!fh1Trials){
806       //      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
807       return kFALSE;
808     }
809     PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
810     fh1Xsec->Fill("<#sigma>",xsection);
811     fh1Trials->Fill("#sum{ntrials}",ftrials);
812   } 
813   return kTRUE;
814 }
815
816 //________________________________________________________________________
817 AliGenPythiaEventHeader*  AliPWG4HighPtSpectra::GetPythiaEventHeader(AliMCEvent *mcEvent)
818 {
819   
820   if(!mcEvent)return 0;
821   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
822   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
823   if(!pythiaGenHeader){
824     // cocktail ??
825     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
826     
827     if (!genCocktailHeader) {
828       //      AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
829       //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
830       return 0;
831     }
832     TList* headerList = genCocktailHeader->GetHeaders();
833     for (Int_t i=0; i<headerList->GetEntries(); i++) {
834       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
835       if (pythiaGenHeader)
836         break;
837     }
838     if(!pythiaGenHeader){
839       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
840       return 0;
841     }
842   }
843   return pythiaGenHeader;
844
845 }
846
847 //________________________________________________________________________
848 AliGenHijingEventHeader*  AliPWG4HighPtSpectra::GetHijingEventHeader(AliMCEvent *mcEvent)
849 {
850   
851   if(!mcEvent)return 0;
852   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
853   AliGenHijingEventHeader* hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
854   if(!hijingGenHeader){
855     // cocktail ??
856     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
857     
858     if (!genCocktailHeader) {
859       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
860       //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
861       return 0;
862     }
863     TList* headerList = genCocktailHeader->GetHeaders();
864     for (Int_t i=0; i<headerList->GetEntries(); i++) {
865       hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(headerList->At(i));
866       if (hijingGenHeader)
867         break;
868     }
869     if(!hijingGenHeader){
870       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Hijing event header not found");
871       return 0;
872     }
873   }
874   return hijingGenHeader;
875
876 }
877
878 //________________________________________________________________________
879 AliGenHijingEventHeader*  AliPWG4HighPtSpectra::GetHijingEventHeaderAOD(AliAODEvent *aodEvent)
880 {
881   AliGenHijingEventHeader* hijingGenHeader = 0x0;
882   if(aodEvent) {
883     AliAODMCHeader* mcHeader = (AliAODMCHeader*) aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
884     TList* headerList = mcHeader->GetCocktailHeaders();
885     for (Int_t i=0; i<headerList->GetEntries(); i++) {
886       hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(headerList->At(i));
887       if (hijingGenHeader)
888         break;
889     }
890   }
891   return hijingGenHeader;
892 }
893
894 //___________________________________________________________________________
895 void AliPWG4HighPtSpectra::Terminate(Option_t*)
896 {
897   // The Terminate() function is the last function to be called during
898   // a query. It always runs on the client, it can be used to present
899   // the results graphically or save the results to file.
900
901 }
902
903 //___________________________________________________________________________
904 void AliPWG4HighPtSpectra::CreateOutputObjects() 
905 {
906   //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
907   //TO BE SET BEFORE THE EXECUTION OF THE TASK
908   //
909   AliDebug(2,Form("CreateOutputObjects CreateOutputObjects of task %s", GetName()));
910
911   Bool_t oldStatus = TH1::AddDirectoryStatus();
912   TH1::AddDirectory(kFALSE); 
913
914   //slot #1
915   OpenFile(0);
916   fHistList = new TList();
917   fHistList->SetOwner(kTRUE);
918   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
919   fHistList->Add(fNEventAll);
920   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
921   fHistList->Add(fNEventSel);
922
923   fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
924   //Set labels
925   fNEventReject->Fill("noESD/AOD",0);
926   fNEventReject->Fill("Trigger",0);
927   fNEventReject->Fill("NTracks<2",0);
928   fNEventReject->Fill("noVTX",0);
929   fNEventReject->Fill("VtxStatus",0);
930   fNEventReject->Fill("NCont<2",0);
931   fNEventReject->Fill("ZVTX>10",0);
932   fNEventReject->Fill("cent",0);
933   fNEventReject->Fill("cent>90",0);
934   fHistList->Add(fNEventReject);
935
936   fh1Centrality = new TH1F("fh1Centrality","fh1Centrality; Centrality %",100,0,100);
937   fHistList->Add(fh1Centrality);
938
939   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
940   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
941   fHistList->Add(fh1Xsec);
942
943   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
944   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
945   fHistList->Add(fh1Trials);
946
947   fh1PtHard       = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
948   fHistList->Add(fh1PtHard);
949   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
950   fHistList->Add(fh1PtHardTrials);
951
952   Int_t fgkNPtBins = 100;
953   Float_t kMinPt   = 0.;
954   Float_t kMaxPt   = 100.;
955   Double_t *binsPt = new Double_t[fgkNPtBins+1];
956   for(Int_t i=0; i<=fgkNPtBins; i++) binsPt[i]=(Double_t)kMinPt + (kMaxPt-kMinPt)/fgkNPtBins*(Double_t)i ;
957
958   Int_t fgkNRel1PtUncertaintyBins=50;
959   Float_t fgkRel1PtUncertaintyMin = 0.;
960   Float_t fgkRel1PtUncertaintyMax = 1.;
961
962   Double_t *binsRel1PtUncertainty=new Double_t[fgkNRel1PtUncertaintyBins+1];
963   for(Int_t i=0; i<=fgkNRel1PtUncertaintyBins; i++) binsRel1PtUncertainty[i]=(Double_t)fgkRel1PtUncertaintyMin + (fgkRel1PtUncertaintyMax-fgkRel1PtUncertaintyMin)/fgkNRel1PtUncertaintyBins*(Double_t)i ;
964
965   fPtRelUncertainty1PtPrim = new TH2F("fPtRelUncertainty1PtPrim","fPtRelUncertainty1PtPrim",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
966   fHistList->Add(fPtRelUncertainty1PtPrim);
967
968   fPtRelUncertainty1PtSec = new TH2F("fPtRelUncertainty1PtSec","fPtRelUncertainty1PtSec",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
969   fHistList->Add(fPtRelUncertainty1PtSec);
970
971   TH1::AddDirectory(oldStatus);   
972
973   OpenFile(1);
974   OpenFile(2);
975
976   PostData(0,fHistList);
977   PostData(1,fCFManagerPos->GetParticleContainer());
978   PostData(2,fCFManagerNeg->GetParticleContainer());
979
980   if(binsPt)                delete [] binsPt;
981   if(binsRel1PtUncertainty) delete [] binsRel1PtUncertainty;
982
983 }
984
985 //________________________________________________________________________
986 Bool_t AliPWG4HighPtSpectra::IsHIJINGParticle(Int_t label)
987 {
988   //
989   // Return kTRUE in case particle is from HIJING event
990   //
991
992   AliGenHijingEventHeader* hijingHeader;
993   if(fReadAODData)
994     hijingHeader = GetHijingEventHeaderAOD(fAOD);
995   else
996     hijingHeader = GetHijingEventHeader(fMC);
997
998   Int_t nproduced = hijingHeader->NProduced();
999
1000   if(fReadAODData) {
1001     AliAODMCParticle* mom = (AliAODMCParticle*) fArrayMCAOD->At(label);
1002     Int_t iMom = label;
1003     Int_t iParent = mom->GetMother();
1004     while(iParent!=-1){
1005       iMom = iParent;
1006       mom = (AliAODMCParticle*) fArrayMCAOD->At(iMom);
1007       iParent = mom->GetMother();
1008     }
1009
1010     if(iMom<nproduced) {
1011       return kTRUE;
1012     } else {
1013       return kFALSE;
1014     }
1015   }
1016   else { //ESD analysis
1017     TParticle * mom = fStack->Particle(label);
1018     Int_t iMom = label;
1019     Int_t iParent = mom->GetFirstMother();
1020     while(iParent!=-1){
1021       iMom = iParent;
1022       mom = fStack->Particle(iMom);
1023       iParent = mom->GetFirstMother();
1024     }
1025   
1026     if(iMom<nproduced) {
1027       return kTRUE;
1028     } else {
1029       return kFALSE;
1030     }
1031   }
1032 }
1033
1034 #endif