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