[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 "AliOADBContainer.h"

#include "AliEventplane.h"

ClassImp(AliEPSelectionTask)

//________________________________________________________________________
AliEPSelectionTask::AliEPSelectionTask():
AliAnalysisTaskSE(),
  fAnalysisInput("ESD"),
  fTrackType("TPC"),
  fUseMCRP(kFALSE),
  fUsePhiWeight(kFALSE),
  fUsePtWeight(kFALSE),
  fSaveTrackContribution(kFALSE),
  fUserphidist(kFALSE),
  fUsercuts(kFALSE),
  fRunNumber(-15),
  fESDtrackCuts(0),
  fEPContainer(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),
  fAnalysisInput("ESD"),
  fTrackType("TPC"),
  fUseMCRP(kFALSE),
  fUsePhiWeight(kFALSE),
  fUsePtWeight(kFALSE),  
  fSaveTrackContribution(kFALSE),
  fUserphidist(kFALSE),
  fUsercuts(kFALSE),
  fRunNumber(-15),
  fESDtrackCuts(0),
  fEPContainer(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 (fUserphidist) {
    if (fPhiDist) {
      delete fPhiDist;
      fPhiDist = 0;
    }
  }
  if (fEPContainer){
      delete fEPContainer;
      fEPContainer = 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); 
  
    if(!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user

    TString oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath()));

    TFile foadb(oadbfilename); 
    if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));

    AliInfo("Using Standard OADB");
    fEPContainer = (AliOADBContainer*) foadb.Get("epphidist");    
    if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection");
    foadb.Close();
    }
}

//________________________________________________________________________
void AliEPSelectionTask::UserExec(Option_t */*option*/)
{ 
  // Execute analysis for current event:
  if (fDebug>1) printf(" **** AliEPSelectionTask::UserExec() \n");
  
//   fRunNumber = -15;
 
  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 (esd){    
      if (!(fRunNumber == esd->GetRunNumber())) {
	  fRunNumber = esd->GetRunNumber();
	    SetPhiDist();      
      }
      
      
      if (fUseMCRP) {
	  AliMCEvent* mcEvent  = MCEvent();      
	  if (mcEvent && mcEvent->GenEventHeader()) {
	      AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
	      if (headerH) fRP = headerH->ReactionPlaneAngle();
	  }
      }
      
      esdEP = esd->GetEventplane();
      if (fSaveTrackContribution) {
	esdEP->GetQContributionXArray()->Set(esd->GetNumberOfTracks());
	esdEP->GetQContributionYArray()->Set(esd->GetNumberOfTracks());
      }
      
      TObjArray* tracklist = new TObjArray;
      if (fTrackType.CompareTo("GLOBAL")==0) tracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE);
      if (fTrackType.CompareTo("TPC")==0) tracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE);
      const int nt = tracklist->GetEntries();
      
      if (nt>4){
	fQVector = new TVector2(GetQ(esdEP,tracklist));
	fEventplaneQ = fQVector->Phi()/2; 
	GetQsub(qq1, qq2, tracklist);
	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*> (tracklist->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*> (tracklist->At(iter));
	  if (track && (track->Pt() > trmax->Pt())) trmax = track;
	}
	fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);      
      }
      delete tracklist;
      tracklist = 0;
    }
  }
  
  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)
{
// Get the Q vector
  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)
{
// Get Qsub
  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::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
  
  if(fESDtrackCuts){ 
    delete fESDtrackCuts;
    fESDtrackCuts = 0;
  }
    
  fUsercuts = kTRUE;
  fESDtrackCuts = trackcuts;
}

//_____________________________________________________________________________
void AliEPSelectionTask::SetTrackType(TString tracktype){
// Set the track type
  fTrackType = tracktype;
  if (!fUsercuts) {
  if (fTrackType.CompareTo("GLOBAL")==0) fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
  if (fTrackType.CompareTo("TPC")==0)    fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
  fESDtrackCuts->SetPtRange(0.15,20);
  fESDtrackCuts->SetEtaRange(-0.8,0.8);
  }
}

//________________________________________________________________________
Double_t AliEPSelectionTask::GetWeight(AliESDtrack* track)
{
// Get weight for 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)
{
// Get phi weight for 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;
}

//__________________________________________________________________________
void AliEPSelectionTask::SetPhiDist() 
{
// Set the phi distribution
  if(!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user

    fPhiDist = (TH1F*) fEPContainer->GetObject(fRunNumber, "Default");
    if (!fPhiDist) AliFatal(Form("Cannot find OADB phi distribution for run %d", fRunNumber));

  } 
  else {
    AliInfo("Using Custom Phi Distribution");
  }
    
  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!!!");
  }
}

//__________________________________________________________________________
void AliEPSelectionTask::SetPersonalPhiDistribution(const char* infilename, char* listname)
{
    // Set a personal phi distribution
  fUserphidist = kTRUE;
  
  TFile f(infilename);
  TObject* list = f.Get(listname);
  fPhiDist = (TH1F*)list->FindObject("fHOutPhi");
  if (!fPhiDist) AliFatal("Phi Distribution not found!!!");

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