99b7d61e68a5d0a85bf0928d97ef042441ea82f4
[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 if (fSplitMethod == AliEPSelectionTask::kCharge) {
535      
536     for (Int_t i = 0; i < nt; i++) {
537       weight = 1;
538       track = dynamic_cast<AliVTrack*> (tracklist->At(i));
539       if (!track) continue;
540       weight = GetWeight(track);
541       Short_t cha = track->Charge();
542       idtemp = track->GetID(); 
543       if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
544
545       if (cha > 0) {  
546         mQx1 += (weight*cos(2*track->Phi()));
547         mQy1 += (weight*sin(2*track->Phi()));
548         if (fSaveTrackContribution){
549           EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
550           EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
551         }
552       } else if (cha < 0) {
553         mQx2 += (weight*cos(2*track->Phi()));
554         mQy2 += (weight*sin(2*track->Phi()));
555         if (fSaveTrackContribution){
556           EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
557           EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
558         }
559       }
560     }
561   } else {
562     printf("plane resolution determination method not available!\n\n ");
563     return;
564   }
565      
566   mQ[0].Set(mQx1,mQy1);
567   mQ[1].Set(mQx2,mQy2);
568   Q1 = mQ[0];
569   Q2 = mQ[1];
570 }
571
572 //________________________________________________________________________
573 void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
574   
575   if(fESDtrackCuts){ 
576     delete fESDtrackCuts;
577     fESDtrackCuts = 0;
578   }
579   if (fAnalysisInput.CompareTo("AOD")==0){
580     AliInfo("ESD track cuts not possible for AOD analysis; please use SetPersonalAODtrackCuts(); using TPC only track cuts");  
581     fUsercuts = kFALSE;
582     SetTrackType("TPC");
583     return;
584   } 
585   fUsercuts = kTRUE;
586   fESDtrackCuts = trackcuts;
587 }
588
589 //________________________________________________________________________
590 void AliEPSelectionTask::SetPersonalAODtrackCuts(UInt_t filterbit, Float_t etalow, Float_t etaup, Float_t ptlow, Float_t ptup){
591   
592   if(fESDtrackCuts){ 
593     delete fESDtrackCuts;
594     fESDtrackCuts = 0;
595   }
596   if (fAnalysisInput.CompareTo("ESD")==0){
597     AliInfo("AOD track cuts not possible for ESD analysis; please use SetPersonalESDtrackCuts(); using TPC only track cuts");  
598     fUsercuts = kFALSE;
599     SetTrackType("TPC");
600     return;
601   }
602   fUsercuts = kTRUE;
603   fESDtrackCuts = new AliESDtrackCuts();
604   fESDtrackCuts->SetPtRange(ptlow,ptup);
605   fESDtrackCuts->SetEtaRange(etalow,etaup);
606   fAODfilterbit = filterbit;
607 }
608
609 //_____________________________________________________________________________
610
611 void AliEPSelectionTask::SetTrackType(TString tracktype){
612   fTrackType = tracktype;
613   if (!fUsercuts) {
614   if (fTrackType.CompareTo("GLOBAL")==0){ 
615     fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
616     fAODfilterbit = 32;
617   }     
618   if (fTrackType.CompareTo("TPC")==0){  
619     fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
620     fAODfilterbit = 128;
621   }
622   fESDtrackCuts->SetPtRange(0.15,20.);
623   fESDtrackCuts->SetEtaRange(-0.8,0.8);
624   }
625 }
626
627 //________________________________________________________________________
628 Double_t AliEPSelectionTask::GetWeight(TObject* track1)
629 {
630   Double_t ptweight=1;
631   AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
632   if (fUsePtWeight && track) {      
633     if (track->Pt()<2) ptweight=track->Pt();
634     else ptweight=2;
635   }
636   return ptweight*GetPhiWeight(track);
637 }
638
639 //________________________________________________________________________
640 Double_t AliEPSelectionTask::GetPhiWeight(TObject* track1)
641 {
642   Double_t phiweight=1;
643   AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
644
645   TH1F *phiDist = 0x0;
646   if(track) phiDist = SelectPhiDist(track);
647   
648   if (fUsePhiWeight && phiDist && track) {
649     Double_t nParticles = phiDist->Integral();
650     Double_t nPhibins = phiDist->GetNbinsX();
651   
652     Double_t Phi = track->Phi();
653     
654     while (Phi<0) Phi += TMath::TwoPi();
655     while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
656       
657     Double_t PhiDistValue = phiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
658     
659     if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
660   }
661   return phiweight;
662 }
663
664 //__________________________________________________________________________
665 void AliEPSelectionTask::SetPhiDist() 
666 {
667   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
668
669     if (fPeriod.CompareTo("LHC10h")==0)
670        {
671         fPhiDist[0] = (TH1F*) fEPContainer->GetObject(fRunNumber, "Default");}
672         else if(fPeriod.CompareTo("LHC11h")==0){
673             Int_t runbin=fHruns->FindBin(fRunNumber);
674             if (fHruns->GetBinContent(runbin) > 1){
675               fSparseDist->GetAxis(0)->SetRange(runbin,runbin);
676               }
677             else if(fHruns->GetBinContent(runbin) < 2){
678                fSparseDist->GetAxis(0)->SetRange(1,2901); // not calibrated run, use integrated phi-weights
679                AliInfo("Using integrated Phi-weights for this run");
680                }
681             for (Int_t i = 0; i<4 ;i++)
682             {
683              if(fPhiDist[i]){
684                delete fPhiDist[i];
685                fPhiDist[i] = 0x0;
686                }
687              if(i == 0){
688                fSparseDist->GetAxis(1)->SetRange(1,1);  // neg charge
689                fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta 
690              if(i == 1){
691                fSparseDist->GetAxis(1)->SetRange(2,2);  // pos charge
692                fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta
693              if(i == 2){
694                fSparseDist->GetAxis(1)->SetRange(1,1);  // neg charge
695                fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta 
696              if(i == 3){
697                fSparseDist->GetAxis(1)->SetRange(2,2);  // pos charge
698                fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta
699              fPhiDist[i] = (TH1F*)fSparseDist->Projection(3); // Projection on Phi
700              fPhiDist[i]->SetName(Form("phidist%d%d",i,fRunNumber));
701              fSparseDist->GetAxis(1)->SetRange(1,2); // reset axes
702              fSparseDist->GetAxis(2)->SetRange(1,2);
703              }
704              fSparseDist->GetAxis(0)->SetRange(1,2901);// reset run axis
705        }
706
707     if (!fPhiDist[0]) AliFatal(Form("Cannot find OADB phi distribution for run %d", fRunNumber));
708
709   } 
710   
711     
712   if (fPeriod.CompareTo("LHC10h")==0 || fUserphidist){
713      Bool_t emptybins;
714
715      int iter = 0;  
716      while (iter<3){
717          emptybins = kFALSE;
718    
719          for (int i=1; i<fPhiDist[0]->GetNbinsX(); i++){
720            if (!((fPhiDist[0]->GetBinContent(i))>0)) {
721              emptybins = kTRUE;
722            }
723          }  
724          if (emptybins) {
725            cout << "empty bins - rebinning!" << endl;
726            fPhiDist[0]->Rebin();
727            iter++;
728          }      
729          else iter = 3;
730      }
731   
732      if (emptybins) {
733        AliError("After Maximum of rebinning still empty Phi-bins!!!");
734      }
735   }
736   if (fPeriod.CompareTo("LHC10h") != 0 && fPeriod.CompareTo("LHC11h") != 0 && !fUserphidist){
737   AliInfo("No Phi-weights available. All Phi weights set to 1");
738   SetUsePhiWeight(kFALSE);
739   }
740 }
741
742 //__________________________________________________________________________
743 void AliEPSelectionTask::SetPersonalPhiDistribution(const char* infilename, char* listname)
744 {
745   
746   fUserphidist = kTRUE;
747   
748   TFile f(infilename);
749   TObject* list = f.Get(listname);
750   fPhiDist[0] = (TH1F*)list->FindObject("fHOutPhi");
751   if (!fPhiDist[0]) AliFatal("Phi Distribution not found!!!");
752
753   f.Close();
754
755
756
757 //_________________________________________________________________________
758 TObjArray* AliEPSelectionTask::GetAODTracksAndMaxID(AliAODEvent* aod, Int_t& maxid)
759 {
760   TObjArray *acctracks = new TObjArray();
761   
762   AliAODTrack *tr = 0;
763   Int_t maxid1 = 0;
764   Int_t maxidtemp = -1;
765   Float_t ptlow = 0;
766   Float_t ptup = 0;
767   Float_t etalow = 0;
768   Float_t etaup = 0;
769   fESDtrackCuts->GetPtRange(ptlow,ptup);
770   fESDtrackCuts->GetEtaRange(etalow,etaup);
771   
772   for (Int_t i = 0; i < aod->GetNumberOfTracks() ; i++){
773      tr = aod->GetTrack(i);
774      maxidtemp = tr->GetID(); 
775      if(maxidtemp < 0 && fAODfilterbit != 128) continue;
776      if(maxidtemp > -1 && fAODfilterbit == 128) continue;
777      if (fAODfilterbit == 128) maxidtemp = maxidtemp*(-1) - 1;
778      if (maxidtemp > maxid1) maxid1 = maxidtemp;
779      if(tr->TestFilterBit(fAODfilterbit) && tr->Pt() < ptup && tr->Pt() > ptlow && tr->Eta() < etaup && tr->Eta() > etalow){
780      acctracks->Add(tr);
781      }
782   }
783   
784   maxid = maxid1;
785   
786   return acctracks;
787   
788 }
789
790 //_________________________________________________________________________
791 void AliEPSelectionTask::SetOADBandPeriod()
792
793   TString oadbfilename; 
794
795   if (fRunNumber >= 136851 && fRunNumber <= 139515) // LHC10h
796     {fPeriod = "LHC10h";
797      if (!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
798         
799         if (fAnalysisInput.CompareTo("AOD")==0){
800            oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.aod.root", AliAnalysisManager::GetOADBPath()));
801            } else if (fAnalysisInput.CompareTo("ESD")==0){
802            oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath()));
803            }
804  
805        TFile foadb(oadbfilename); 
806        if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
807
808        AliInfo("Using Standard OADB");
809        fEPContainer = (AliOADBContainer*) foadb.Get("epphidist");    
810        if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection");
811        foadb.Close();
812        }
813      }
814   
815   if (fRunNumber >= 166529 && fRunNumber <= 170593) // LHC11h
816      {fPeriod = "LHC11h";
817       if (!fUserphidist) {
818       // if it's already set and custom class is required, we use the one provided by the user
819       
820       oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist2011.root", AliAnalysisManager::GetOADBPath()));
821       TFile *foadb = TFile::Open(oadbfilename); 
822       if(!foadb->IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
823
824       AliInfo("Using Standard OADB");
825       fSparseDist = (THnSparse*) foadb->Get("Default");    
826       if (!fSparseDist) AliFatal("Cannot fetch OADB container for EP selection");
827       foadb->Close();
828       if(!fHruns){ 
829            fHruns = (TH1F*)fSparseDist->Projection(0); //projection on run axis;
830            fHruns->SetName("runsHisto");
831       }
832       } 
833      } 
834 }
835
836 //_________________________________________________________________________
837 TH1F* AliEPSelectionTask::SelectPhiDist(AliVTrack *track)
838
839   if (fPeriod.CompareTo("LHC10h")==0  || fUserphidist) return fPhiDist[0];
840   else if(fPeriod.CompareTo("LHC11h")==0)
841     {
842      if (track->Charge() < 0)
843        {
844         if(track->Eta() < 0.)       return fPhiDist[0];
845         else if (track->Eta() > 0.) return fPhiDist[2];
846        }
847       else if (track->Charge() > 0)
848        {
849         if(track->Eta() < 0.)       return fPhiDist[1];
850         else if (track->Eta() > 0.) return fPhiDist[3];
851        }
852        
853     }
854   return 0;   
855 }
856
857 TObjArray* AliEPSelectionTask::GetTracksForLHC11h(AliESDEvent* esd)
858 {
859   // 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
860   TObjArray *acctracks = new TObjArray();
861   acctracks->SetOwner(kTRUE);
862   
863   const AliESDVertex *vtxSPD = esd->GetPrimaryVertexSPD();
864
865   Float_t ptlow = 0;
866   Float_t ptup = 0;
867   Float_t etalow = 0;
868   Float_t etaup = 0;
869   fESDtrackCuts->GetPtRange(ptlow,ptup);
870   fESDtrackCuts->GetEtaRange(etalow,etaup);
871   
872   Double_t pDCA[3] = { 0. }; // momentum at DCA
873   Double_t rDCA[3] = { 0. }; // position at DCA
874   Float_t  dDCA[2] = {0.};    // DCA to the vertex d and z
875   Float_t  cDCA[3] = {0.};    // covariance of impact parameters
876
877  
878   
879   for (Int_t i = 0; i < esd->GetNumberOfTracks() ; i++){
880   
881     AliESDtrack* esdTrack = esd->GetTrack(i); //carefull do not modify it othwise  need to work with a copy 
882     //
883     // Track selection
884     if (!fESDtrackCuts->AcceptTrack(esdTrack)) continue;
885     
886     // create a tpc only tracl
887     AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(esd),esdTrack->GetID());
888     if(!track) continue;
889     
890     if(track->Pt()>0.)
891     {
892       // only constrain tracks above threshold
893       AliExternalTrackParam exParam;
894       // take the B-field from the ESD, no 3D fieldMap available at this point
895       Bool_t relate = false;
896       relate = track->RelateToVertexTPC(vtxSPD,esd->GetMagneticField(),kVeryBig,&exParam);
897       if(!relate){
898         delete track;
899         continue;
900       }
901       // fetch the track parameters at the DCA (unconstraint)
902       if(track->GetTPCInnerParam()){
903         track->GetTPCInnerParam()->GetPxPyPz(pDCA);
904         track->GetTPCInnerParam()->GetXYZ(rDCA);
905       }
906       // get the DCA to the vertex:
907       track->GetImpactParametersTPC(dDCA,cDCA);
908       // set the constrained parameters to the track
909       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
910     }
911     
912     
913     Float_t eta = track->Eta();
914     Float_t pT = track->Pt();
915
916     if(pT < ptlow || pT > ptup || eta < etalow || eta > etaup){
917       delete track;
918       continue;
919     }
920
921     acctracks->Add(track);
922    }
923   
924   
925   return acctracks;
926   
927 }
928