I did some modifications on the Event plane task, the diff to today's trunk is attach...
[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 <TProfile.h>
29 #include <TFile.h>
30 #include <TObjString.h>
31 #include <TString.h>
32 #include <TCanvas.h>
33 #include <TROOT.h>
34 #include <TDirectory.h>
35 #include <TSystem.h>
36 #include <iostream>
37 #include <TRandom2.h>
38 #include <TArrayF.h>
39
40 #include "AliAnalysisManager.h"
41 #include "AliVEvent.h"
42 #include "AliESD.h"
43 #include "AliESDEvent.h"
44 #include "AliMCEvent.h"
45 #include "AliESDtrack.h"
46 #include "AliESDtrackCuts.h"
47 #include "AliESDHeader.h"
48 #include "AliESDInputHandler.h"
49 #include "AliAODHandler.h"
50 #include "AliAODEvent.h"
51 #include "AliHeader.h"
52 #include "AliGenHijingEventHeader.h"
53 #include "AliAnalysisTaskSE.h"
54 #include "AliPhysicsSelectionTask.h"
55 #include "AliPhysicsSelection.h"
56 #include "AliBackgroundSelection.h"
57 #include "AliESDUtils.h"
58 #include "AliOADBContainer.h"
59
60 #include "AliEventplane.h"
61
62 ClassImp(AliEPSelectionTask)
63
64 //________________________________________________________________________
65 AliEPSelectionTask::AliEPSelectionTask():
66 AliAnalysisTaskSE(),
67   fDebug(0),
68   fAnalysisInput("ESD"),
69   fTrackType("TPC"),
70   fUseMCRP(kFALSE),
71   fUsePhiWeight(kFALSE),
72   fUsePtWeight(kFALSE),
73   fSaveTrackContribution(kFALSE),
74   fuserphidist(kFALSE),
75   fusercuts(kFALSE),
76   frunNumber(-15),
77   fESDtrackCuts(0),
78   ftracklist(0),
79   fEPContainer(0),
80   fPhiDist(0),
81   fQVector(0),
82   fQContributionX(0),
83   fQContributionY(0),
84   fEventplaneQ(0),
85   fQsub1(0),
86   fQsub2(0),
87   fQsubRes(0),  
88   fOutputList(0),
89   fHOutEventplaneQ(0),
90   fHOutPhi(0),
91   fHOutPhiCorr(0),
92   fHOutsub1sub2(0),
93   fHOutNTEPRes(0),
94   fHOutPTPsi(0),
95   fHOutDiff(0),
96   fHOutleadPTPsi(0)
97 {   
98   // Default constructor
99   AliInfo("Event Plane Selection enabled.");
100 }   
101
102 //________________________________________________________________________
103 AliEPSelectionTask::AliEPSelectionTask(const char *name):
104   AliAnalysisTaskSE(name),
105   fDebug(0),
106   fAnalysisInput("ESD"),
107   fTrackType("TPC"),
108   fUseMCRP(kFALSE),
109   fUsePhiWeight(kFALSE),
110   fUsePtWeight(kFALSE),  
111   fSaveTrackContribution(kFALSE),
112   fuserphidist(kFALSE),
113   fusercuts(kFALSE),
114   frunNumber(-15),
115   fESDtrackCuts(0),
116   ftracklist(0),
117   fEPContainer(0),
118   fPhiDist(0),
119   fQVector(0),
120   fQContributionX(0),
121   fQContributionY(0),
122   fEventplaneQ(0),
123   fQsub1(0),
124   fQsub2(0),
125   fQsubRes(0),
126   fOutputList(0),
127   fHOutEventplaneQ(0),
128   fHOutPhi(0),
129   fHOutPhiCorr(0),
130   fHOutsub1sub2(0),
131   fHOutNTEPRes(0),
132   fHOutPTPsi(0),
133   fHOutDiff(0),
134   fHOutleadPTPsi(0)
135 {
136   // Default constructor
137   AliInfo("Event Plane Selection enabled.");
138   DefineOutput(1, TList::Class());
139 }
140  
141 //________________________________________________________________________
142 AliEPSelectionTask::~AliEPSelectionTask()
143 {
144   // Destructor  
145   if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()){
146       delete fOutputList;
147       fOutputList = 0;
148   }
149   if (fESDtrackCuts){
150       delete fESDtrackCuts;
151       fESDtrackCuts = 0;
152   }
153   if (fQVector){
154       delete fQVector;
155       fQVector = 0;
156   }
157   if (fQsub1){
158       delete fQsub1;
159       fQsub1 = 0;
160   }
161   if (fQsub2){
162       delete fQsub2;
163       fQsub2 = 0;
164   }
165   if (fPhiDist){
166       delete fPhiDist;
167       fPhiDist = 0;
168   }
169   if (ftracklist){
170       delete ftracklist;
171       ftracklist = 0;
172   }
173   if (fEPContainer){
174       delete fEPContainer;
175       fEPContainer = 0;
176   }
177 }  
178
179 //________________________________________________________________________
180 void AliEPSelectionTask::UserCreateOutputObjects()
181 {  
182   // Create the output containers
183   if (fDebug>1) printf("AliEPSelectionTask::UserCreateOutputObjects() \n");
184   AliLog::SetClassDebugLevel("AliEPSelectionTask", AliLog::kInfo);
185
186   fOutputList = new TList();
187   fOutputList->SetOwner();
188   fHOutEventplaneQ = new TH1F("fHOutEventplaneQ","fHOutEventplaneQ; Event Plane Q",100,0,TMath::Pi());
189   fHOutPhi = new TH1F("fHOutPhi","fHOutPhi; Phi Distribution",100,0,TMath::TwoPi());
190   fHOutPhiCorr = new TH1F("fHOutPhiCorr","fHOutPhiCorr; Corrected Phi Distribution",100,0,TMath::TwoPi());
191   fHOutsub1sub2 = new TH2F("fHOutsub1sub2","fHOutsub1sub2; EP1; EP2",100,0,TMath::Pi(),100,0,TMath::Pi());
192   fHOutNTEPRes = new TH2F("fHOutNTEPRes","fHOutNTEPRes; Number of Tracks; Event Plane Resolution",100,0,5000,100,-TMath::Pi(),TMath::Pi());
193   fHOutPTPsi = new TH2F("fHOutPTPsi","fHOutPTPsi; PT; Phi-EP",100,0,20,100,0,TMath::Pi());
194   fHOutDiff = new TH2F("fHOutDiff","fHOutDiff; EP; MCEP",100,0,TMath::Pi(),100,0,TMath::Pi());
195   fHOutleadPTPsi = new TH2F("fHOutleadPTPsi","fHOutleadPTPsi; leadPT; EP",100,0,TMath::Pi(),100,0,TMath::Pi());
196
197   fOutputList->Add(fHOutEventplaneQ);
198   fOutputList->Add(fHOutPhi);
199   fOutputList->Add(fHOutPhiCorr);
200   fOutputList->Add(fHOutsub1sub2);
201   fOutputList->Add(fHOutNTEPRes);
202   fOutputList->Add(fHOutPTPsi);
203   fOutputList->Add(fHOutleadPTPsi);
204   fOutputList->Add(fHOutDiff);
205   
206   PostData(1, fOutputList); 
207   
208     if(!fuserphidist) { // if it's already set and custom class is required, we use the one provided by the user
209
210     TString oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath()));
211
212     TFile foadb(oadbfilename); 
213     if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
214
215     AliInfo("Using Standard OADB");
216     fEPContainer = (AliOADBContainer*) foadb.Get("epphidist");    
217     if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection");
218     foadb.Close();
219     }
220 }
221
222 //________________________________________________________________________
223 void AliEPSelectionTask::UserExec(Option_t */*option*/)
224
225   // Execute analysis for current event:
226   if (fDebug>1) printf(" **** AliEPSelectionTask::UserExec() \n");
227   
228 //   frunNumber = -15;
229  
230   AliEventplane* esdEP = 0;
231   TVector2 QQ1;
232   TVector2 QQ2;
233   Double_t fRP = 0.; // the monte carlo reaction plane angle
234     
235   if (fAnalysisInput.CompareTo("ESD")==0){
236     
237     AliVEvent* event = InputEvent();
238     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
239     if (!(frunNumber == esd->GetRunNumber())) {
240       frunNumber = esd->GetRunNumber();
241       SetPhiDist();      
242     }
243     
244     
245     if (fUseMCRP) {
246       AliMCEvent* mcEvent  = MCEvent();      
247       if (mcEvent && mcEvent->GenEventHeader()) {
248         AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
249         if (headerH) fRP = headerH->ReactionPlaneAngle();
250       }
251     }
252     if (esd){
253       esdEP = esd->GetEventplane();
254       if (fSaveTrackContribution) {
255         esdEP->GetQContributionXArray()->Set(esd->GetNumberOfTracks());
256         esdEP->GetQContributionYArray()->Set(esd->GetNumberOfTracks());
257       }
258       
259       if (fTrackType.CompareTo("GLOBAL")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE);
260       if (fTrackType.CompareTo("TPC")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE);
261       const int NT = ftracklist->GetEntries();
262       
263       if (NT>4){
264         fQVector = new TVector2(GetQ(esdEP,ftracklist));
265         fEventplaneQ = fQVector->Phi()/2; 
266         GetQsub(QQ1, QQ2, ftracklist);
267         fQsub1 = new TVector2(QQ1);
268         fQsub2 = new TVector2(QQ2);
269         fQsubRes = (fQsub1->Phi()/2 - fQsub2->Phi()/2);
270         
271         esdEP->SetQVector(fQVector);
272         esdEP->SetEventplaneQ(fEventplaneQ);
273         esdEP->SetQsub(fQsub1,fQsub2);
274         esdEP->SetQsubRes(fQsubRes);
275         
276         fHOutEventplaneQ->Fill(fEventplaneQ);
277         fHOutsub1sub2->Fill(fQsub1->Phi()/2,fQsub2->Phi()/2);
278         fHOutNTEPRes->Fill(NT,fQsubRes);
279
280         if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
281         
282         for (int iter = 0; iter<NT;iter++){
283           AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter));
284           if (track) {
285             float delta = track->Phi()-fEventplaneQ;
286             while (delta < 0) delta += TMath::Pi();
287             while (delta > TMath::Pi()) delta -= TMath::Pi();
288             fHOutPTPsi->Fill(track->Pt(),delta);
289             fHOutPhi->Fill(track->Phi());
290             fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
291           }
292         }
293         
294         AliESDtrack* trmax = esd->GetTrack(0);
295         for (int iter = 1; iter<NT;iter++){
296           AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter));
297           if (track && (track->Pt() > trmax->Pt())) trmax = track;
298         }
299         fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);      
300       }
301       delete ftracklist;
302     }
303   }
304   
305   else if (fAnalysisInput.CompareTo("AOD")==0){
306     //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
307     // to be implemented
308     printf("  AOD analysis not yet implemented!!!\n\n");
309     return;
310   }  
311   else {
312     printf(" Analysis Input not known!\n\n ");
313     return;
314   }
315   PostData(1, fOutputList); 
316 }
317
318 //________________________________________________________________________
319 void AliEPSelectionTask::Terminate(Option_t */*option*/)
320 {
321   // Terminate analysis
322 }
323
324 //__________________________________________________________________________
325 TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
326 {
327   TVector2 mQ;
328   float mQx=0, mQy=0;
329   AliESDtrack* track;
330   Double_t weight;
331   
332   int NT = tracklist->GetEntries();
333
334   for (int i=0; i<NT; i++){
335     weight = 1;
336     track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
337     if (track) {
338       weight = GetWeight(track);
339       if (fSaveTrackContribution){
340         EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi()),track->GetID());
341         EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi()),track->GetID());
342       }
343       mQx += (weight*cos(2*track->Phi()));
344       mQy += (weight*sin(2*track->Phi()));
345     }
346   }
347   mQ.Set(mQx,mQy);
348   return mQ;
349 }
350   
351   //________________________________________________________________________
352 void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist)
353 {
354   TVector2 mQ[2];
355   float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
356   Double_t weight;
357
358   AliESDtrack* track;
359   TRandom2 RN = 0;
360   
361   int NT = tracklist->GetEntries();
362   int trackcounter1=0, trackcounter2=0;
363   
364   for (Int_t i = 0; i < NT; i++) {
365     weight = 1;
366     track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
367     if (!track) continue;
368     weight = GetWeight(track);
369     
370     // This loop splits the track set into 2 random subsets
371     if( trackcounter1 < int(NT/2.) && trackcounter2 < int(NT/2.)){
372       float random = RN.Rndm();
373       if(random < .5){
374         mQx1 += (weight*cos(2*track->Phi()));
375         mQy1 += (weight*sin(2*track->Phi()));
376         trackcounter1++;
377       }
378       else {
379         mQx2 += (weight*cos(2*track->Phi()));
380         mQy2 += (weight*sin(2*track->Phi()));
381         trackcounter2++;
382       }
383     }
384     else if( trackcounter1 >= int(NT/2.)){
385       mQx2 += (weight*cos(2*track->Phi()));
386       mQy2 += (weight*sin(2*track->Phi()));
387       trackcounter2++;
388     }
389     else {
390       mQx1 += (weight*cos(2*track->Phi()));
391       mQy1 += (weight*sin(2*track->Phi()));
392       trackcounter1++;
393     }
394   }
395   mQ[0].Set(mQx1,mQy1);
396   mQ[1].Set(mQx2,mQy2);
397   Q1 = mQ[0];
398   Q2 = mQ[1];
399 }
400
401 //________________________________________________________________________
402 void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
403   
404   fusercuts = kTRUE;
405   fESDtrackCuts = trackcuts;
406 }
407
408 //_____________________________________________________________________________
409 void AliEPSelectionTask::SetTrackType(TString tracktype){
410   fTrackType = tracktype;
411   if (!fusercuts) {
412   if (fTrackType.CompareTo("GLOBAL")==0) fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
413   if (fTrackType.CompareTo("TPC")==0)    fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
414   fESDtrackCuts->SetPtRange(0.15,20);
415   fESDtrackCuts->SetEtaRange(-0.8,0.8);
416   }
417 }
418
419 //________________________________________________________________________
420 Double_t AliEPSelectionTask::GetWeight(AliESDtrack* track)
421 {
422   Double_t ptweight=1;
423
424   if (fUsePtWeight) {      
425     if (track->Pt()<2) ptweight=track->Pt();
426     else ptweight=2;
427   }
428   return ptweight*GetPhiWeight(track);
429 }
430
431 //________________________________________________________________________
432 Double_t AliEPSelectionTask::GetPhiWeight(AliESDtrack* track)
433 {
434   Double_t phiweight=1;
435   
436   if (fUsePhiWeight && fPhiDist && track) {
437     Double_t nParticles = fPhiDist->Integral();
438     Double_t nPhibins = fPhiDist->GetNbinsX();
439   
440     Double_t Phi = track->Phi();
441     
442     while (Phi<0) Phi += TMath::TwoPi();
443     while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
444       
445     Double_t PhiDistValue = fPhiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
446     
447     if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
448   }
449   return phiweight;
450 }
451
452 //__________________________________________________________________________
453 void AliEPSelectionTask::SetPhiDist() 
454 {
455   if(!fuserphidist) { // if it's already set and custom class is required, we use the one provided by the user
456
457     fPhiDist = (TH1F*) fEPContainer->GetObject(frunNumber, "Default");
458     if (!fPhiDist) AliFatal(Form("Cannot find OADB phi distribution for run %d", frunNumber));
459
460   } 
461   else {
462     AliInfo("Using Custom Phi Distribution");
463   }
464     
465   Bool_t emptybins;
466
467   int iter = 0;  
468   while (iter<3){
469       emptybins = kFALSE;
470    
471       for (int i=1; i<fPhiDist->GetNbinsX(); i++){
472         if (!((fPhiDist->GetBinContent(i))>0)) {
473           emptybins = kTRUE;
474         }
475       }  
476       if (emptybins) {
477         cout << "empty bins - rebinning!" << endl;
478         fPhiDist->Rebin();
479         iter++;
480       }      
481       else iter = 3;
482   }
483   
484   if (emptybins) {
485     AliError("After Maximum of rebinning still empty Phi-bins!!!");
486   }
487 }
488
489 //__________________________________________________________________________
490 void AliEPSelectionTask::SetPersonalPhiDistribution(char* infilename, char* listname)
491 {
492   
493   fuserphidist = kTRUE;
494   
495   TFile f(infilename);
496   TObject* list = f.Get(listname);
497   fPhiDist = (TH1F*)list->FindObject("fHOutPhi");
498   if (!fPhiDist) AliFatal("Phi Distribution not found!!!");
499
500   f.Close();
501