[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(),
  fDebug(0),
  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),
  fDebug(0),
  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 (fQVector){
      delete fQVector;
      fQVector = 0;
  }
  if (fQsub1){
      delete fQsub1;
      fQsub1 = 0;
  }
  if (fQsub2){
      delete fQsub2;
      fQsub2 = 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 (!(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();
      }
    }
    if (esd){
      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)
{
  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::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
  
  fusercuts = kTRUE;
  fESDtrackCuts = trackcuts;
}

//_____________________________________________________________________________
void AliEPSelectionTask::SetTrackType(TString tracktype){
  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)
{
  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;
}

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