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