[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());
	if (headerH) 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));
	  if (track) {
	    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 && (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));
    if (track) {
      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));
    if (!track) continue;
    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;
    return;
  }
  
  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 && fPhiDist && track) {
    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());
b1a983b4	208	if (headerH) fRP = headerH->ReactionPlaneAngle();
ce7adfe9	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));
b1a983b4	242	if (track) {
	243	float delta = track->Phi()-fEventplaneQ;
	244	while (delta < 0) delta += TMath::Pi();
	245	while (delta > TMath::Pi()) delta -= TMath::Pi();
	246	fHOutPTPsi->Fill(track->Pt(),delta);
	247	fHOutPhi->Fill(track->Phi());
	248	fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
	249	}
ce7adfe9	250	}
	251
	252	AliESDtrack* trmax = esd->GetTrack(0);
	253	for (int iter = 1; iter<NT;iter++){
	254	AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter));
b1a983b4	255	if (track && (track->Pt() > trmax->Pt())) trmax = track;
ce7adfe9	256	}
	257	fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);
	258	}
	259	delete ftracklist;
	260	}
	261	}
	262
	263	else if (fAnalysisInput.CompareTo("AOD")==0){
	264	//AliAODEvent aod = dynamic_cast<AliAODEvent> (InputEvent());
	265	// to be implemented
	266	printf(" AOD analysis not yet implemented!!!\n\n");
	267	return;
	268	}
	269	else {
	270	printf(" Analysis Input not known!\n\n ");
	271	return;
	272	}
	273	PostData(1, fOutputList);
	274	}
	275
	276	//________________________________________________________________________
	277	void AliEPSelectionTask::Terminate(Option_t /option*/)
	278	{
	279	// Terminate analysis
	280	}
	281
	282	//__________________________________________________________________________
	283	TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
	284	{
	285	TVector2 mQ;
	286	float mQx=0, mQy=0;
	287	AliESDtrack* track;
	288	Double_t weight;
	289
	290	int NT = tracklist->GetEntries();
	291
	292	for (int i=0; i<NT; i++){
	293	weight = 1;
	294	track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
b1a983b4	295	if (track) {
	296	weight = GetWeight(track);
	297	if (fSaveTrackContribution){
	298	EP->GetQContributionXArray()->AddAt(weightcos(2track->Phi()),track->GetID());
	299	EP->GetQContributionYArray()->AddAt(weightsin(2track->Phi()),track->GetID());
	300	}
	301	mQx += (weightcos(2track->Phi()));
	302	mQy += (weightsin(2track->Phi()));
ce7adfe9	303	}
ce7adfe9	304	}
	305	mQ.Set(mQx,mQy);
	306	return mQ;
	307	}
	308
	309	//________________________________________________________________________
	310	void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist)
	311	{
	312	TVector2 mQ[2];
	313	float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
	314	Double_t weight;
	315
	316	AliESDtrack* track;
	317	TRandom2 RN = 0;
	318
	319	int NT = tracklist->GetEntries();
	320	int trackcounter1=0, trackcounter2=0;
	321
	322	for (Int_t i = 0; i < NT; i++) {
	323	weight = 1;
	324	track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
b1a983b4	325	if (!track) continue;
ce7adfe9	326	weight = GetWeight(track);
	327
	328	// This loop splits the track set into 2 random subsets
	329	if( trackcounter1 < int(NT/2.) && trackcounter2 < int(NT/2.)){
	330	float random = RN.Rndm();
	331	if(random < .5){
	332	mQx1 += (weightcos(2track->Phi()));
	333	mQy1 += (weightsin(2track->Phi()));
	334	trackcounter1++;
	335	}
	336	else {
	337	mQx2 += (weightcos(2track->Phi()));
	338	mQy2 += (weightsin(2track->Phi()));
	339	trackcounter2++;
	340	}
	341	}
	342	else if( trackcounter1 >= int(NT/2.)){
	343	mQx2 += (weightcos(2track->Phi()));
	344	mQy2 += (weightsin(2track->Phi()));
	345	trackcounter2++;
	346	}
	347	else {
	348	mQx1 += (weightcos(2track->Phi()));
	349	mQy1 += (weightsin(2track->Phi()));
	350	trackcounter1++;
	351	}
	352	}
	353	mQ[0].Set(mQx1,mQy1);
	354	mQ[1].Set(mQx2,mQy2);
	355	Q1 = mQ[0];
	356	Q2 = mQ[1];
	357	}
	358
	359	//________________________________________________________________________
	360	void AliEPSelectionTask::SetESDtrackCuts(TString status){
	361
	362	fStatus = status;
	363	if (fStatus.CompareTo("GLOBAL")==0) fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
	364	if (fStatus.CompareTo("TPC")==0) fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
	365	fESDtrackCuts->SetPtRange(0.15,20);
	366	fESDtrackCuts->SetEtaRange(-0.8,0.8);
	367	}
	368
	369	//__________________________________________________________________________
08c34a44	370	void AliEPSelectionTask::SetPhiDistribution(char* infilename, char* listname)
08c34a44	371	{
ce7adfe9	372	TFile f(infilename);
	373	TObject* list = f.Get(listname);
	374	fPhiDist = (TH1F*)list->FindObject("fHOutPhi");
e4f1ed0c	375	if (!fPhiDist) {
	376	cout << "Phi Distribution not found!!!" << endl;
	377	return;
	378	}
ce7adfe9	379
	380	Bool_t emptybins;
	381
	382	int iter = 0;
	383	while (iter<3){
	384	emptybins = kFALSE;
	385
	386	for (int i=1; i<fPhiDist->GetNbinsX(); i++){
	387	if (!((fPhiDist->GetBinContent(i))>0)) {
	388	emptybins = kTRUE;
	389	}
	390	}
	391	if (emptybins) {
	392	cout << "empty bins - rebinning!" << endl;
	393	fPhiDist->Rebin();
	394	iter++;
	395	}
	396	else iter = 3;
	397	}
	398
	399	if (emptybins) {
	400	AliError("After Maximum of rebinning still empty Phi-bins!!!");
	401	}
	402	f.Close();
	403	}
	404
	405	//________________________________________________________________________
	406	Double_t AliEPSelectionTask::GetWeight(AliESDtrack* track)
	407	{
	408	Double_t ptweight=1;
	409
	410	if (fUsePtWeight) {
	411	if (track->Pt()<2) ptweight=track->Pt();
	412	else ptweight=2;
	413	}
	414	return ptweight*GetPhiWeight(track);
	415	}
	416
	417	//________________________________________________________________________
	418	Double_t AliEPSelectionTask::GetPhiWeight(AliESDtrack* track)
	419	{
	420	Double_t phiweight=1;
	421
e4f1ed0c	422	if (fUsePhiWeight && fPhiDist && track) {
ce7adfe9	423	Double_t nParticles = fPhiDist->Integral();
	424	Double_t nPhibins = fPhiDist->GetNbinsX();
	425
	426	Double_t Phi = track->Phi();
	427
	428	while (Phi<0) Phi += TMath::TwoPi();
	429	while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
	430
08c34a44	431	Double_t PhiDistValue = fPhiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
ce7adfe9	432
	433	if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
	434	}
	435	return phiweight;
	436	}