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