T0 ESD convert:only time from triggered interaction used for T0 raw vertex calculations
[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       if(!esdEP) return; // protection against missing EP branch (nanoAODs)
336       esdEP->Reset(); 
337      
338       Int_t maxID = 0;
339       TObjArray* tracklist = GetAODTracksAndMaxID(aod,maxID);
340         
341       if (fSaveTrackContribution) {
342         esdEP->GetQContributionXArray()->Set(maxID+1);
343         esdEP->GetQContributionYArray()->Set(maxID+1);
344         esdEP->GetQContributionXArraysub1()->Set(maxID+1);
345         esdEP->GetQContributionYArraysub1()->Set(maxID+1);
346         esdEP->GetQContributionXArraysub2()->Set(maxID+1);
347         esdEP->GetQContributionYArraysub2()->Set(maxID+1);
348       }
349         
350       const int NT = tracklist->GetEntries();
351       
352       if (NT>4){
353         fQVector = new TVector2(GetQ(esdEP,tracklist));
354         fEventplaneQ = fQVector->Phi()/2.; 
355         GetQsub(qq1, qq2, tracklist, esdEP);
356         fQsub1 = new TVector2(qq1);
357         fQsub2 = new TVector2(qq2);
358         fQsubRes = (fQsub1->Phi()/2. - fQsub2->Phi()/2.);
359         
360         esdEP->SetQVector(fQVector);
361         esdEP->SetEventplaneQ(fEventplaneQ);
362         esdEP->SetQsub(fQsub1,fQsub2);
363         esdEP->SetQsubRes(fQsubRes);
364         
365         fHOutEventplaneQ->Fill(fEventplaneQ);
366         fHOutsub1sub2->Fill(fQsub1->Phi()/2.,fQsub2->Phi()/2.);
367         fHOutNTEPRes->Fill(NT,fQsubRes);
368         
369         if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
370         
371         for (int iter = 0; iter<NT;iter++){
372           AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
373           if (track) {
374             float delta = track->Phi()-fEventplaneQ;
375             while (delta < 0) delta += TMath::Pi();
376             while (delta > TMath::Pi()) delta -= TMath::Pi();
377             fHOutPTPsi->Fill(track->Pt(),delta);
378             fHOutPhi->Fill(track->Phi());
379             fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
380           }
381         }
382         
383         AliAODTrack* trmax = aod->GetTrack(0);
384         for (int iter = 1; iter<NT;iter++){
385           AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
386           if (track && (track->Pt() > trmax->Pt())) trmax = track;
387         }
388         fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);      
389       }     
390       delete tracklist;
391       tracklist = 0;
392     }   
393         
394     
395   }  
396
397   
398   else {
399     printf(" Analysis Input not known!\n\n ");
400     return;
401   }
402   PostData(1, fOutputList); 
403 }
404
405 //________________________________________________________________________
406 void AliEPSelectionTask::Terminate(Option_t */*option*/)
407 {
408   // Terminate analysis
409 }
410
411 //__________________________________________________________________________
412 TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
413 {
414 // Get the Q vector
415   TVector2 mQ;
416   float mQx=0, mQy=0;
417   AliVTrack* track;
418   Double_t weight;
419   Int_t idtemp = -1;
420   
421   int nt = tracklist->GetEntries();
422
423   for (int i=0; i<nt; i++){
424     weight = 1;
425     track = dynamic_cast<AliVTrack*> (tracklist->At(i));
426     if (track) {
427       weight = GetWeight(track);
428     if (fSaveTrackContribution){
429       idtemp = track->GetID(); 
430       if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
431       EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi()),idtemp);
432       EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi()),idtemp);
433      }
434      mQx += (weight*cos(2*track->Phi()));
435      mQy += (weight*sin(2*track->Phi()));
436     }
437   }
438   mQ.Set(mQx,mQy);
439   return mQ;
440 }
441   
442   //________________________________________________________________________
443 void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist,AliEventplane* EP)
444 {
445 // Get Qsub
446   TVector2 mQ[2];
447   float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
448   Double_t weight;
449
450   AliVTrack* track;
451   TRandom2 rn = 0;
452   
453   int nt = tracklist->GetEntries();
454   int trackcounter1=0, trackcounter2=0;
455   int idtemp = 0;
456
457   if (fSplitMethod == AliEPSelectionTask::kRandom){
458     
459     for (Int_t i = 0; i < nt; i++) {
460       weight = 1;
461       track = dynamic_cast<AliVTrack*> (tracklist->At(i));
462       if (!track) continue;
463       weight = GetWeight(track);
464       idtemp = track->GetID(); 
465       if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
466     
467       // This loop splits the track set into 2 random subsets
468       if( trackcounter1 < int(nt/2.) && trackcounter2 < int(nt/2.)){
469         float random = rn.Rndm();
470         if(random < .5){
471           mQx1 += (weight*cos(2*track->Phi()));
472           mQy1 += (weight*sin(2*track->Phi()));
473           if (fSaveTrackContribution){
474             EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
475             EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
476           }
477           trackcounter1++;
478         }
479         else {
480           mQx2 += (weight*cos(2*track->Phi()));
481           mQy2 += (weight*sin(2*track->Phi()));
482           if (fSaveTrackContribution){
483             EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
484             EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
485           }
486           trackcounter2++;
487         }
488       }
489       else if( trackcounter1 >= int(nt/2.)){
490         mQx2 += (weight*cos(2*track->Phi()));
491         mQy2 += (weight*sin(2*track->Phi()));
492         if (fSaveTrackContribution){
493           EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
494           EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
495         }
496         trackcounter2++;
497       }
498       else {
499         mQx1 += (weight*cos(2*track->Phi()));
500         mQy1 += (weight*sin(2*track->Phi()));
501         if (fSaveTrackContribution){
502           EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
503           EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
504         }
505         trackcounter1++;
506       }
507     }
508   } else if (fSplitMethod == AliEPSelectionTask::kEta) {
509      
510     for (Int_t i = 0; i < nt; i++) {
511       weight = 1;
512       track = dynamic_cast<AliVTrack*> (tracklist->At(i));
513       if (!track) continue;
514       weight = GetWeight(track);
515       Double_t eta = track->Eta();
516       idtemp = track->GetID(); 
517       if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
518
519       if (eta > fEtaGap/2.) {  
520         mQx1 += (weight*cos(2*track->Phi()));
521         mQy1 += (weight*sin(2*track->Phi()));
522         if (fSaveTrackContribution){
523           EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
524           EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
525         }
526       } else if (eta < -1.*fEtaGap/2.) {
527         mQx2 += (weight*cos(2*track->Phi()));
528         mQy2 += (weight*sin(2*track->Phi()));
529         if (fSaveTrackContribution){
530           EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
531           EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
532         }
533       }
534     }
535   } else if (fSplitMethod == AliEPSelectionTask::kCharge) {
536      
537     for (Int_t i = 0; i < nt; i++) {
538       weight = 1;
539       track = dynamic_cast<AliVTrack*> (tracklist->At(i));
540       if (!track) continue;
541       weight = GetWeight(track);
542       Short_t cha = track->Charge();
543       idtemp = track->GetID(); 
544       if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
545
546       if (cha > 0) {  
547         mQx1 += (weight*cos(2*track->Phi()));
548         mQy1 += (weight*sin(2*track->Phi()));
549         if (fSaveTrackContribution){
550           EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
551           EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
552         }
553       } else if (cha < 0) {
554         mQx2 += (weight*cos(2*track->Phi()));
555         mQy2 += (weight*sin(2*track->Phi()));
556         if (fSaveTrackContribution){
557           EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
558           EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
559         }
560       }
561     }
562   } else {
563     printf("plane resolution determination method not available!\n\n ");
564     return;
565   }
566      
567   mQ[0].Set(mQx1,mQy1);
568   mQ[1].Set(mQx2,mQy2);
569   Q1 = mQ[0];
570   Q2 = mQ[1];
571 }
572
573 //________________________________________________________________________
574 void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
575   
576   if(fESDtrackCuts){ 
577     delete fESDtrackCuts;
578     fESDtrackCuts = 0;
579   }
580   if (fAnalysisInput.CompareTo("AOD")==0){
581     AliInfo("ESD track cuts not possible for AOD analysis; please use SetPersonalAODtrackCuts(); using TPC only track cuts");  
582     fUsercuts = kFALSE;
583     SetTrackType("TPC");
584     return;
585   } 
586   fUsercuts = kTRUE;
587   fESDtrackCuts = trackcuts;
588 }
589
590 //________________________________________________________________________
591 void AliEPSelectionTask::SetPersonalAODtrackCuts(UInt_t filterbit, Float_t etalow, Float_t etaup, Float_t ptlow, Float_t ptup, Int_t ntpc){
592   
593   if(fESDtrackCuts){ 
594     delete fESDtrackCuts;
595     fESDtrackCuts = 0;
596   }
597   if (fAnalysisInput.CompareTo("ESD")==0){
598     AliInfo("AOD track cuts not possible for ESD analysis; please use SetPersonalESDtrackCuts(); using TPC only track cuts");  
599     fUsercuts = kFALSE;
600     SetTrackType("TPC");
601     return;
602   }
603   fUsercuts = kTRUE;
604   fESDtrackCuts = new AliESDtrackCuts();
605   fESDtrackCuts->SetPtRange(ptlow,ptup);
606   fESDtrackCuts->SetMinNClustersTPC(ntpc);
607   fESDtrackCuts->SetEtaRange(etalow,etaup);
608   fAODfilterbit = filterbit;
609 }
610
611 //_____________________________________________________________________________
612
613 void AliEPSelectionTask::SetTrackType(TString tracktype){
614   fTrackType = tracktype;
615   if (!fUsercuts) {
616   if (fTrackType.CompareTo("GLOBAL")==0){ 
617     fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
618     fAODfilterbit = 32;
619   }     
620   if (fTrackType.CompareTo("TPC")==0){  
621     fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
622     fAODfilterbit = 128;
623   }
624   fESDtrackCuts->SetPtRange(0.15,20.);
625   fESDtrackCuts->SetEtaRange(-0.8,0.8);
626   }
627 }
628
629 //________________________________________________________________________
630 Double_t AliEPSelectionTask::GetWeight(TObject* track1)
631 {
632   Double_t ptweight=1;
633   AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
634   if (fUsePtWeight && track) {      
635     if (track->Pt()<2) ptweight=track->Pt();
636     else ptweight=2;
637   }
638   return ptweight*GetPhiWeight(track);
639 }
640
641 //________________________________________________________________________
642 Double_t AliEPSelectionTask::GetPhiWeight(TObject* track1)
643 {
644   Double_t phiweight=1;
645   AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
646
647   TH1F *phiDist = 0x0;
648   if(track) phiDist = SelectPhiDist(track);
649   
650   if (fUsePhiWeight && phiDist && track) {
651     Double_t nParticles = phiDist->Integral();
652     Double_t nPhibins = phiDist->GetNbinsX();
653   
654     Double_t Phi = track->Phi();
655     
656     while (Phi<0) Phi += TMath::TwoPi();
657     while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
658       
659     Double_t PhiDistValue = phiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
660     
661     if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
662   }
663   return phiweight;
664 }
665
666 //__________________________________________________________________________
667 void AliEPSelectionTask::SetPhiDist() 
668 {
669   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
670
671     if (fPeriod.CompareTo("LHC10h")==0)
672        {
673         fPhiDist[0] = (TH1F*) fEPContainer->GetObject(fRunNumber, "Default");}
674         else if(fPeriod.CompareTo("LHC11h")==0){
675             Int_t runbin=fHruns->FindBin(fRunNumber);
676             if (fHruns->GetBinContent(runbin) > 1){
677               fSparseDist->GetAxis(0)->SetRange(runbin,runbin);
678               }
679             else if(fHruns->GetBinContent(runbin) < 2){
680                fSparseDist->GetAxis(0)->SetRange(1,2901); // not calibrated run, use integrated phi-weights
681                AliInfo("Using integrated Phi-weights for this run");
682                }
683             for (Int_t i = 0; i<4 ;i++)
684             {
685              if(fPhiDist[i]){
686                delete fPhiDist[i];
687                fPhiDist[i] = 0x0;
688                }
689              if(i == 0){
690                fSparseDist->GetAxis(1)->SetRange(1,1);  // neg charge
691                fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta 
692              if(i == 1){
693                fSparseDist->GetAxis(1)->SetRange(2,2);  // pos charge
694                fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta
695              if(i == 2){
696                fSparseDist->GetAxis(1)->SetRange(1,1);  // neg charge
697                fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta 
698              if(i == 3){
699                fSparseDist->GetAxis(1)->SetRange(2,2);  // pos charge
700                fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta
701              fPhiDist[i] = (TH1F*)fSparseDist->Projection(3); // Projection on Phi
702              fPhiDist[i]->SetName(Form("phidist%d%d",i,fRunNumber));
703              fSparseDist->GetAxis(1)->SetRange(1,2); // reset axes
704              fSparseDist->GetAxis(2)->SetRange(1,2);
705              }
706              fSparseDist->GetAxis(0)->SetRange(1,2901);// reset run axis
707        }
708
709     if (!fPhiDist[0]) AliFatal(Form("Cannot find OADB phi distribution for run %d", fRunNumber));
710
711   } 
712   
713     
714   if (fPeriod.CompareTo("LHC10h")==0 || fUserphidist){
715      Bool_t emptybins;
716
717      int iter = 0;  
718      while (iter<3){
719          emptybins = kFALSE;
720    
721          for (int i=1; i<fPhiDist[0]->GetNbinsX(); i++){
722            if (!((fPhiDist[0]->GetBinContent(i))>0)) {
723              emptybins = kTRUE;
724            }
725          }  
726          if (emptybins) {
727            cout << "empty bins - rebinning!" << endl;
728            fPhiDist[0]->Rebin();
729            iter++;
730          }      
731          else iter = 3;
732      }
733   
734      if (emptybins) {
735        AliError("After Maximum of rebinning still empty Phi-bins!!!");
736      }
737   }
738   if (fPeriod.CompareTo("LHC10h") != 0 && fPeriod.CompareTo("LHC11h") != 0 && !fUserphidist){
739   AliInfo("No Phi-weights available. All Phi weights set to 1");
740   SetUsePhiWeight(kFALSE);
741   }
742 }
743
744 //__________________________________________________________________________
745 void AliEPSelectionTask::SetPersonalPhiDistribution(const char* infilename, char* listname)
746 {
747   
748   fUserphidist = kTRUE;
749   
750   TFile f(infilename);
751   TObject* list = f.Get(listname);
752   fPhiDist[0] = (TH1F*)list->FindObject("fHOutPhi");
753   if (!fPhiDist[0]) AliFatal("Phi Distribution not found!!!");
754
755   f.Close();
756
757
758
759 //_________________________________________________________________________
760 TObjArray* AliEPSelectionTask::GetAODTracksAndMaxID(AliAODEvent* aod, Int_t& maxid)
761 {
762   TObjArray *acctracks = new TObjArray();
763   
764   AliAODTrack *tr = 0;
765   Int_t maxid1 = 0;
766   Int_t maxidtemp = -1;
767   Float_t ptlow = 0;
768   Float_t ptup = 0;
769   Float_t etalow = 0;
770   Float_t etaup = 0;
771   fESDtrackCuts->GetPtRange(ptlow,ptup);
772   fESDtrackCuts->GetEtaRange(etalow,etaup);
773   Int_t ntpc = fESDtrackCuts->GetMinNClusterTPC(); 
774   
775   for (Int_t i = 0; i < aod->GetNumberOfTracks() ; i++){
776      tr = aod->GetTrack(i);
777      maxidtemp = tr->GetID(); 
778      if(maxidtemp < 0 && fAODfilterbit != 128) continue;
779      if(maxidtemp > -1 && fAODfilterbit == 128) continue;
780      if (fAODfilterbit == 128) maxidtemp = maxidtemp*(-1) - 1;
781      if (maxidtemp > maxid1) maxid1 = maxidtemp;
782      if(tr->TestFilterBit(fAODfilterbit) && tr->Pt() < ptup && tr->Pt() > ptlow && tr->Eta() < etaup && tr->Eta() > etalow && tr->GetTPCNcls() > ntpc){
783      acctracks->Add(tr);
784      }
785   }
786   
787   maxid = maxid1;
788   
789   return acctracks;
790   
791 }
792
793 //_________________________________________________________________________
794 void AliEPSelectionTask::SetOADBandPeriod()
795
796   TString oadbfilename; 
797
798   if (fRunNumber >= 136851 && fRunNumber <= 139515) // LHC10h
799     {fPeriod = "LHC10h";
800      if (!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
801         
802         if (fAnalysisInput.CompareTo("AOD")==0){
803            oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.aod.root", AliAnalysisManager::GetOADBPath()));
804            } else if (fAnalysisInput.CompareTo("ESD")==0){
805            oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath()));
806            }
807  
808        TFile foadb(oadbfilename); 
809        if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
810
811        AliInfo("Using Standard OADB");
812        fEPContainer = (AliOADBContainer*) foadb.Get("epphidist");    
813        if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection");
814        foadb.Close();
815        }
816      }
817   
818   if (fRunNumber >= 166529 && fRunNumber <= 170593) // LHC11h
819      {fPeriod = "LHC11h";
820       if (!fUserphidist) {
821       // if it's already set and custom class is required, we use the one provided by the user
822       
823       oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist2011.root", AliAnalysisManager::GetOADBPath()));
824       TFile *foadb = TFile::Open(oadbfilename); 
825       if(!foadb->IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
826
827       AliInfo("Using Standard OADB");
828       fSparseDist = (THnSparse*) foadb->Get("Default");    
829       if (!fSparseDist) AliFatal("Cannot fetch OADB container for EP selection");
830       foadb->Close();
831       if(!fHruns){ 
832            fHruns = (TH1F*)fSparseDist->Projection(0); //projection on run axis;
833            fHruns->SetName("runsHisto");
834       }
835       } 
836      } 
837 }
838
839 //_________________________________________________________________________
840 TH1F* AliEPSelectionTask::SelectPhiDist(AliVTrack *track)
841
842   if (fPeriod.CompareTo("LHC10h")==0  || fUserphidist) return fPhiDist[0];
843   else if(fPeriod.CompareTo("LHC11h")==0)
844     {
845      if (track->Charge() < 0)
846        {
847         if(track->Eta() < 0.)       return fPhiDist[0];
848         else if (track->Eta() > 0.) return fPhiDist[2];
849        }
850       else if (track->Charge() > 0)
851        {
852         if(track->Eta() < 0.)       return fPhiDist[1];
853         else if (track->Eta() > 0.) return fPhiDist[3];
854        }
855        
856     }
857   return 0;   
858 }
859
860 TObjArray* AliEPSelectionTask::GetTracksForLHC11h(AliESDEvent* esd)
861 {
862   // 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
863   TObjArray *acctracks = new TObjArray();
864   acctracks->SetOwner(kTRUE);
865   
866   const AliESDVertex *vtxSPD = esd->GetPrimaryVertexSPD();
867
868   Float_t ptlow = 0;
869   Float_t ptup = 0;
870   Float_t etalow = 0;
871   Float_t etaup = 0;
872   fESDtrackCuts->GetPtRange(ptlow,ptup);
873   fESDtrackCuts->GetEtaRange(etalow,etaup);
874   
875   Double_t pDCA[3] = { 0. }; // momentum at DCA
876   Double_t rDCA[3] = { 0. }; // position at DCA
877   Float_t  dDCA[2] = {0.};    // DCA to the vertex d and z
878   Float_t  cDCA[3] = {0.};    // covariance of impact parameters
879
880  
881   
882   for (Int_t i = 0; i < esd->GetNumberOfTracks() ; i++){
883   
884     AliESDtrack* esdTrack = esd->GetTrack(i); //carefull do not modify it othwise  need to work with a copy 
885     //
886     // Track selection
887     if (!fESDtrackCuts->AcceptTrack(esdTrack)) continue;
888     
889     // create a tpc only tracl
890     AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(esd),esdTrack->GetID());
891     if(!track) continue;
892     
893     if(track->Pt()>0.)
894     {
895       // only constrain tracks above threshold
896       AliExternalTrackParam exParam;
897       // take the B-field from the ESD, no 3D fieldMap available at this point
898       Bool_t relate = false;
899       relate = track->RelateToVertexTPC(vtxSPD,esd->GetMagneticField(),kVeryBig,&exParam);
900       if(!relate){
901         delete track;
902         continue;
903       }
904       // fetch the track parameters at the DCA (unconstraint)
905       if(track->GetTPCInnerParam()){
906         track->GetTPCInnerParam()->GetPxPyPz(pDCA);
907         track->GetTPCInnerParam()->GetXYZ(rDCA);
908       }
909       // get the DCA to the vertex:
910       track->GetImpactParametersTPC(dDCA,cDCA);
911       // set the constrained parameters to the track
912       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
913     }
914     
915     
916     Float_t eta = track->Eta();
917     Float_t pT = track->Pt();
918
919     if(pT < ptlow || pT > ptup || eta < etalow || eta > etaup){
920       delete track;
921       continue;
922     }
923
924     acctracks->Add(track);
925    }
926   
927   
928   return acctracks;
929   
930 }
931