Compatibility with ROOT trunk
[u/mrichter/AliRoot.git] / ANALYSIS / AliEPSelectionTask.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2008, 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 //   Class AliEPSelectionTask
18 //   Class to determine event plane            
19 //   author: Alberica Toia, Johanna Gramling
20 //*****************************************************
21
22 #include "AliEPSelectionTask.h"
23
24 #include <TTree.h>
25 #include <TList.h>
26 #include <TH1F.h>
27 #include <TH2F.h>
28 #include <THnSparse.h>
29 #include <TProfile.h>
30 #include <TFile.h>
31 #include <TObjString.h>
32 #include <TString.h>
33 #include <TCanvas.h>
34 #include <TROOT.h>
35 #include <TDirectory.h>
36 #include <TSystem.h>
37 #include <iostream>
38 #include <TRandom2.h>
39 #include <TArrayF.h>
40
41 #include "AliAnalysisManager.h"
42 #include "AliVEvent.h"
43 #include "AliESD.h"
44 #include "AliESDEvent.h"
45 #include "AliMCEvent.h"
46 #include "AliESDtrack.h"
47 #include "AliESDtrackCuts.h"
48 #include "AliESDHeader.h"
49 #include "AliESDInputHandler.h"
50 #include "AliAODHandler.h"
51 #include "AliAODEvent.h"
52 #include "AliHeader.h"
53 #include "AliGenHijingEventHeader.h"
54 #include "AliAnalysisTaskSE.h"
55 #include "AliPhysicsSelectionTask.h"
56 #include "AliPhysicsSelection.h"
57 #include "AliBackgroundSelection.h"
58 #include "AliESDUtils.h"
59 #include "AliOADBContainer.h"
60 #include "AliAODMCHeader.h"
61 #include "AliAODTrack.h"
62 #include "AliVTrack.h"
63 #include "AliEventplane.h"
64
65 using std::cout;
66 using std::endl;
67 ClassImp(AliEPSelectionTask)
68
69 //________________________________________________________________________
70 AliEPSelectionTask::AliEPSelectionTask():
71 AliAnalysisTaskSE(),
72   fAnalysisInput("ESD"),
73   fTrackType("TPC"),
74   fPeriod(""),
75   fUseMCRP(kFALSE),
76   fUsePhiWeight(kFALSE),
77   fUsePtWeight(kFALSE),
78   fSaveTrackContribution(kFALSE),
79   fUserphidist(kFALSE),
80   fUsercuts(kFALSE),
81   fRunNumber(-15),
82   fAODfilterbit(1),
83   fEtaGap(0.),
84   fSplitMethod(0),
85   fESDtrackCuts(0),
86   fEPContainer(0),
87   fSparseDist(0),
88   fHruns(0),
89   fQVector(0),
90   fQContributionX(0),
91   fQContributionY(0),
92   fEventplaneQ(0),
93   fQsub1(0),
94   fQsub2(0),
95   fQsubRes(0),  
96   fOutputList(0),
97   fHOutEventplaneQ(0),
98   fHOutPhi(0),
99   fHOutPhiCorr(0),
100   fHOutsub1sub2(0),
101   fHOutNTEPRes(0),
102   fHOutPTPsi(0),
103   fHOutDiff(0),
104   fHOutleadPTPsi(0)
105 {   
106   // Default constructor
107   AliInfo("Event Plane Selection enabled.");
108   for(Int_t i = 0; i < 4; ++i) {
109      fPhiDist[i] = 0;
110   }
111 }   
112
113 //________________________________________________________________________
114 AliEPSelectionTask::AliEPSelectionTask(const char *name):
115   AliAnalysisTaskSE(name),
116   fAnalysisInput("ESD"),
117   fTrackType("TPC"),
118   fPeriod(""),
119   fUseMCRP(kFALSE),
120   fUsePhiWeight(kFALSE),
121   fUsePtWeight(kFALSE),  
122   fSaveTrackContribution(kFALSE),
123   fUserphidist(kFALSE),
124   fUsercuts(kFALSE),
125   fRunNumber(-15),
126   fAODfilterbit(1),
127   fEtaGap(0.),
128   fSplitMethod(0),
129   fESDtrackCuts(0),
130   fEPContainer(0),
131   fSparseDist(0),
132   fHruns(0),
133   fQVector(0),
134   fQContributionX(0),
135   fQContributionY(0),
136   fEventplaneQ(0),
137   fQsub1(0),
138   fQsub2(0),
139   fQsubRes(0),
140   fOutputList(0),
141   fHOutEventplaneQ(0),
142   fHOutPhi(0),
143   fHOutPhiCorr(0),
144   fHOutsub1sub2(0),
145   fHOutNTEPRes(0),
146   fHOutPTPsi(0),
147   fHOutDiff(0),
148   fHOutleadPTPsi(0)
149 {
150   // Default constructor
151   AliInfo("Event Plane Selection enabled.");
152   DefineOutput(1, TList::Class());
153   for(Int_t i = 0; i < 4; i++) {
154      fPhiDist[i] = 0;
155   }
156 }
157  
158 //________________________________________________________________________
159 AliEPSelectionTask::~AliEPSelectionTask()
160 {
161   // Destructor  
162   if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()){
163       delete fOutputList;
164       fOutputList = 0;
165   }
166   if (fESDtrackCuts){
167       delete fESDtrackCuts;
168       fESDtrackCuts = 0;
169   }
170   if (fUserphidist) {
171     if (fPhiDist[0]) {
172       delete fPhiDist[0];
173       fPhiDist[0] = 0;
174     }
175   }
176   if (fEPContainer){
177       delete fEPContainer;
178       fEPContainer = 0;
179   }
180   if (fPeriod.CompareTo("LHC11h")==0){
181       for(Int_t i = 0; i < 4; i++) {
182         if(fPhiDist[i]){
183           delete fPhiDist[i];
184           fPhiDist[i] = 0;
185         }
186       }
187       if(fHruns) delete fHruns;
188   }
189 }  
190
191 //________________________________________________________________________
192 void AliEPSelectionTask::UserCreateOutputObjects()
193 {  
194   // Create the output containers
195   if (fDebug>1) printf("AliEPSelectionTask::UserCreateOutputObjects() \n");
196   AliLog::SetClassDebugLevel("AliEPSelectionTask", AliLog::kInfo);
197
198   fOutputList = new TList();
199   fOutputList->SetOwner();
200   fHOutEventplaneQ = new TH1F("fHOutEventplaneQ","fHOutEventplaneQ; Event Plane Q",100,0,TMath::Pi());
201   fHOutPhi = new TH1F("fHOutPhi","fHOutPhi; Phi Distribution",100,0,TMath::TwoPi());
202   fHOutPhiCorr = new TH1F("fHOutPhiCorr","fHOutPhiCorr; Corrected Phi Distribution",100,0,TMath::TwoPi());
203   fHOutsub1sub2 = new TH2F("fHOutsub1sub2","fHOutsub1sub2; EP1; EP2",100,0,TMath::Pi(),100,0,TMath::Pi());
204   fHOutNTEPRes = new TH2F("fHOutNTEPRes","fHOutNTEPRes; Number of Tracks; Event Plane Resolution",100,0,5000,100,-TMath::Pi(),TMath::Pi());
205   fHOutPTPsi = new TH2F("fHOutPTPsi","fHOutPTPsi; PT; Phi-EP",100,0,20,100,0,TMath::Pi());
206   fHOutDiff = new TH2F("fHOutDiff","fHOutDiff; EP; MCEP",100,0,TMath::Pi(),100,0,TMath::Pi());
207   fHOutleadPTPsi = new TH2F("fHOutleadPTPsi","fHOutleadPTPsi; leadPT; EP",100,0,TMath::Pi(),100,0,TMath::Pi());
208
209   fOutputList->Add(fHOutEventplaneQ);
210   fOutputList->Add(fHOutPhi);
211   fOutputList->Add(fHOutPhiCorr);
212   fOutputList->Add(fHOutsub1sub2);
213   fOutputList->Add(fHOutNTEPRes);
214   fOutputList->Add(fHOutPTPsi);
215   fOutputList->Add(fHOutleadPTPsi);
216   fOutputList->Add(fHOutDiff);
217   
218   PostData(1, fOutputList); 
219
220
221 }
222
223 //________________________________________________________________________
224 void AliEPSelectionTask::UserExec(Option_t */*option*/)
225
226   // Execute analysis for current event:
227   if (fDebug>1) printf(" **** AliEPSelectionTask::UserExec() \n");
228   
229 //   fRunNumber = -15;
230  
231   AliEventplane *esdEP;
232   TVector2 qq1;
233   TVector2 qq2;
234   Double_t fRP = 0.; // monte carlo reaction plane angle
235     
236   if (fAnalysisInput.CompareTo("ESD")==0){
237
238     AliVEvent* event = InputEvent();
239     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
240     if (esd){    
241       if (!(fRunNumber == esd->GetRunNumber())) {
242           fRunNumber = esd->GetRunNumber();
243           AliInfo(Form("Changing Phi-distribution to run %d",fRunNumber));
244           SetOADBandPeriod();
245           SetPhiDist();      
246       }
247       
248       
249       if (fUseMCRP) {
250           AliMCEvent* mcEvent  = MCEvent();      
251           if (mcEvent && mcEvent->GenEventHeader()) {
252               AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
253               if (headerH) fRP = headerH->ReactionPlaneAngle();
254           }
255       }
256       
257       esdEP = esd->GetEventplane();
258       if (fSaveTrackContribution) {
259         esdEP->GetQContributionXArray()->Set(esd->GetNumberOfTracks());
260         esdEP->GetQContributionYArray()->Set(esd->GetNumberOfTracks());
261         esdEP->GetQContributionXArraysub1()->Set(esd->GetNumberOfTracks());
262         esdEP->GetQContributionYArraysub1()->Set(esd->GetNumberOfTracks());
263         esdEP->GetQContributionXArraysub2()->Set(esd->GetNumberOfTracks());
264         esdEP->GetQContributionYArraysub2()->Set(esd->GetNumberOfTracks());
265       }
266       
267       TObjArray* tracklist = new TObjArray;
268       if (fTrackType.CompareTo("GLOBAL")==0) tracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE);
269       if (fTrackType.CompareTo("TPC")==0 && fPeriod.CompareTo("LHC10h")==0) tracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE);
270       else if (fTrackType.CompareTo("TPC")==0 && fPeriod.CompareTo("LHC11h")==0) tracklist = GetTracksForLHC11h(esd);
271       const int nt = tracklist->GetEntries();
272       
273       if (nt>4){
274         fQVector = new TVector2(GetQ(esdEP,tracklist));
275         fEventplaneQ = fQVector->Phi()/2; 
276         GetQsub(qq1, qq2, tracklist, esdEP);
277         fQsub1 = new TVector2(qq1);
278         fQsub2 = new TVector2(qq2);
279         fQsubRes = (fQsub1->Phi()/2 - fQsub2->Phi()/2);
280         
281         esdEP->SetQVector(fQVector);
282         esdEP->SetEventplaneQ(fEventplaneQ);
283         esdEP->SetQsub(fQsub1,fQsub2);
284         esdEP->SetQsubRes(fQsubRes);
285         
286         fHOutEventplaneQ->Fill(fEventplaneQ);
287         fHOutsub1sub2->Fill(fQsub1->Phi()/2,fQsub2->Phi()/2);
288         fHOutNTEPRes->Fill(nt,fQsubRes);
289
290         if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
291         
292         for (int iter = 0; iter<nt;iter++){
293           AliESDtrack* track = dynamic_cast<AliESDtrack*> (tracklist->At(iter));
294           if (track) {
295             float delta = track->Phi()-fEventplaneQ;
296             while (delta < 0) delta += TMath::Pi();
297             while (delta > TMath::Pi()) delta -= TMath::Pi();
298             fHOutPTPsi->Fill(track->Pt(),delta);
299             fHOutPhi->Fill(track->Phi());
300             fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track)); 
301           }
302         }
303         
304         AliESDtrack* trmax = esd->GetTrack(0);
305         for (int iter = 1; iter<nt;iter++){
306           AliESDtrack* track = dynamic_cast<AliESDtrack*> (tracklist->At(iter));
307           if (track && (track->Pt() > trmax->Pt())) trmax = track;
308         }
309         fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);      
310       }
311       tracklist->Clear();
312       delete tracklist;
313       tracklist = 0;
314     }
315   }
316   
317     else if (fAnalysisInput.CompareTo("AOD")==0){
318     AliVEvent* event = InputEvent();
319     AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
320
321     if (aod){
322       if (!(fRunNumber == aod->GetRunNumber())) {
323         fRunNumber = aod->GetRunNumber();
324         AliInfo(Form("Changing Phi-distribution to run %d",fRunNumber));
325         SetOADBandPeriod();
326         SetPhiDist();      
327       }
328
329       if (fUseMCRP) {
330         AliAODMCHeader *headerH = dynamic_cast<AliAODMCHeader*>(aod->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
331         if (headerH) fRP = headerH->GetReactionPlaneAngle();
332       }
333   
334       esdEP = aod->GetHeader()->GetEventplaneP();
335       esdEP->Reset(); 
336      
337       Int_t maxID = 0;
338       TObjArray* tracklist = GetAODTracksAndMaxID(aod,maxID);
339         
340       if (fSaveTrackContribution) {
341         esdEP->GetQContributionXArray()->Set(maxID+1);
342         esdEP->GetQContributionYArray()->Set(maxID+1);
343         esdEP->GetQContributionXArraysub1()->Set(maxID+1);
344         esdEP->GetQContributionYArraysub1()->Set(maxID+1);
345         esdEP->GetQContributionXArraysub2()->Set(maxID+1);
346         esdEP->GetQContributionYArraysub2()->Set(maxID+1);
347       }
348         
349       const int NT = tracklist->GetEntries();
350       
351       if (NT>4){
352         fQVector = new TVector2(GetQ(esdEP,tracklist));
353         fEventplaneQ = fQVector->Phi()/2.; 
354         GetQsub(qq1, qq2, tracklist, esdEP);
355         fQsub1 = new TVector2(qq1);
356         fQsub2 = new TVector2(qq2);
357         fQsubRes = (fQsub1->Phi()/2. - fQsub2->Phi()/2.);
358         
359         esdEP->SetQVector(fQVector);
360         esdEP->SetEventplaneQ(fEventplaneQ);
361         esdEP->SetQsub(fQsub1,fQsub2);
362         esdEP->SetQsubRes(fQsubRes);
363         
364         fHOutEventplaneQ->Fill(fEventplaneQ);
365         fHOutsub1sub2->Fill(fQsub1->Phi()/2.,fQsub2->Phi()/2.);
366         fHOutNTEPRes->Fill(NT,fQsubRes);
367         
368         if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
369         
370         for (int iter = 0; iter<NT;iter++){
371           AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
372           if (track) {
373             float delta = track->Phi()-fEventplaneQ;
374             while (delta < 0) delta += TMath::Pi();
375             while (delta > TMath::Pi()) delta -= TMath::Pi();
376             fHOutPTPsi->Fill(track->Pt(),delta);
377             fHOutPhi->Fill(track->Phi());
378             fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
379           }
380         }
381         
382         AliAODTrack* trmax = aod->GetTrack(0);
383         for (int iter = 1; iter<NT;iter++){
384           AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
385           if (track && (track->Pt() > trmax->Pt())) trmax = track;
386         }
387         fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);      
388       }     
389       delete tracklist;
390       tracklist = 0;
391     }   
392         
393     
394   }  
395
396   
397   else {
398     printf(" Analysis Input not known!\n\n ");
399     return;
400   }
401   PostData(1, fOutputList); 
402 }
403
404 //________________________________________________________________________
405 void AliEPSelectionTask::Terminate(Option_t */*option*/)
406 {
407   // Terminate analysis
408 }
409
410 //__________________________________________________________________________
411 TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
412 {
413 // Get the Q vector
414   TVector2 mQ;
415   float mQx=0, mQy=0;
416   AliVTrack* track;
417   Double_t weight;
418   Int_t idtemp = -1;
419   
420   int nt = tracklist->GetEntries();
421
422   for (int i=0; i<nt; i++){
423     weight = 1;
424     track = dynamic_cast<AliVTrack*> (tracklist->At(i));
425     if (track) {
426       weight = GetWeight(track);
427     if (fSaveTrackContribution){
428       idtemp = track->GetID(); 
429       if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
430       EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi()),idtemp);
431       EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi()),idtemp);
432      }
433      mQx += (weight*cos(2*track->Phi()));
434      mQy += (weight*sin(2*track->Phi()));
435     }
436   }
437   mQ.Set(mQx,mQy);
438   return mQ;
439 }
440   
441   //________________________________________________________________________
442 void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist,AliEventplane* EP)
443 {
444 // Get Qsub
445   TVector2 mQ[2];
446   float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
447   Double_t weight;
448
449   AliVTrack* track;
450   TRandom2 rn = 0;
451   
452   int nt = tracklist->GetEntries();
453   int trackcounter1=0, trackcounter2=0;
454   int idtemp = 0;
455
456   if (fSplitMethod == AliEPSelectionTask::kRandom){
457     
458     for (Int_t i = 0; i < nt; i++) {
459       weight = 1;
460       track = dynamic_cast<AliVTrack*> (tracklist->At(i));
461       if (!track) continue;
462       weight = GetWeight(track);
463       idtemp = track->GetID(); 
464       if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
465     
466       // This loop splits the track set into 2 random subsets
467       if( trackcounter1 < int(nt/2.) && trackcounter2 < int(nt/2.)){
468         float random = rn.Rndm();
469         if(random < .5){
470           mQx1 += (weight*cos(2*track->Phi()));
471           mQy1 += (weight*sin(2*track->Phi()));
472           if (fSaveTrackContribution){
473             EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
474             EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
475           }
476           trackcounter1++;
477         }
478         else {
479           mQx2 += (weight*cos(2*track->Phi()));
480           mQy2 += (weight*sin(2*track->Phi()));
481           if (fSaveTrackContribution){
482             EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
483             EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
484           }
485           trackcounter2++;
486         }
487       }
488       else if( trackcounter1 >= int(nt/2.)){
489         mQx2 += (weight*cos(2*track->Phi()));
490         mQy2 += (weight*sin(2*track->Phi()));
491         if (fSaveTrackContribution){
492           EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
493           EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
494         }
495         trackcounter2++;
496       }
497       else {
498         mQx1 += (weight*cos(2*track->Phi()));
499         mQy1 += (weight*sin(2*track->Phi()));
500         if (fSaveTrackContribution){
501           EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
502           EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
503         }
504         trackcounter1++;
505       }
506     }
507   } else if (fSplitMethod == AliEPSelectionTask::kEta) {
508      
509     for (Int_t i = 0; i < nt; i++) {
510       weight = 1;
511       track = dynamic_cast<AliVTrack*> (tracklist->At(i));
512       if (!track) continue;
513       weight = GetWeight(track);
514       Double_t eta = track->Eta();
515       idtemp = track->GetID(); 
516       if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
517
518       if (eta > fEtaGap/2.) {  
519         mQx1 += (weight*cos(2*track->Phi()));
520         mQy1 += (weight*sin(2*track->Phi()));
521         if (fSaveTrackContribution){
522           EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
523           EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
524         }
525       } else if (eta < -1.*fEtaGap/2.) {
526         mQx2 += (weight*cos(2*track->Phi()));
527         mQy2 += (weight*sin(2*track->Phi()));
528         if (fSaveTrackContribution){
529           EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
530           EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
531         }
532       }
533     }
534   } else {
535     printf("plane resolution determination method not available!\n\n ");
536     return;
537   }
538      
539   mQ[0].Set(mQx1,mQy1);
540   mQ[1].Set(mQx2,mQy2);
541   Q1 = mQ[0];
542   Q2 = mQ[1];
543 }
544
545 //________________________________________________________________________
546 void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
547   
548   if(fESDtrackCuts){ 
549     delete fESDtrackCuts;
550     fESDtrackCuts = 0;
551   }
552   if (fAnalysisInput.CompareTo("AOD")==0){
553     AliInfo("ESD track cuts not possible for AOD analysis; please use SetPersonalAODtrackCuts(); using TPC only track cuts");  
554     fUsercuts = kFALSE;
555     SetTrackType("TPC");
556     return;
557   } 
558   fUsercuts = kTRUE;
559   fESDtrackCuts = trackcuts;
560 }
561
562 //________________________________________________________________________
563 void AliEPSelectionTask::SetPersonalAODtrackCuts(UInt_t filterbit, Float_t etalow, Float_t etaup, Float_t ptlow, Float_t ptup){
564   
565   if(fESDtrackCuts){ 
566     delete fESDtrackCuts;
567     fESDtrackCuts = 0;
568   }
569   if (fAnalysisInput.CompareTo("ESD")==0){
570     AliInfo("AOD track cuts not possible for ESD analysis; please use SetPersonalESDtrackCuts(); using TPC only track cuts");  
571     fUsercuts = kFALSE;
572     SetTrackType("TPC");
573     return;
574   }
575   fUsercuts = kTRUE;
576   fESDtrackCuts = new AliESDtrackCuts();
577   fESDtrackCuts->SetPtRange(ptlow,ptup);
578   fESDtrackCuts->SetEtaRange(etalow,etaup);
579   fAODfilterbit = filterbit;
580 }
581
582 //_____________________________________________________________________________
583
584 void AliEPSelectionTask::SetTrackType(TString tracktype){
585   fTrackType = tracktype;
586   if (!fUsercuts) {
587   if (fTrackType.CompareTo("GLOBAL")==0){ 
588     fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
589     fAODfilterbit = 32;
590   }     
591   if (fTrackType.CompareTo("TPC")==0){  
592     fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
593     fAODfilterbit = 128;
594   }
595   fESDtrackCuts->SetPtRange(0.15,20.);
596   fESDtrackCuts->SetEtaRange(-0.8,0.8);
597   }
598 }
599
600 //________________________________________________________________________
601 Double_t AliEPSelectionTask::GetWeight(TObject* track1)
602 {
603   Double_t ptweight=1;
604   AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
605   if (fUsePtWeight && track) {      
606     if (track->Pt()<2) ptweight=track->Pt();
607     else ptweight=2;
608   }
609   return ptweight*GetPhiWeight(track);
610 }
611
612 //________________________________________________________________________
613 Double_t AliEPSelectionTask::GetPhiWeight(TObject* track1)
614 {
615   Double_t phiweight=1;
616   AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
617
618   TH1F *phiDist = 0x0;
619   if(track) phiDist = SelectPhiDist(track);
620   
621   if (fUsePhiWeight && phiDist && track) {
622     Double_t nParticles = phiDist->Integral();
623     Double_t nPhibins = phiDist->GetNbinsX();
624   
625     Double_t Phi = track->Phi();
626     
627     while (Phi<0) Phi += TMath::TwoPi();
628     while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
629       
630     Double_t PhiDistValue = phiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
631     
632     if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
633   }
634   return phiweight;
635 }
636
637 //__________________________________________________________________________
638 void AliEPSelectionTask::SetPhiDist() 
639 {
640   if(!fUserphidist && (fPeriod.CompareTo("LHC10h") == 0 || fPeriod.CompareTo("LHC11h") == 0)) { // if it's already set and custom class is required, we use the one provided by the user
641
642     if (fPeriod.CompareTo("LHC10h")==0)
643        {
644         fPhiDist[0] = (TH1F*) fEPContainer->GetObject(fRunNumber, "Default");}
645         else if(fPeriod.CompareTo("LHC11h")==0){
646             Int_t runbin=fHruns->FindBin(fRunNumber);
647             if (fHruns->GetBinContent(runbin) > 1){
648               fSparseDist->GetAxis(0)->SetRange(runbin,runbin);
649               }
650             else if(fHruns->GetBinContent(runbin) < 2){
651                fSparseDist->GetAxis(0)->SetRange(1,2901); // not calibrated run, use integrated phi-weights
652                AliInfo("Using integrated Phi-weights for this run");
653                }
654             for (Int_t i = 0; i<4 ;i++)
655             {
656              if(fPhiDist[i]){
657                delete fPhiDist[i];
658                fPhiDist[i] = 0x0;
659                }
660              if(i == 0){
661                fSparseDist->GetAxis(1)->SetRange(1,1);  // neg charge
662                fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta 
663              if(i == 1){
664                fSparseDist->GetAxis(1)->SetRange(2,2);  // pos charge
665                fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta
666              if(i == 2){
667                fSparseDist->GetAxis(1)->SetRange(1,1);  // neg charge
668                fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta 
669              if(i == 3){
670                fSparseDist->GetAxis(1)->SetRange(2,2);  // pos charge
671                fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta
672              fPhiDist[i] = (TH1F*)fSparseDist->Projection(3); // Projection on Phi
673              fPhiDist[i]->SetName(Form("phidist%d%d",i,fRunNumber));
674              fSparseDist->GetAxis(1)->SetRange(1,2); // reset axes
675              fSparseDist->GetAxis(2)->SetRange(1,2);
676              }
677              fSparseDist->GetAxis(0)->SetRange(1,2901);// reset run axis
678        }
679
680     if (!fPhiDist[0]) AliFatal(Form("Cannot find OADB phi distribution for run %d", fRunNumber));
681
682   } 
683   
684     
685   if (fPeriod.CompareTo("LHC10h")==0 || fUserphidist){
686      Bool_t emptybins;
687
688      int iter = 0;  
689      while (iter<3){
690          emptybins = kFALSE;
691    
692          for (int i=1; i<fPhiDist[0]->GetNbinsX(); i++){
693            if (!((fPhiDist[0]->GetBinContent(i))>0)) {
694              emptybins = kTRUE;
695            }
696          }  
697          if (emptybins) {
698            cout << "empty bins - rebinning!" << endl;
699            fPhiDist[0]->Rebin();
700            iter++;
701          }      
702          else iter = 3;
703      }
704   
705      if (emptybins) {
706        AliError("After Maximum of rebinning still empty Phi-bins!!!");
707      }
708   }
709   if (fPeriod.CompareTo("LHC10h") != 0 && fPeriod.CompareTo("LHC11h") != 0 && !fUserphidist){
710   AliInfo("No Phi-weights available. All Phi weights set to 1");
711   SetUsePhiWeight(kFALSE);
712   }
713 }
714
715 //__________________________________________________________________________
716 void AliEPSelectionTask::SetPersonalPhiDistribution(const char* infilename, char* listname)
717 {
718   
719   fUserphidist = kTRUE;
720   
721   TFile f(infilename);
722   TObject* list = f.Get(listname);
723   fPhiDist[0] = (TH1F*)list->FindObject("fHOutPhi");
724   if (!fPhiDist[0]) AliFatal("Phi Distribution not found!!!");
725
726   f.Close();
727
728
729
730 //_________________________________________________________________________
731 TObjArray* AliEPSelectionTask::GetAODTracksAndMaxID(AliAODEvent* aod, Int_t& maxid)
732 {
733   TObjArray *acctracks = new TObjArray();
734   
735   AliAODTrack *tr = 0;
736   Int_t maxid1 = 0;
737   Int_t maxidtemp = -1;
738   Float_t ptlow = 0;
739   Float_t ptup = 0;
740   Float_t etalow = 0;
741   Float_t etaup = 0;
742   fESDtrackCuts->GetPtRange(ptlow,ptup);
743   fESDtrackCuts->GetEtaRange(etalow,etaup);
744   
745   for (Int_t i = 0; i < aod->GetNumberOfTracks() ; i++){
746      tr = aod->GetTrack(i);
747      maxidtemp = tr->GetID(); 
748      if(maxidtemp < 0 && fAODfilterbit != 128) continue;
749      if(maxidtemp > -1 && fAODfilterbit == 128) continue;
750      if (fAODfilterbit == 128) maxidtemp = maxidtemp*(-1) - 1;
751      if (maxidtemp > maxid1) maxid1 = maxidtemp;
752      if(tr->TestFilterBit(fAODfilterbit) && tr->Pt() < ptup && tr->Pt() > ptlow && tr->Eta() < etaup && tr->Eta() > etalow){
753      acctracks->Add(tr);
754      }
755   }
756   
757   maxid = maxid1;
758   
759   return acctracks;
760   
761 }
762
763 //_________________________________________________________________________
764 void AliEPSelectionTask::SetOADBandPeriod()
765
766   TString oadbfilename; 
767
768   if (fRunNumber >= 136851 && fRunNumber <= 139515) // LHC10h
769     {fPeriod = "LHC10h";
770      if (!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
771         
772         if (fAnalysisInput.CompareTo("AOD")==0){
773            oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.aod.root", AliAnalysisManager::GetOADBPath()));
774            } else if (fAnalysisInput.CompareTo("ESD")==0){
775            oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath()));
776            }
777  
778        TFile foadb(oadbfilename); 
779        if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
780
781        AliInfo("Using Standard OADB");
782        fEPContainer = (AliOADBContainer*) foadb.Get("epphidist");    
783        if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection");
784        foadb.Close();
785        }
786      }
787   
788   if (fRunNumber >= 166529 && fRunNumber <= 170593) // LHC11h
789      {fPeriod = "LHC11h";
790       if (!fUserphidist) {
791       // if it's already set and custom class is required, we use the one provided by the user
792       
793       oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist2011.root", AliAnalysisManager::GetOADBPath()));
794       TFile *foadb = TFile::Open(oadbfilename); 
795       if(!foadb->IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
796
797       AliInfo("Using Standard OADB");
798       fSparseDist = (THnSparse*) foadb->Get("Default");    
799       if (!fSparseDist) AliFatal("Cannot fetch OADB container for EP selection");
800       foadb->Close();
801       if(!fHruns){ 
802            fHruns = (TH1F*)fSparseDist->Projection(0); //projection on run axis;
803            fHruns->SetName("runsHisto");
804       }
805       } 
806      } 
807 }
808
809 //_________________________________________________________________________
810 TH1F* AliEPSelectionTask::SelectPhiDist(AliVTrack *track)
811
812   if (fPeriod.CompareTo("LHC10h")==0  || fUserphidist) return fPhiDist[0];
813   else if(fPeriod.CompareTo("LHC11h")==0)
814     {
815      if (track->Charge() < 0)
816        {
817         if(track->Eta() < 0.)       return fPhiDist[0];
818         else if (track->Eta() > 0.) return fPhiDist[2];
819        }
820       else if (track->Charge() > 0)
821        {
822         if(track->Eta() < 0.)       return fPhiDist[1];
823         else if (track->Eta() > 0.) return fPhiDist[3];
824        }
825        
826     }
827   return 0;   
828 }
829
830 TObjArray* AliEPSelectionTask::GetTracksForLHC11h(AliESDEvent* esd)
831 {
832   // Need to do this hack beacuse only this type of TPC only tracks in AOD is available and one type of Phi-weights is used
833   TObjArray *acctracks = new TObjArray();
834   acctracks->SetOwner(kTRUE);
835   
836   const AliESDVertex *vtxSPD = esd->GetPrimaryVertexSPD();
837
838   Float_t ptlow = 0;
839   Float_t ptup = 0;
840   Float_t etalow = 0;
841   Float_t etaup = 0;
842   fESDtrackCuts->GetPtRange(ptlow,ptup);
843   fESDtrackCuts->GetEtaRange(etalow,etaup);
844   
845   Double_t pDCA[3] = { 0. }; // momentum at DCA
846   Double_t rDCA[3] = { 0. }; // position at DCA
847   Float_t  dDCA[2] = {0.};    // DCA to the vertex d and z
848   Float_t  cDCA[3] = {0.};    // covariance of impact parameters
849
850  
851   
852   for (Int_t i = 0; i < esd->GetNumberOfTracks() ; i++){
853   
854     AliESDtrack* esdTrack = esd->GetTrack(i); //carefull do not modify it othwise  need to work with a copy 
855     //
856     // Track selection
857     if (!fESDtrackCuts->AcceptTrack(esdTrack)) continue;
858     
859     // create a tpc only tracl
860     AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(esd),esdTrack->GetID());
861     if(!track) continue;
862     
863     if(track->Pt()>0.)
864     {
865       // only constrain tracks above threshold
866       AliExternalTrackParam exParam;
867       // take the B-field from the ESD, no 3D fieldMap available at this point
868       Bool_t relate = false;
869       relate = track->RelateToVertexTPC(vtxSPD,esd->GetMagneticField(),kVeryBig,&exParam);
870       if(!relate){
871         delete track;
872         continue;
873       }
874       // fetch the track parameters at the DCA (unconstraint)
875       if(track->GetTPCInnerParam()){
876         track->GetTPCInnerParam()->GetPxPyPz(pDCA);
877         track->GetTPCInnerParam()->GetXYZ(rDCA);
878       }
879       // get the DCA to the vertex:
880       track->GetImpactParametersTPC(dDCA,cDCA);
881       // set the constrained parameters to the track
882       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
883     }
884     
885     
886     Float_t eta = track->Eta();
887     Float_t pT = track->Pt();
888
889     if(pT < ptlow || pT > ptup || eta < etalow || eta > etaup){
890       delete track;
891       continue;
892     }
893
894     acctracks->Add(track);
895    }
896   
897   
898   return acctracks;
899   
900 }
901