[u/mrichter/AliRoot.git] / ANALYSIS / AliEPSelectionTask.cxx

/**************************************************************************
 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//*****************************************************
//   Class AliEPSelectionTask
//   Class to determine event plane            
//   author: Alberica Toia, Johanna Gramling
//*****************************************************

#include "AliEPSelectionTask.h"

#include <TTree.h>
#include <TList.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TProfile.h>
#include <TFile.h>
#include <TObjString.h>
#include <TString.h>
#include <TCanvas.h>
#include <TROOT.h>
#include <TDirectory.h>
#include <TSystem.h>
#include <iostream>
#include <TRandom2.h>
#include <TArrayF.h>

#include "AliAnalysisManager.h"
#include "AliVEvent.h"
#include "AliESD.h"
#include "AliESDEvent.h"
#include "AliMCEvent.h"
#include "AliESDtrack.h"
#include "AliESDtrackCuts.h"
#include "AliESDHeader.h"
#include "AliESDInputHandler.h"
#include "AliAODHandler.h"
#include "AliAODEvent.h"
#include "AliHeader.h"
#include "AliGenHijingEventHeader.h"
#include "AliAnalysisTaskSE.h"
#include "AliPhysicsSelectionTask.h"
#include "AliPhysicsSelection.h"
#include "AliBackgroundSelection.h"
#include "AliESDUtils.h"

#include "AliEventplane.h"

ClassImp(AliEPSelectionTask)

//________________________________________________________________________
AliEPSelectionTask::AliEPSelectionTask():
AliAnalysisTaskSE(),
  fDebug(0),
  fAnalysisInput("ESD"),
  fStatus("TPC"),
  fUseMCRP(kFALSE),
  fUsePhiWeight(kFALSE),
  fUsePtWeight(kFALSE),
  fSaveTrackContribution(kFALSE),
  fESDtrackCuts(0),
  ftracklist(0),
  fPhiDist(0),
  fQVector(0),
  fQContributionX(0),
  fQContributionY(0),
  fEventplaneQ(0),
  fQsub1(0),
  fQsub2(0),
  fQsubRes(0),  
  fOutputList(0),
  fHOutEventplaneQ(0),
  fHOutPhi(0),
  fHOutPhiCorr(0),
  fHOutsub1sub2(0),
  fHOutNTEPRes(0),
  fHOutPTPsi(0),
  fHOutDiff(0),
  fHOutleadPTPsi(0)
{   
  // Default constructor
  AliInfo("Event Plane Selection enabled.");
}   

//________________________________________________________________________
AliEPSelectionTask::AliEPSelectionTask(const char *name):
  AliAnalysisTaskSE(name),
  fDebug(0),
  fAnalysisInput("ESD"),
  fStatus("TPC"),
  fUseMCRP(kFALSE),
  fUsePhiWeight(kFALSE),
  fUsePtWeight(kFALSE),  
  fSaveTrackContribution(kFALSE),
  fESDtrackCuts(0),
  ftracklist(0),
  fPhiDist(0),
  fQVector(0),
  fQContributionX(0),
  fQContributionY(0),
  fEventplaneQ(0),
  fQsub1(0),
  fQsub2(0),
  fQsubRes(0),
  fOutputList(0),
  fHOutEventplaneQ(0),
  fHOutPhi(0),
  fHOutPhiCorr(0),
  fHOutsub1sub2(0),
  fHOutNTEPRes(0),
  fHOutPTPsi(0),
  fHOutDiff(0),
  fHOutleadPTPsi(0)
{
  // Default constructor
  AliInfo("Event Plane Selection enabled.");
  DefineOutput(1, TList::Class());
}
 
//________________________________________________________________________
AliEPSelectionTask::~AliEPSelectionTask()
{
  // Destructor  
  if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()){
      delete fOutputList;
      fOutputList = 0;
  }
  if (fESDtrackCuts){
      delete fESDtrackCuts;
      fESDtrackCuts = 0;
  }
  if (fQVector){
      delete fQVector;
      fQVector = 0;
  }
  if (fQsub1){
      delete fQsub1;
      fQsub1 = 0;
  }
    if (fQsub2){
      delete fQsub2;
      fQsub2 = 0;
  }
}  

//________________________________________________________________________
void AliEPSelectionTask::UserCreateOutputObjects()
{  
  // Create the output containers
  if (fDebug>1) printf("AliEPSelectionTask::UserCreateOutputObjects() \n");
  AliLog::SetClassDebugLevel("AliEPSelectionTask", AliLog::kInfo);

  fOutputList = new TList();
  fOutputList->SetOwner();
  fHOutEventplaneQ = new TH1F("fHOutEventplaneQ","fHOutEventplaneQ; Event Plane Q",100,0,TMath::Pi());
  fHOutPhi = new TH1F("fHOutPhi","fHOutPhi; Phi Distribution",100,0,TMath::TwoPi());
  fHOutPhiCorr = new TH1F("fHOutPhiCorr","fHOutPhiCorr; Corrected Phi Distribution",100,0,TMath::TwoPi());
  fHOutsub1sub2 = new TH2F("fHOutsub1sub2","fHOutsub1sub2; EP1; EP2",100,0,TMath::Pi(),100,0,TMath::Pi());
  fHOutNTEPRes = new TH2F("fHOutNTEPRes","fHOutNTEPRes; Number of Tracks; Event Plane Resolution",100,0,5000,100,-TMath::Pi(),TMath::Pi());
  fHOutPTPsi = new TH2F("fHOutPTPsi","fHOutPTPsi; PT; Phi-EP",100,0,20,100,0,TMath::Pi());
  fHOutDiff = new TH2F("fHOutDiff","fHOutDiff; EP; MCEP",100,0,TMath::Pi(),100,0,TMath::Pi());
  fHOutleadPTPsi = new TH2F("fHOutleadPTPsi","fHOutleadPTPsi; leadPT; EP",100,0,TMath::Pi(),100,0,TMath::Pi());

  fOutputList->Add(fHOutEventplaneQ);
  fOutputList->Add(fHOutPhi);
  fOutputList->Add(fHOutPhiCorr);
  fOutputList->Add(fHOutsub1sub2);
  fOutputList->Add(fHOutNTEPRes);
  fOutputList->Add(fHOutPTPsi);
  fOutputList->Add(fHOutleadPTPsi);
  fOutputList->Add(fHOutDiff);
  
  PostData(1, fOutputList); 
}

//________________________________________________________________________
void AliEPSelectionTask::UserExec(Option_t */*option*/)
{ 
  // Execute analysis for current event:
  if (fDebug>1) printf(" **** AliEPSelectionTask::UserExec() \n");
 
  AliEventplane* esdEP = 0;
  TVector2 QQ1;
  TVector2 QQ2;
  Double_t fRP = 0.; // the monte carlo reaction plane angle
    
  if (fAnalysisInput.CompareTo("ESD")==0){
    
    AliVEvent* event = InputEvent();
    AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
    
    if (fUseMCRP) {
      AliMCEvent* mcEvent  = MCEvent();      
      if (mcEvent && mcEvent->GenEventHeader()) {
	AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
	fRP = headerH->ReactionPlaneAngle();
      }
    }
    if (esd){
      esdEP = esd->GetEventplane();
      if (fSaveTrackContribution) {
	esdEP->GetQContributionXArray()->Set(esd->GetNumberOfTracks());
	esdEP->GetQContributionYArray()->Set(esd->GetNumberOfTracks());
      }
      
      if (fStatus.CompareTo("GLOBAL")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE);
      if (fStatus.CompareTo("TPC")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE);
      const int NT = ftracklist->GetEntries();
      
      if (NT>4){
	fQVector = new TVector2(GetQ(esdEP,ftracklist));
	fEventplaneQ = fQVector->Phi()/2; 
	GetQsub(QQ1, QQ2, ftracklist);
	fQsub1 = new TVector2(QQ1);
	fQsub2 = new TVector2(QQ2);
	fQsubRes = (fQsub1->Phi()/2 - fQsub2->Phi()/2);
	esdEP->SetQVector(fQVector);
	esdEP->SetEventplaneQ(fEventplaneQ);
	esdEP->SetQsub(fQsub1,fQsub2);
	esdEP->SetQsubRes(fQsubRes);
	
	fHOutEventplaneQ->Fill(fEventplaneQ);
	fHOutsub1sub2->Fill(fQsub1->Phi()/2,fQsub2->Phi()/2);
	fHOutNTEPRes->Fill(NT,fQsubRes);

	if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
	
	for (int iter = 0; iter<NT;iter++){
	  AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter));
	  float delta = track->Phi()-fEventplaneQ;
	  while (delta < 0) delta += TMath::Pi();
	  while (delta > TMath::Pi()) delta -= TMath::Pi();
	  fHOutPTPsi->Fill(track->Pt(),delta);
	  fHOutPhi->Fill(track->Phi());
	  fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
	}
	
	AliESDtrack* trmax = esd->GetTrack(0);
	for (int iter = 1; iter<NT;iter++){
	  AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter));
	  if (track->Pt() > trmax->Pt()) trmax = track;
	}
	fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);      
      }
      delete ftracklist;
    }
  }
  
  else if (fAnalysisInput.CompareTo("AOD")==0){
    //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
    // to be implemented
    printf("  AOD analysis not yet implemented!!!\n\n");
    return;
  }  
  else {
    printf(" Analysis Input not known!\n\n ");
    return;
  }
  PostData(1, fOutputList); 
}

//________________________________________________________________________
void AliEPSelectionTask::Terminate(Option_t */*option*/)
{
  // Terminate analysis
}

//__________________________________________________________________________
TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
{
  TVector2 mQ;
  float mQx=0, mQy=0;
  AliESDtrack* track;
  Double_t weight;
  
  int NT = tracklist->GetEntries();

  for (int i=0; i<NT; i++){
    weight = 1;
    track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
    weight = GetWeight(track);
    if (fSaveTrackContribution){
      EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi()),track->GetID());
      EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi()),track->GetID());
    }
    mQx += (weight*cos(2*track->Phi()));
    mQy += (weight*sin(2*track->Phi()));
  }
  mQ.Set(mQx,mQy);
  return mQ;
}
  
  //________________________________________________________________________
void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist)
{
  TVector2 mQ[2];
  float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
  Double_t weight;

  AliESDtrack* track;
  TRandom2 RN = 0;
  
  int NT = tracklist->GetEntries();
  int trackcounter1=0, trackcounter2=0;
  
  for (Int_t i = 0; i < NT; i++) {
    weight = 1;
    track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
    weight = GetWeight(track);
    
    // This loop splits the track set into 2 random subsets
    if( trackcounter1 < int(NT/2.) && trackcounter2 < int(NT/2.)){
      float random = RN.Rndm();
      if(random < .5){
        mQx1 += (weight*cos(2*track->Phi()));
        mQy1 += (weight*sin(2*track->Phi()));
        trackcounter1++;
      }
      else {
        mQx2 += (weight*cos(2*track->Phi()));
        mQy2 += (weight*sin(2*track->Phi()));
        trackcounter2++;
      }
    }
    else if( trackcounter1 >= int(NT/2.)){
      mQx2 += (weight*cos(2*track->Phi()));
      mQy2 += (weight*sin(2*track->Phi()));
      trackcounter2++;
    }
    else {
      mQx1 += (weight*cos(2*track->Phi()));
      mQy1 += (weight*sin(2*track->Phi()));
      trackcounter1++;
    }
  }
  mQ[0].Set(mQx1,mQy1);
  mQ[1].Set(mQx2,mQy2);
  Q1 = mQ[0];
  Q2 = mQ[1];
}

//________________________________________________________________________
void AliEPSelectionTask::SetESDtrackCuts(TString status){
  
  fStatus = status;
  if (fStatus.CompareTo("GLOBAL")==0) fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
  if (fStatus.CompareTo("TPC")==0)    fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
  fESDtrackCuts->SetPtRange(0.15,20);
  fESDtrackCuts->SetEtaRange(-0.8,0.8);
}

//__________________________________________________________________________
void AliEPSelectionTask::SetPhiDistribution(char* infilename, char* listname)
{
  TFile f(infilename);
  TObject* list = f.Get(listname);
  fPhiDist = (TH1F*)list->FindObject("fHOutPhi");
  if (!fPhiDist) cout << "Phi Distribution not found!!!" << endl;
  
  Bool_t emptybins;

  int iter = 0;  
  while (iter<3){
      emptybins = kFALSE;
   
      for (int i=1; i<fPhiDist->GetNbinsX(); i++){
	if (!((fPhiDist->GetBinContent(i))>0)) {
	  emptybins = kTRUE;
	}
      }  
      if (emptybins) {
	cout << "empty bins - rebinning!" << endl;
	fPhiDist->Rebin();
	iter++;
      }      
      else iter = 3;
  }
  
  if (emptybins) {
    AliError("After Maximum of rebinning still empty Phi-bins!!!");
  }
  f.Close();
}

//________________________________________________________________________
Double_t AliEPSelectionTask::GetWeight(AliESDtrack* track)
{
  Double_t ptweight=1;

  if (fUsePtWeight) {      
    if (track->Pt()<2) ptweight=track->Pt();
    else ptweight=2;
  }
  return ptweight*GetPhiWeight(track);
}

//________________________________________________________________________
Double_t AliEPSelectionTask::GetPhiWeight(AliESDtrack* track)
{
  Double_t phiweight=1;
  
  if (fUsePhiWeight) {
    Double_t nParticles = fPhiDist->Integral();
    Double_t nPhibins = fPhiDist->GetNbinsX();
  
    Double_t Phi = track->Phi();
    
    while (Phi<0) Phi += TMath::TwoPi();
    while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
      
    Double_t PhiDistValue = fPhiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
    
    if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
  }
  return phiweight;
}
Commit	Line	Data
ce7adfe9	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
	59	#include "AliEventplane.h"
	60
	61	ClassImp(AliEPSelectionTask)
	62
	63	//________________________________________________________________________
	64	AliEPSelectionTask::AliEPSelectionTask():
65	AliAnalysisTaskSE(),
66	fDebug(0),
67	fAnalysisInput("ESD"),
68	fStatus("TPC"),
69	fUseMCRP(kFALSE),
70	fUsePhiWeight(kFALSE),
71	fUsePtWeight(kFALSE),
72	fSaveTrackContribution(kFALSE),
73	fESDtrackCuts(0),
74	ftracklist(0),
75	fPhiDist(0),
76	fQVector(0),
77	fQContributionX(0),
78	fQContributionY(0),
79	fEventplaneQ(0),
80	fQsub1(0),
81	fQsub2(0),
82	fQsubRes(0),
83	fOutputList(0),
84	fHOutEventplaneQ(0),
85	fHOutPhi(0),
86	fHOutPhiCorr(0),
87	fHOutsub1sub2(0),
88	fHOutNTEPRes(0),
89	fHOutPTPsi(0),
90	fHOutDiff(0),
91	fHOutleadPTPsi(0)
92	{
93	// Default constructor
94	AliInfo("Event Plane Selection enabled.");
95	}
96
97	//________________________________________________________________________
98	AliEPSelectionTask::AliEPSelectionTask(const char *name):
99	AliAnalysisTaskSE(name),
100	fDebug(0),
101	fAnalysisInput("ESD"),
102	fStatus("TPC"),
103	fUseMCRP(kFALSE),
104	fUsePhiWeight(kFALSE),
105	fUsePtWeight(kFALSE),
106	fSaveTrackContribution(kFALSE),
107	fESDtrackCuts(0),
108	ftracklist(0),
109	fPhiDist(0),
110	fQVector(0),
111	fQContributionX(0),
112	fQContributionY(0),
113	fEventplaneQ(0),
114	fQsub1(0),
115	fQsub2(0),
116	fQsubRes(0),
117	fOutputList(0),
118	fHOutEventplaneQ(0),
119	fHOutPhi(0),
120	fHOutPhiCorr(0),
121	fHOutsub1sub2(0),
122	fHOutNTEPRes(0),
123	fHOutPTPsi(0),
124	fHOutDiff(0),
125	fHOutleadPTPsi(0)
126	{
127	// Default constructor
128	AliInfo("Event Plane Selection enabled.");
129	DefineOutput(1, TList::Class());
130	}
131
132	//________________________________________________________________________
133	AliEPSelectionTask::~AliEPSelectionTask()
134	{
135	// Destructor
136	if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()){
137	delete fOutputList;
138	fOutputList = 0;
139	}
140	if (fESDtrackCuts){
141	delete fESDtrackCuts;
142	fESDtrackCuts = 0;
143	}
144	if (fQVector){
145	delete fQVector;
146	fQVector = 0;
147	}
148	if (fQsub1){
149	delete fQsub1;
150	fQsub1 = 0;
151	}
152	if (fQsub2){
153	delete fQsub2;
154	fQsub2 = 0;
155	}
156	}
157
158	//________________________________________________________________________
159	void AliEPSelectionTask::UserCreateOutputObjects()
160	{
161	// Create the output containers
162	if (fDebug>1) printf("AliEPSelectionTask::UserCreateOutputObjects() \n");
163	AliLog::SetClassDebugLevel("AliEPSelectionTask", AliLog::kInfo);
164
165	fOutputList = new TList();
166	fOutputList->SetOwner();
167	fHOutEventplaneQ = new TH1F("fHOutEventplaneQ","fHOutEventplaneQ; Event Plane Q",100,0,TMath::Pi());
168	fHOutPhi = new TH1F("fHOutPhi","fHOutPhi; Phi Distribution",100,0,TMath::TwoPi());
169	fHOutPhiCorr = new TH1F("fHOutPhiCorr","fHOutPhiCorr; Corrected Phi Distribution",100,0,TMath::TwoPi());
170	fHOutsub1sub2 = new TH2F("fHOutsub1sub2","fHOutsub1sub2; EP1; EP2",100,0,TMath::Pi(),100,0,TMath::Pi());
171	fHOutNTEPRes = new TH2F("fHOutNTEPRes","fHOutNTEPRes; Number of Tracks; Event Plane Resolution",100,0,5000,100,-TMath::Pi(),TMath::Pi());
172	fHOutPTPsi = new TH2F("fHOutPTPsi","fHOutPTPsi; PT; Phi-EP",100,0,20,100,0,TMath::Pi());
173	fHOutDiff = new TH2F("fHOutDiff","fHOutDiff; EP; MCEP",100,0,TMath::Pi(),100,0,TMath::Pi());
174	fHOutleadPTPsi = new TH2F("fHOutleadPTPsi","fHOutleadPTPsi; leadPT; EP",100,0,TMath::Pi(),100,0,TMath::Pi());
175
176	fOutputList->Add(fHOutEventplaneQ);
177	fOutputList->Add(fHOutPhi);
178	fOutputList->Add(fHOutPhiCorr);
179	fOutputList->Add(fHOutsub1sub2);
180	fOutputList->Add(fHOutNTEPRes);
181	fOutputList->Add(fHOutPTPsi);
182	fOutputList->Add(fHOutleadPTPsi);
183	fOutputList->Add(fHOutDiff);
184
185	PostData(1, fOutputList);
186	}
187
188	//________________________________________________________________________
189	void AliEPSelectionTask::UserExec(Option_t /option*/)
190	{
191	// Execute analysis for current event:
192	if (fDebug>1) printf(" **** AliEPSelectionTask::UserExec() \n");
193
194	AliEventplane* esdEP = 0;
195	TVector2 QQ1;
196	TVector2 QQ2;
197	Double_t fRP = 0.; // the monte carlo reaction plane angle
198
199	if (fAnalysisInput.CompareTo("ESD")==0){
200
201	AliVEvent* event = InputEvent();
202	AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
203
204	if (fUseMCRP) {
205	AliMCEvent* mcEvent = MCEvent();
206	if (mcEvent && mcEvent->GenEventHeader()) {
207	AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
208	fRP = headerH->ReactionPlaneAngle();
209	}
210	}
211	if (esd){
212	esdEP = esd->GetEventplane();
213	if (fSaveTrackContribution) {
214	esdEP->GetQContributionXArray()->Set(esd->GetNumberOfTracks());
215	esdEP->GetQContributionYArray()->Set(esd->GetNumberOfTracks());
216	}
217
218	if (fStatus.CompareTo("GLOBAL")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE);
219	if (fStatus.CompareTo("TPC")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE);
220	const int NT = ftracklist->GetEntries();
221
222	if (NT>4){
223	fQVector = new TVector2(GetQ(esdEP,ftracklist));
224	fEventplaneQ = fQVector->Phi()/2;
225	GetQsub(QQ1, QQ2, ftracklist);
226	fQsub1 = new TVector2(QQ1);
227	fQsub2 = new TVector2(QQ2);
228	fQsubRes = (fQsub1->Phi()/2 - fQsub2->Phi()/2);
ce7adfe9	229	esdEP->SetQVector(fQVector);
	230	esdEP->SetEventplaneQ(fEventplaneQ);
	231	esdEP->SetQsub(fQsub1,fQsub2);
	232	esdEP->SetQsubRes(fQsubRes);
	233
	234	fHOutEventplaneQ->Fill(fEventplaneQ);
	235	fHOutsub1sub2->Fill(fQsub1->Phi()/2,fQsub2->Phi()/2);
	236	fHOutNTEPRes->Fill(NT,fQsubRes);
	237
	238	if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
	239
	240	for (int iter = 0; iter<NT;iter++){
	241	AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter));
	242	float delta = track->Phi()-fEventplaneQ;
	243	while (delta < 0) delta += TMath::Pi();
	244	while (delta > TMath::Pi()) delta -= TMath::Pi();
	245	fHOutPTPsi->Fill(track->Pt(),delta);
	246	fHOutPhi->Fill(track->Phi());
	247	fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
	248	}
	249
	250	AliESDtrack* trmax = esd->GetTrack(0);
	251	for (int iter = 1; iter<NT;iter++){
	252	AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter));
	253	if (track->Pt() > trmax->Pt()) trmax = track;
	254	}
	255	fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);
	256	}
	257	delete ftracklist;
	258	}
	259	}
	260
	261	else if (fAnalysisInput.CompareTo("AOD")==0){
	262	//AliAODEvent aod = dynamic_cast<AliAODEvent> (InputEvent());
	263	// to be implemented
	264	printf(" AOD analysis not yet implemented!!!\n\n");
	265	return;
	266	}
	267	else {
	268	printf(" Analysis Input not known!\n\n ");
	269	return;
	270	}
	271	PostData(1, fOutputList);
	272	}
	273
	274	//________________________________________________________________________
	275	void AliEPSelectionTask::Terminate(Option_t /option*/)
	276	{
	277	// Terminate analysis
	278	}
	279
	280	//__________________________________________________________________________
	281	TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
	282	{
	283	TVector2 mQ;
	284	float mQx=0, mQy=0;
	285	AliESDtrack* track;
	286	Double_t weight;
	287
	288	int NT = tracklist->GetEntries();
	289
	290	for (int i=0; i<NT; i++){
	291	weight = 1;
	292	track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
293	weight = GetWeight(track);
294	if (fSaveTrackContribution){
295	EP->GetQContributionXArray()->AddAt(weightcos(2track->Phi()),track->GetID());
296	EP->GetQContributionYArray()->AddAt(weightsin(2track->Phi()),track->GetID());
297	}
298	mQx += (weightcos(2track->Phi()));
299	mQy += (weightsin(2track->Phi()));
300	}
301	mQ.Set(mQx,mQy);
302	return mQ;
303	}
304
305	//________________________________________________________________________
306	void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist)
307	{
308	TVector2 mQ[2];
309	float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
310	Double_t weight;
311
312	AliESDtrack* track;
313	TRandom2 RN = 0;
314
315	int NT = tracklist->GetEntries();
316	int trackcounter1=0, trackcounter2=0;
317
318	for (Int_t i = 0; i < NT; i++) {
319	weight = 1;
320	track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
321	weight = GetWeight(track);
322
323	// This loop splits the track set into 2 random subsets
324	if( trackcounter1 < int(NT/2.) && trackcounter2 < int(NT/2.)){
325	float random = RN.Rndm();
326	if(random < .5){
327	mQx1 += (weightcos(2track->Phi()));
328	mQy1 += (weightsin(2track->Phi()));
329	trackcounter1++;
330	}
331	else {
332	mQx2 += (weightcos(2track->Phi()));
333	mQy2 += (weightsin(2track->Phi()));
334	trackcounter2++;
335	}
336	}
337	else if( trackcounter1 >= int(NT/2.)){
338	mQx2 += (weightcos(2track->Phi()));
339	mQy2 += (weightsin(2track->Phi()));
340	trackcounter2++;
341	}
342	else {
343	mQx1 += (weightcos(2track->Phi()));
344	mQy1 += (weightsin(2track->Phi()));
345	trackcounter1++;
346	}
347	}
348	mQ[0].Set(mQx1,mQy1);
349	mQ[1].Set(mQx2,mQy2);
350	Q1 = mQ[0];
351	Q2 = mQ[1];
352	}
353
354	//________________________________________________________________________
355	void AliEPSelectionTask::SetESDtrackCuts(TString status){
356
357	fStatus = status;
358	if (fStatus.CompareTo("GLOBAL")==0) fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
359	if (fStatus.CompareTo("TPC")==0) fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
360	fESDtrackCuts->SetPtRange(0.15,20);
361	fESDtrackCuts->SetEtaRange(-0.8,0.8);
362	}
363
364	//__________________________________________________________________________
08c34a44	365	void AliEPSelectionTask::SetPhiDistribution(char* infilename, char* listname)
08c34a44	366	{
ce7adfe9	367	TFile f(infilename);
	368	TObject* list = f.Get(listname);
	369	fPhiDist = (TH1F*)list->FindObject("fHOutPhi");
	370	if (!fPhiDist) cout << "Phi Distribution not found!!!" << endl;
	371
	372	Bool_t emptybins;
	373
	374	int iter = 0;
	375	while (iter<3){
	376	emptybins = kFALSE;
	377
	378	for (int i=1; i<fPhiDist->GetNbinsX(); i++){
	379	if (!((fPhiDist->GetBinContent(i))>0)) {
	380	emptybins = kTRUE;
	381	}
	382	}
	383	if (emptybins) {
	384	cout << "empty bins - rebinning!" << endl;
	385	fPhiDist->Rebin();
	386	iter++;
	387	}
	388	else iter = 3;
	389	}
	390
	391	if (emptybins) {
	392	AliError("After Maximum of rebinning still empty Phi-bins!!!");
	393	}
	394	f.Close();
	395	}
	396
	397	//________________________________________________________________________
	398	Double_t AliEPSelectionTask::GetWeight(AliESDtrack* track)
	399	{
	400	Double_t ptweight=1;
	401
	402	if (fUsePtWeight) {
	403	if (track->Pt()<2) ptweight=track->Pt();
	404	else ptweight=2;
	405	}
	406	return ptweight*GetPhiWeight(track);
	407	}
	408
	409	//________________________________________________________________________
	410	Double_t AliEPSelectionTask::GetPhiWeight(AliESDtrack* track)
	411	{
	412	Double_t phiweight=1;
	413
	414	if (fUsePhiWeight) {
	415	Double_t nParticles = fPhiDist->Integral();
	416	Double_t nPhibins = fPhiDist->GetNbinsX();
	417
	418	Double_t Phi = track->Phi();
	419
	420	while (Phi<0) Phi += TMath::TwoPi();
	421	while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
	422
08c34a44	423	Double_t PhiDistValue = fPhiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
ce7adfe9	424
	425	if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
	426	}
	427	return phiweight;
	428	}