1 /**************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //*****************************************************
17 // Class AliEPSelectionTask
18 // Class to determine event plane
19 // author: Alberica Toia, Johanna Gramling
20 //*****************************************************
22 #include "AliEPSelectionTask.h"
30 #include <TObjString.h>
34 #include <TDirectory.h>
40 #include "AliAnalysisManager.h"
41 #include "AliVEvent.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"
60 #include "AliEventplane.h"
62 ClassImp(AliEPSelectionTask)
64 //________________________________________________________________________
65 AliEPSelectionTask::AliEPSelectionTask():
68 fAnalysisInput("ESD"),
71 fUsePhiWeight(kFALSE),
73 fSaveTrackContribution(kFALSE),
98 // Default constructor
99 AliInfo("Event Plane Selection enabled.");
102 //________________________________________________________________________
103 AliEPSelectionTask::AliEPSelectionTask(const char *name):
104 AliAnalysisTaskSE(name),
106 fAnalysisInput("ESD"),
109 fUsePhiWeight(kFALSE),
110 fUsePtWeight(kFALSE),
111 fSaveTrackContribution(kFALSE),
112 fuserphidist(kFALSE),
136 // Default constructor
137 AliInfo("Event Plane Selection enabled.");
138 DefineOutput(1, TList::Class());
141 //________________________________________________________________________
142 AliEPSelectionTask::~AliEPSelectionTask()
145 if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()){
150 delete fESDtrackCuts;
179 //________________________________________________________________________
180 void AliEPSelectionTask::UserCreateOutputObjects()
182 // Create the output containers
183 if (fDebug>1) printf("AliEPSelectionTask::UserCreateOutputObjects() \n");
184 AliLog::SetClassDebugLevel("AliEPSelectionTask", AliLog::kInfo);
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());
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);
206 PostData(1, fOutputList);
208 if(!fuserphidist) { // if it's already set and custom class is required, we use the one provided by the user
210 TString oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath()));
212 TFile foadb(oadbfilename);
213 if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
215 AliInfo("Using Standard OADB");
216 fEPContainer = (AliOADBContainer*) foadb.Get("epphidist");
217 if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection");
222 //________________________________________________________________________
223 void AliEPSelectionTask::UserExec(Option_t */*option*/)
225 // Execute analysis for current event:
226 if (fDebug>1) printf(" **** AliEPSelectionTask::UserExec() \n");
230 AliEventplane* esdEP = 0;
233 Double_t fRP = 0.; // the monte carlo reaction plane angle
235 if (fAnalysisInput.CompareTo("ESD")==0){
237 AliVEvent* event = InputEvent();
238 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
240 if (!(frunNumber == esd->GetRunNumber())) {
241 frunNumber = esd->GetRunNumber();
247 AliMCEvent* mcEvent = MCEvent();
248 if (mcEvent && mcEvent->GenEventHeader()) {
249 AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
250 if (headerH) fRP = headerH->ReactionPlaneAngle();
254 esdEP = esd->GetEventplane();
255 if (fSaveTrackContribution) {
256 esdEP->GetQContributionXArray()->Set(esd->GetNumberOfTracks());
257 esdEP->GetQContributionYArray()->Set(esd->GetNumberOfTracks());
260 if (fTrackType.CompareTo("GLOBAL")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE);
261 if (fTrackType.CompareTo("TPC")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE);
262 const int NT = ftracklist->GetEntries();
265 fQVector = new TVector2(GetQ(esdEP,ftracklist));
266 fEventplaneQ = fQVector->Phi()/2;
267 GetQsub(QQ1, QQ2, ftracklist);
268 fQsub1 = new TVector2(QQ1);
269 fQsub2 = new TVector2(QQ2);
270 fQsubRes = (fQsub1->Phi()/2 - fQsub2->Phi()/2);
272 esdEP->SetQVector(fQVector);
273 esdEP->SetEventplaneQ(fEventplaneQ);
274 esdEP->SetQsub(fQsub1,fQsub2);
275 esdEP->SetQsubRes(fQsubRes);
277 fHOutEventplaneQ->Fill(fEventplaneQ);
278 fHOutsub1sub2->Fill(fQsub1->Phi()/2,fQsub2->Phi()/2);
279 fHOutNTEPRes->Fill(NT,fQsubRes);
281 if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
283 for (int iter = 0; iter<NT;iter++){
284 AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter));
286 float delta = track->Phi()-fEventplaneQ;
287 while (delta < 0) delta += TMath::Pi();
288 while (delta > TMath::Pi()) delta -= TMath::Pi();
289 fHOutPTPsi->Fill(track->Pt(),delta);
290 fHOutPhi->Fill(track->Phi());
291 fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
295 AliESDtrack* trmax = esd->GetTrack(0);
296 for (int iter = 1; iter<NT;iter++){
297 AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter));
298 if (track && (track->Pt() > trmax->Pt())) trmax = track;
300 fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);
306 else if (fAnalysisInput.CompareTo("AOD")==0){
307 //AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
309 printf(" AOD analysis not yet implemented!!!\n\n");
313 printf(" Analysis Input not known!\n\n ");
316 PostData(1, fOutputList);
319 //________________________________________________________________________
320 void AliEPSelectionTask::Terminate(Option_t */*option*/)
322 // Terminate analysis
325 //__________________________________________________________________________
326 TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
333 int NT = tracklist->GetEntries();
335 for (int i=0; i<NT; i++){
337 track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
339 weight = GetWeight(track);
340 if (fSaveTrackContribution){
341 EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi()),track->GetID());
342 EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi()),track->GetID());
344 mQx += (weight*cos(2*track->Phi()));
345 mQy += (weight*sin(2*track->Phi()));
352 //________________________________________________________________________
353 void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist)
356 float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
362 int NT = tracklist->GetEntries();
363 int trackcounter1=0, trackcounter2=0;
365 for (Int_t i = 0; i < NT; i++) {
367 track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
368 if (!track) continue;
369 weight = GetWeight(track);
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();
375 mQx1 += (weight*cos(2*track->Phi()));
376 mQy1 += (weight*sin(2*track->Phi()));
380 mQx2 += (weight*cos(2*track->Phi()));
381 mQy2 += (weight*sin(2*track->Phi()));
385 else if( trackcounter1 >= int(NT/2.)){
386 mQx2 += (weight*cos(2*track->Phi()));
387 mQy2 += (weight*sin(2*track->Phi()));
391 mQx1 += (weight*cos(2*track->Phi()));
392 mQy1 += (weight*sin(2*track->Phi()));
396 mQ[0].Set(mQx1,mQy1);
397 mQ[1].Set(mQx2,mQy2);
402 //________________________________________________________________________
403 void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
406 fESDtrackCuts = trackcuts;
409 //_____________________________________________________________________________
410 void AliEPSelectionTask::SetTrackType(TString tracktype){
411 fTrackType = tracktype;
413 if (fTrackType.CompareTo("GLOBAL")==0) fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
414 if (fTrackType.CompareTo("TPC")==0) fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
415 fESDtrackCuts->SetPtRange(0.15,20);
416 fESDtrackCuts->SetEtaRange(-0.8,0.8);
420 //________________________________________________________________________
421 Double_t AliEPSelectionTask::GetWeight(AliESDtrack* track)
426 if (track->Pt()<2) ptweight=track->Pt();
429 return ptweight*GetPhiWeight(track);
432 //________________________________________________________________________
433 Double_t AliEPSelectionTask::GetPhiWeight(AliESDtrack* track)
435 Double_t phiweight=1;
437 if (fUsePhiWeight && fPhiDist && track) {
438 Double_t nParticles = fPhiDist->Integral();
439 Double_t nPhibins = fPhiDist->GetNbinsX();
441 Double_t Phi = track->Phi();
443 while (Phi<0) Phi += TMath::TwoPi();
444 while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
446 Double_t PhiDistValue = fPhiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
448 if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
453 //__________________________________________________________________________
454 void AliEPSelectionTask::SetPhiDist()
456 if(!fuserphidist) { // if it's already set and custom class is required, we use the one provided by the user
458 fPhiDist = (TH1F*) fEPContainer->GetObject(frunNumber, "Default");
459 if (!fPhiDist) AliFatal(Form("Cannot find OADB phi distribution for run %d", frunNumber));
463 AliInfo("Using Custom Phi Distribution");
472 for (int i=1; i<fPhiDist->GetNbinsX(); i++){
473 if (!((fPhiDist->GetBinContent(i))>0)) {
478 cout << "empty bins - rebinning!" << endl;
486 AliError("After Maximum of rebinning still empty Phi-bins!!!");
490 //__________________________________________________________________________
491 void AliEPSelectionTask::SetPersonalPhiDistribution(char* infilename, char* listname)
494 fuserphidist = kTRUE;
497 TObject* list = f.Get(listname);
498 fPhiDist = (TH1F*)list->FindObject("fHOutPhi");
499 if (!fPhiDist) AliFatal("Phi Distribution not found!!!");