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