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