[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaChargedParticles.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, 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.                  *
 **************************************************************************/
/* $Id: $ */

//_________________________________________________________________________
//
// Class for track selection and identification (not done now)
// Tracks from the CTS are kept in the AOD.
// Few histograms produced.
//
//-- Author: Gustavo Conesa (INFN-LNF)
//_________________________________________________________________________


// --- ROOT system ---
#include "TParticle.h"
#include "TH2F.h"
#include "TH3D.h"
//---- AliRoot system ----
#include "AliAnaChargedParticles.h"
#include "AliCaloTrackReader.h"
#include "AliAODPWG4Particle.h"
#include "AliStack.h"
#include "AliFiducialCut.h"
#include "AliAODTrack.h"
#include "AliAODMCParticle.h"

ClassImp(AliAnaChargedParticles)
  
//____________________________________________________________________________
  AliAnaChargedParticles::AliAnaChargedParticles() : 
    AliAnaPartCorrBaseClass(),fPdg(0), fhNtracks(0),fhVertex(0), fhPt(0),fhPhi(0),fhEta(0), 
    fhPtEtaPhiPos(0), fhPtEtaPhiNeg(0),
    fhPtPion(0),fhPhiPion(0),fhEtaPion(0),
    fhPtProton(0),fhPhiProton(0),fhEtaProton(0),
    fhPtElectron(0),fhPhiElectron(0),fhEtaElectron(0),
    fhPtKaon(0),fhPhiKaon(0),fhEtaKaon(0),
    fhPtUnknown(0),fhPhiUnknown(0),fhEtaUnknown(0)
      
{
  //Default Ctor

  //Initialize parameters
  InitParameters();
}

/*
//____________________________________________________________________________
AliAnaChargedParticles::AliAnaChargedParticles(const AliAnaChargedParticles & ch) :   
  AliAnaPartCorrBaseClass(ch), fPdg(ch.fPdg),  fhPt(ch.fhPt),  fhPhi(ch.fhPhi),fhEta(ch.fhEta), 
  fhPtPion(ch. fhPtPion),fhPhiPion(ch.fhPhiPion),fhEtaPion(ch.fhEtaPion),
  fhPtProton(ch.fhPtProton),fhPhiProton(ch.fhPhiProton),fhEtaProton(ch.fhEtaProton),
  fhPtElectron(ch. fhPtElectron),fhPhiElectron(ch.fhPhiElectron),fhEtaElectron(ch.fhEtaElectron),
  fhPtKaon(ch. fhPtKaon),fhPhiKaon(ch.fhPhiKaon),fhEtaKaon(ch.fhEtaKaon),
  fhPtUnknown(ch.fhPtUnknown),fhPhiUnknown(ch.fhPhiUnknown),fhEtaUnknown(ch.fhEtaUnknown)
  
{
  // cpy ctor
  
}

//_________________________________________________________________________
AliAnaChargedParticles & AliAnaChargedParticles::operator = (const AliAnaChargedParticles & ch)
{
  // assignment operator
  
  if(this == &ch)return *this;
  ((AliAnaPartCorrBaseClass *)this)->operator=(ch);
 
  fPdg = ch.fPdg;
  fhPt = ch.fhPt;
  fhPhi = ch.fhPhi;
  fhEta = ch.fhEta;
  
  fhPtPion = ch. fhPtPion; fhPhiPion = ch.fhPhiPion; fhEtaPion = ch.fhEtaPion; 
  fhPtProton = ch.fhPtProton; fhPhiProton = ch.fhPhiProton; fhEtaProton = ch.fhEtaProton; 
  fhPtElectron = ch. fhPtElectron; fhPhiElectron = ch.fhPhiElectron; fhEtaElectron = ch.fhEtaElectron; 
  fhPtKaon = ch. fhPtKaon; fhPhiKaon = ch.fhPhiKaon; fhEtaKaon = ch.fhEtaKaon; 
  fhPtUnknown = ch.fhPtUnknown; fhPhiUnknown = ch.fhPhiUnknown; fhEtaUnknown = ch.fhEtaUnknown;

  return *this;
  
}
*/

//________________________________________________________________________
TList *  AliAnaChargedParticles::GetCreateOutputObjects()
{  
  // Create histograms to be saved in output file and 
  // store them in fOutputContainer
  
  
  TList * outputContainer = new TList() ; 
  outputContainer->SetName("ExampleHistos") ; 
  
  Int_t nptbins  = GetHistoPtBins();
  Int_t nphibins = GetHistoPhiBins();
  Int_t netabins = GetHistoEtaBins();
  Float_t ptmax  = GetHistoPtMax();
  Float_t phimax = GetHistoPhiMax();
  Float_t etamax = GetHistoEtaMax();
  Float_t ptmin  = GetHistoPtMin();
  Float_t phimin = GetHistoPhiMin();
  Float_t etamin = GetHistoEtaMin();	

  fhNtracks  = new TH1F ("hNtracks","# of tracks", 1000,0,1000); 
  fhNtracks->SetXTitle("# of tracks");
  outputContainer->Add(fhNtracks);
  
  fhVertex  = new TH3D ("Vertex","vertex position", 100,-50.,50., 100,-50.,50., 100,-50.,50.); 
  fhVertex->SetXTitle("X");
  fhVertex->SetYTitle("Y");
  fhVertex->SetZTitle("Z");
  outputContainer->Add(fhVertex);
  
  fhPt  = new TH1F ("hPtCharged","p_T distribution", nptbins,ptmin,ptmax); 
  fhPt->SetXTitle("p_{T} (GeV/c)");
  outputContainer->Add(fhPt);
  
  fhPhi  = new TH2F ("hPhiCharged","#phi distribution",nptbins,ptmin,ptmax, nphibins,phimin,phimax); 
  fhPhi->SetXTitle("#phi (rad)");
  outputContainer->Add(fhPhi);
  
  fhEta  = new TH2F ("hEtaCharged","#eta distribution",nptbins,ptmin,ptmax, netabins,etamin,etamax); 
  fhEta->SetXTitle("#eta ");
  outputContainer->Add(fhEta);
  
  fhPtEtaPhiPos  = new TH3D ("hPtEtaPhiPos","pt/eta/phi of positive charge",nptbins,ptmin,ptmax, netabins,etamin,etamax, nphibins,phimin,phimax); 
  fhPtEtaPhiPos->SetXTitle("p_{T}^{h^{+}} (GeV/c)");
  fhPtEtaPhiPos->SetYTitle("#eta ");
  fhPtEtaPhiPos->SetZTitle("#phi (rad)");  
  outputContainer->Add(fhPtEtaPhiPos);
  
  fhPtEtaPhiNeg  = new TH3D ("hPtEtaPhiNeg","pt/eta/phi of negative charge",nptbins,ptmin,ptmax, netabins,etamin,etamax, nphibins,phimin,phimax); 
  fhPtEtaPhiNeg->SetXTitle("p_{T}^{h^{-}} (GeV/c)");
  fhPtEtaPhiNeg->SetYTitle("#eta ");
  fhPtEtaPhiNeg->SetZTitle("#phi (rad)");  
  outputContainer->Add(fhPtEtaPhiNeg);
  
  if(IsDataMC()){
    
    fhPtPion  = new TH1F ("hPtMCPion","p_T distribution from #pi", nptbins,ptmin,ptmax); 
    fhPtPion->SetXTitle("p_{T} (GeV/c)");
    outputContainer->Add(fhPtPion);
    
    fhPhiPion  = new TH2F ("hPhiMCPion","#phi distribution from #pi",nptbins,ptmin,ptmax, nphibins,phimin,phimax); 
    fhPhiPion->SetXTitle("#phi (rad)");
    outputContainer->Add(fhPhiPion);
    
    fhEtaPion  = new TH2F ("hEtaMCPion","#eta distribution from #pi",nptbins,ptmin,ptmax, netabins,etamin,etamax); 
    fhEtaPion->SetXTitle("#eta ");
    outputContainer->Add(fhEtaPion);
    
    fhPtProton  = new TH1F ("hPtMCProton","p_T distribution from proton", nptbins,ptmin,ptmax); 
    fhPtProton->SetXTitle("p_{T} (GeV/c)");
    outputContainer->Add(fhPtProton);
    
    fhPhiProton  = new TH2F ("hPhiMCProton","#phi distribution from proton",nptbins,ptmin,ptmax, nphibins,phimin,phimax); 
    fhPhiProton->SetXTitle("#phi (rad)");
    outputContainer->Add(fhPhiProton);
    
    fhEtaProton  = new TH2F ("hEtaMCProton","#eta distribution from proton",nptbins,ptmin,ptmax, netabins,etamin,etamax); 
    fhEtaProton->SetXTitle("#eta ");
    outputContainer->Add(fhEtaProton);
    
    fhPtKaon  = new TH1F ("hPtMCKaon","p_T distribution from kaon", nptbins,ptmin,ptmax); 
    fhPtKaon->SetXTitle("p_{T} (GeV/c)");
    outputContainer->Add(fhPtKaon);
    
    fhPhiKaon  = new TH2F ("hPhiMCKaon","#phi distribution from kaon",nptbins,ptmin,ptmax, nphibins,phimin,phimax); 
    fhPhiKaon->SetXTitle("#phi (rad)");
    outputContainer->Add(fhPhiKaon);
    
    fhEtaKaon  = new TH2F ("hEtaMCKaon","#eta distribution from kaon",nptbins,ptmin,ptmax, netabins,etamin,etamax); 
    fhEtaKaon->SetXTitle("#eta ");
    outputContainer->Add(fhEtaKaon);
    
    fhPtElectron  = new TH1F ("hPtMCElectron","p_T distribution from electron", nptbins,ptmin,ptmax); 
    fhPtElectron->SetXTitle("p_{T} (GeV/c)");
    outputContainer->Add(fhPtElectron);
    
    fhPhiElectron  = new TH2F ("hPhiMCElectron","#phi distribution from electron",nptbins,ptmin,ptmax, nphibins,phimin,phimax); 
    fhPhiElectron->SetXTitle("#phi (rad)");
    outputContainer->Add(fhPhiElectron);
    
    fhEtaElectron  = new TH2F ("hEtaMCElectron","#eta distribution from electron",nptbins,ptmin,ptmax, netabins,etamin,etamax); 
    fhEtaElectron->SetXTitle("#eta ");
    outputContainer->Add(fhEtaElectron);
    
    fhPtUnknown  = new TH1F ("hPtMCUnknown","p_T distribution from unknown", nptbins,ptmin,ptmax); 
    fhPtUnknown->SetXTitle("p_{T} (GeV/c)");
    outputContainer->Add(fhPtUnknown);
    
    fhPhiUnknown  = new TH2F ("hPhiMCUnknown","#phi distribution from unknown",nptbins,ptmin,ptmax, nphibins,phimin,phimax); 
    fhPhiUnknown->SetXTitle("#phi (rad)");
    outputContainer->Add(fhPhiUnknown);
    
    fhEtaUnknown  = new TH2F ("hEtaMCUnknown","#eta distribution from unknown",nptbins,ptmin,ptmax, netabins,etamin,etamax); 
    fhEtaUnknown->SetXTitle("#eta ");
    outputContainer->Add(fhEtaUnknown);
    
  }
  
  return outputContainer;

}

//__________________________________________________
void AliAnaChargedParticles::InitParameters()
{ 
  //Initialize the parameters of the analysis.
  SetOutputAODClassName("AliAODPWG4Particle");
  SetOutputAODName("PWG4Particle");

  AddToHistogramsName("AnaCharged_");

  fPdg = -1; //Select all tracks 
  
}

//__________________________________________________________________
void AliAnaChargedParticles::Print(const Option_t * opt) const
{
  //Print some relevant parameters set for the analysis
  if(! opt)
    return;
  
  printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
  AliAnaPartCorrBaseClass::Print(" ");	
	
  printf("Min Pt = %3.2f\n", GetMinPt());
  printf("Max Pt = %3.2f\n", GetMaxPt());
  printf("Select clusters with pdg %d \n",fPdg);
  
} 

//____________________________________________________________________________
void AliAnaChargedParticles::Init()
{  
  //Init
  //Do some checks
  if(!GetReader()->IsCTSSwitchedOn()){
    printf("AliAnaChargedParticles::Init() - STOP!: You want to use CTS tracks in analysis but not read!! \n!!Check the configuration file!!\n");
    abort();
  }
  
}

//__________________________________________________________________
void  AliAnaChargedParticles::MakeAnalysisFillAOD() 
{
  //Do analysis and fill aods
  if(!GetAODCTS() || GetAODCTS()->GetEntriesFast() == 0) return ;
  Int_t ntracks = GetAODCTS()->GetEntriesFast();
  Double_t vert[3] = {0,0,0}; //vertex ;
  //Some prints
  if(GetDebug() > 0)
    printf("AliAnaChargedParticles::MakeAnalysisFillAOD() - In CTS aod entries %d\n", ntracks);
  
  //Fill AODParticle with CTS aods
  TVector3 p3;
  Int_t evtIndex = 0;
  for(Int_t i = 0; i < ntracks; i++){
    
    AliAODTrack * track =  (AliAODTrack*) (GetAODCTS()->At(i));
    
    //Fill AODParticle after some selection       
    Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
    p3.SetXYZ(mom[0],mom[1],mom[2]);
    
    //Acceptance selection
    Bool_t in =  GetFiducialCut()->IsInFiducialCut(mom,"CTS") ;
    if(GetDebug() > 1) printf("AliAnaChargedParticles::MakeAnalysisFillAOD() - Track pt %2.2f, phi %2.2f, in fiducial cut %d\n",p3.Pt(), p3.Phi(), in);
    if(p3.Pt() > GetMinPt() && in) {
      //Keep only particles identified with fPdg
      //Selection not done for the moment
      //Should be done here.
      if (GetMixedEvent()){
        evtIndex = GetMixedEvent()->EventIndex(track->GetID()) ;
      }        
      GetVertex(vert,evtIndex); 
      if(TMath::Abs(vert[2])> GetZvertexCut()) return;      
      AliAODPWG4Particle tr = AliAODPWG4Particle(mom[0],mom[1],mom[2],0);
      tr.SetDetector("CTS");
      tr.SetLabel(track->GetLabel());
      tr.SetTrackLabel(track->GetID(),-1);
    //  tr.SetChargedBit(track->Charge());
	  //Input from second AOD?
	  //if(GetReader()->GetAODCTSNormalInputEntries() <= i) tr.SetInputFileIndex(1);
		
      AddAODParticle(tr);
    }//selection
  }//loop
  
  if(GetDebug() > 0) 	
    printf("AliAnaChargedParticles::MakeAnalysisFillAOD() - Final aod branch entries %d\n", GetOutputAODBranch()->GetEntriesFast());   
} 

//__________________________________________________________________
void  AliAnaChargedParticles::MakeAnalysisFillHistograms() 
{
  //Do analysis and fill histograms
  
  //Loop on stored AODParticles
  Int_t naod = GetOutputAODBranch()->GetEntriesFast();
  Double_t v[3] = {0,0,0}; //vertex ;
  GetReader()->GetVertex(v);
  fhVertex->Fill(v[0],v[1],v[2]);
  if(TMath::Abs(v[2]) >GetZvertexCut()) return ;
  fhNtracks->Fill(GetReader()->GetTrackMultiplicity()) ;
  if(GetDebug() > 0) printf("AliAnaChargedParticles::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
  for(Int_t iaod = 0; iaod < naod ; iaod++){
    AliAODPWG4Particle* tr =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
        
    fhPt->Fill(tr->Pt());
    fhPhi->Fill(tr->Pt(), tr->Phi());
    fhEta->Fill(tr->Pt(), tr->Eta());
    //for charge information
    AliAODTrack * track =  (AliAODTrack*) (GetAODCTS()->At(iaod));
    if(track->Charge()>0)fhPtEtaPhiPos->Fill(tr->Pt(), tr->Eta(),tr->Phi());
    if(track->Charge()<0)fhPtEtaPhiNeg->Fill(tr->Pt(), tr->Eta(),tr->Phi());
    
    if(IsDataMC()){
      //Play with the MC stack if available		
      Int_t mompdg = -1;
      Int_t label  = TMath::Abs(tr->GetLabel());
      if(GetReader()->ReadStack()){
        TParticle * mom = GetMCStack()->Particle(label);
        mompdg =TMath::Abs(mom->GetPdgCode());
      }
      else if(GetReader()->ReadAODMCParticles()){
        AliAODMCParticle * aodmom = 0;
        //Get the list of MC particles
        aodmom = (AliAODMCParticle*) (GetReader()->GetAODMCParticles(tr->GetInputFileIndex()))->At(label);
        mompdg =TMath::Abs(aodmom->GetPdgCode());
      }
      
      
      if(mompdg==211){
        fhPtPion->Fill(tr->Pt());
        fhPhiPion->Fill(tr->Pt(), tr->Phi());
        fhEtaPion->Fill(tr->Pt(), tr->Eta());
      }
      else if(mompdg==2212){
        fhPtProton->Fill(tr->Pt());
        fhPhiProton->Fill(tr->Pt(), tr->Phi());
        fhEtaProton->Fill(tr->Pt(), tr->Eta());
      }
      else if(mompdg==321){
        fhPtKaon->Fill(tr->Pt());
        fhPhiKaon->Fill(tr->Pt(), tr->Phi());
        fhEtaKaon->Fill(tr->Pt(), tr->Eta());
      }
      else if(mompdg==11){
        fhPtElectron->Fill(tr->Pt());
        fhPhiElectron->Fill(tr->Pt(), tr->Phi());
        fhEtaElectron->Fill(tr->Pt(), tr->Eta());
      }
      else {
        //printf("unknown pdg %d\n",mompdg);
        fhPtUnknown->Fill(tr->Pt());
        fhPhiUnknown->Fill(tr->Pt(), tr->Phi());
        fhEtaUnknown->Fill(tr->Pt(), tr->Eta());
      }
    }//Work with stack also
    
  }// aod branch loop
  
}
Commit	Line	Data
477d6cee	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, 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	/* $Id: $ */
	16
	17	//_________________________________________________________________________
	18	//
	19	// Class for track selection and identification (not done now)
	20	// Tracks from the CTS are kept in the AOD.
	21	// Few histograms produced.
	22	//
	23	//-- Author: Gustavo Conesa (INFN-LNF)
	24	//_________________________________________________________________________
	25
	26
	27	// --- ROOT system ---
477d6cee	28	#include "TParticle.h"
477d6cee	29	#include "TH2F.h"
5025c139	30	#include "TH3D.h"
477d6cee	31	//---- AliRoot system ----
	32	#include "AliAnaChargedParticles.h"
	33	#include "AliCaloTrackReader.h"
	34	#include "AliAODPWG4Particle.h"
	35	#include "AliStack.h"
ff45398a	36	#include "AliFiducialCut.h"
477d6cee	37	#include "AliAODTrack.h"
591cc579	38	#include "AliAODMCParticle.h"
477d6cee	39
	40	ClassImp(AliAnaChargedParticles)
	41
	42	//____________________________________________________________________________
	43	AliAnaChargedParticles::AliAnaChargedParticles() :
5025c139	44	AliAnaPartCorrBaseClass(),fPdg(0), fhNtracks(0),fhVertex(0), fhPt(0),fhPhi(0),fhEta(0),
5025c139	45	fhPtEtaPhiPos(0), fhPtEtaPhiNeg(0),
477d6cee	46	fhPtPion(0),fhPhiPion(0),fhEtaPion(0),
	47	fhPtProton(0),fhPhiProton(0),fhEtaProton(0),
	48	fhPtElectron(0),fhPhiElectron(0),fhEtaElectron(0),
	49	fhPtKaon(0),fhPhiKaon(0),fhEtaKaon(0),
	50	fhPtUnknown(0),fhPhiUnknown(0),fhEtaUnknown(0)
	51
	52	{
	53	//Default Ctor
	54
	55	//Initialize parameters
	56	InitParameters();
	57	}
	58
78219bac	59	/*
477d6cee	60	//____________________________________________________________________________
	61	AliAnaChargedParticles::AliAnaChargedParticles(const AliAnaChargedParticles & ch) :
	62	AliAnaPartCorrBaseClass(ch), fPdg(ch.fPdg), fhPt(ch.fhPt), fhPhi(ch.fhPhi),fhEta(ch.fhEta),
	63	fhPtPion(ch. fhPtPion),fhPhiPion(ch.fhPhiPion),fhEtaPion(ch.fhEtaPion),
	64	fhPtProton(ch.fhPtProton),fhPhiProton(ch.fhPhiProton),fhEtaProton(ch.fhEtaProton),
	65	fhPtElectron(ch. fhPtElectron),fhPhiElectron(ch.fhPhiElectron),fhEtaElectron(ch.fhEtaElectron),
	66	fhPtKaon(ch. fhPtKaon),fhPhiKaon(ch.fhPhiKaon),fhEtaKaon(ch.fhEtaKaon),
	67	fhPtUnknown(ch.fhPtUnknown),fhPhiUnknown(ch.fhPhiUnknown),fhEtaUnknown(ch.fhEtaUnknown)
	68
	69	{
	70	// cpy ctor
	71
	72	}
	73
	74	//_________________________________________________________________________
	75	AliAnaChargedParticles & AliAnaChargedParticles::operator = (const AliAnaChargedParticles & ch)
	76	{
	77	// assignment operator
	78
	79	if(this == &ch)return *this;
	80	((AliAnaPartCorrBaseClass *)this)->operator=(ch);
	81
	82	fPdg = ch.fPdg;
	83	fhPt = ch.fhPt;
	84	fhPhi = ch.fhPhi;
	85	fhEta = ch.fhEta;
	86
	87	fhPtPion = ch. fhPtPion; fhPhiPion = ch.fhPhiPion; fhEtaPion = ch.fhEtaPion;
	88	fhPtProton = ch.fhPtProton; fhPhiProton = ch.fhPhiProton; fhEtaProton = ch.fhEtaProton;
	89	fhPtElectron = ch. fhPtElectron; fhPhiElectron = ch.fhPhiElectron; fhEtaElectron = ch.fhEtaElectron;
	90	fhPtKaon = ch. fhPtKaon; fhPhiKaon = ch.fhPhiKaon; fhEtaKaon = ch.fhEtaKaon;
	91	fhPtUnknown = ch.fhPtUnknown; fhPhiUnknown = ch.fhPhiUnknown; fhEtaUnknown = ch.fhEtaUnknown;
	92
	93	return *this;
	94
	95	}
78219bac	96	*/
477d6cee	97
	98	//________________________________________________________________________
	99	TList * AliAnaChargedParticles::GetCreateOutputObjects()
	100	{
	101	// Create histograms to be saved in output file and
	102	// store them in fOutputContainer
	103
	104
	105	TList * outputContainer = new TList() ;
	106	outputContainer->SetName("ExampleHistos") ;
	107
5a2dbc3c	108	Int_t nptbins = GetHistoPtBins();
	109	Int_t nphibins = GetHistoPhiBins();
	110	Int_t netabins = GetHistoEtaBins();
477d6cee	111	Float_t ptmax = GetHistoPtMax();
	112	Float_t phimax = GetHistoPhiMax();
	113	Float_t etamax = GetHistoEtaMax();
	114	Float_t ptmin = GetHistoPtMin();
	115	Float_t phimin = GetHistoPhiMin();
	116	Float_t etamin = GetHistoEtaMin();
5025c139	117
	118	fhNtracks = new TH1F ("hNtracks","# of tracks", 1000,0,1000);
	119	fhNtracks->SetXTitle("# of tracks");
	120	outputContainer->Add(fhNtracks);
	121
	122	fhVertex = new TH3D ("Vertex","vertex position", 100,-50.,50., 100,-50.,50., 100,-50.,50.);
	123	fhVertex->SetXTitle("X");
	124	fhVertex->SetYTitle("Y");
	125	fhVertex->SetZTitle("Z");
	126	outputContainer->Add(fhVertex);
477d6cee	127
	128	fhPt = new TH1F ("hPtCharged","p_T distribution", nptbins,ptmin,ptmax);
	129	fhPt->SetXTitle("p_{T} (GeV/c)");
	130	outputContainer->Add(fhPt);
	131
	132	fhPhi = new TH2F ("hPhiCharged","#phi distribution",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
	133	fhPhi->SetXTitle("#phi (rad)");
	134	outputContainer->Add(fhPhi);
	135
	136	fhEta = new TH2F ("hEtaCharged","#eta distribution",nptbins,ptmin,ptmax, netabins,etamin,etamax);
	137	fhEta->SetXTitle("#eta ");
	138	outputContainer->Add(fhEta);
	139
5025c139	140	fhPtEtaPhiPos = new TH3D ("hPtEtaPhiPos","pt/eta/phi of positive charge",nptbins,ptmin,ptmax, netabins,etamin,etamax, nphibins,phimin,phimax);
	141	fhPtEtaPhiPos->SetXTitle("p_{T}^{h^{+}} (GeV/c)");
	142	fhPtEtaPhiPos->SetYTitle("#eta ");
	143	fhPtEtaPhiPos->SetZTitle("#phi (rad)");
	144	outputContainer->Add(fhPtEtaPhiPos);
	145
	146	fhPtEtaPhiNeg = new TH3D ("hPtEtaPhiNeg","pt/eta/phi of negative charge",nptbins,ptmin,ptmax, netabins,etamin,etamax, nphibins,phimin,phimax);
	147	fhPtEtaPhiNeg->SetXTitle("p_{T}^{h^{-}} (GeV/c)");
	148	fhPtEtaPhiNeg->SetYTitle("#eta ");
	149	fhPtEtaPhiNeg->SetZTitle("#phi (rad)");
	150	outputContainer->Add(fhPtEtaPhiNeg);
477d6cee	151
	152	if(IsDataMC()){
	153
a3aebfff	154	fhPtPion = new TH1F ("hPtMCPion","p_T distribution from #pi", nptbins,ptmin,ptmax);
477d6cee	155	fhPtPion->SetXTitle("p_{T} (GeV/c)");
	156	outputContainer->Add(fhPtPion);
	157
a3aebfff	158	fhPhiPion = new TH2F ("hPhiMCPion","#phi distribution from #pi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
477d6cee	159	fhPhiPion->SetXTitle("#phi (rad)");
	160	outputContainer->Add(fhPhiPion);
	161
a3aebfff	162	fhEtaPion = new TH2F ("hEtaMCPion","#eta distribution from #pi",nptbins,ptmin,ptmax, netabins,etamin,etamax);
477d6cee	163	fhEtaPion->SetXTitle("#eta ");
	164	outputContainer->Add(fhEtaPion);
	165
a3aebfff	166	fhPtProton = new TH1F ("hPtMCProton","p_T distribution from proton", nptbins,ptmin,ptmax);
477d6cee	167	fhPtProton->SetXTitle("p_{T} (GeV/c)");
	168	outputContainer->Add(fhPtProton);
	169
a3aebfff	170	fhPhiProton = new TH2F ("hPhiMCProton","#phi distribution from proton",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
477d6cee	171	fhPhiProton->SetXTitle("#phi (rad)");
	172	outputContainer->Add(fhPhiProton);
	173
a3aebfff	174	fhEtaProton = new TH2F ("hEtaMCProton","#eta distribution from proton",nptbins,ptmin,ptmax, netabins,etamin,etamax);
477d6cee	175	fhEtaProton->SetXTitle("#eta ");
	176	outputContainer->Add(fhEtaProton);
	177
a3aebfff	178	fhPtKaon = new TH1F ("hPtMCKaon","p_T distribution from kaon", nptbins,ptmin,ptmax);
477d6cee	179	fhPtKaon->SetXTitle("p_{T} (GeV/c)");
	180	outputContainer->Add(fhPtKaon);
	181
a3aebfff	182	fhPhiKaon = new TH2F ("hPhiMCKaon","#phi distribution from kaon",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
477d6cee	183	fhPhiKaon->SetXTitle("#phi (rad)");
	184	outputContainer->Add(fhPhiKaon);
	185
a3aebfff	186	fhEtaKaon = new TH2F ("hEtaMCKaon","#eta distribution from kaon",nptbins,ptmin,ptmax, netabins,etamin,etamax);
477d6cee	187	fhEtaKaon->SetXTitle("#eta ");
	188	outputContainer->Add(fhEtaKaon);
	189
a3aebfff	190	fhPtElectron = new TH1F ("hPtMCElectron","p_T distribution from electron", nptbins,ptmin,ptmax);
477d6cee	191	fhPtElectron->SetXTitle("p_{T} (GeV/c)");
	192	outputContainer->Add(fhPtElectron);
	193
a3aebfff	194	fhPhiElectron = new TH2F ("hPhiMCElectron","#phi distribution from electron",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
477d6cee	195	fhPhiElectron->SetXTitle("#phi (rad)");
	196	outputContainer->Add(fhPhiElectron);
	197
a3aebfff	198	fhEtaElectron = new TH2F ("hEtaMCElectron","#eta distribution from electron",nptbins,ptmin,ptmax, netabins,etamin,etamax);
477d6cee	199	fhEtaElectron->SetXTitle("#eta ");
	200	outputContainer->Add(fhEtaElectron);
	201
a3aebfff	202	fhPtUnknown = new TH1F ("hPtMCUnknown","p_T distribution from unknown", nptbins,ptmin,ptmax);
477d6cee	203	fhPtUnknown->SetXTitle("p_{T} (GeV/c)");
	204	outputContainer->Add(fhPtUnknown);
	205
a3aebfff	206	fhPhiUnknown = new TH2F ("hPhiMCUnknown","#phi distribution from unknown",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
477d6cee	207	fhPhiUnknown->SetXTitle("#phi (rad)");
	208	outputContainer->Add(fhPhiUnknown);
	209
a3aebfff	210	fhEtaUnknown = new TH2F ("hEtaMCUnknown","#eta distribution from unknown",nptbins,ptmin,ptmax, netabins,etamin,etamax);
477d6cee	211	fhEtaUnknown->SetXTitle("#eta ");
	212	outputContainer->Add(fhEtaUnknown);
	213
	214	}
a3aebfff	215
477d6cee	216	return outputContainer;
a3aebfff	217
477d6cee	218	}
	219
	220	//__________________________________________________
	221	void AliAnaChargedParticles::InitParameters()
	222	{
	223	//Initialize the parameters of the analysis.
	224	SetOutputAODClassName("AliAODPWG4Particle");
a3aebfff	225	SetOutputAODName("PWG4Particle");
	226
	227	AddToHistogramsName("AnaCharged_");
	228
477d6cee	229	fPdg = -1; //Select all tracks
	230
	231	}
	232
	233	//__________________________________________________________________
	234	void AliAnaChargedParticles::Print(const Option_t * opt) const
	235	{
	236	//Print some relevant parameters set for the analysis
	237	if(! opt)
	238	return;
	239
a3aebfff	240	printf("** Print %s %s **\n", GetName(), GetTitle() ) ;
	241	AliAnaPartCorrBaseClass::Print(" ");
	242
477d6cee	243	printf("Min Pt = %3.2f\n", GetMinPt());
	244	printf("Max Pt = %3.2f\n", GetMaxPt());
	245	printf("Select clusters with pdg %d \n",fPdg);
	246
	247	}
	248
	249	//____________________________________________________________________________
	250	void AliAnaChargedParticles::Init()
	251	{
	252	//Init
	253	//Do some checks
	254	if(!GetReader()->IsCTSSwitchedOn()){
7cd4e982	255	printf("AliAnaChargedParticles::Init() - STOP!: You want to use CTS tracks in analysis but not read!! \n!!Check the configuration file!!\n");
477d6cee	256	abort();
	257	}
	258
	259	}
	260
	261	//__________________________________________________________________
	262	void AliAnaChargedParticles::MakeAnalysisFillAOD()
	263	{
	264	//Do analysis and fill aods
	265	if(!GetAODCTS() \|\| GetAODCTS()->GetEntriesFast() == 0) return ;
	266	Int_t ntracks = GetAODCTS()->GetEntriesFast();
2244659d	267	Double_t vert[3] = {0,0,0}; //vertex ;
477d6cee	268	//Some prints
477d6cee	269	if(GetDebug() > 0)
a3aebfff	270	printf("AliAnaChargedParticles::MakeAnalysisFillAOD() - In CTS aod entries %d\n", ntracks);
477d6cee	271
	272	//Fill AODParticle with CTS aods
	273	TVector3 p3;
5025c139	274	Int_t evtIndex = 0;
477d6cee	275	for(Int_t i = 0; i < ntracks; i++){
	276
	277	AliAODTrack * track = (AliAODTrack*) (GetAODCTS()->At(i));
	278
	279	//Fill AODParticle after some selection
	280	Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
	281	p3.SetXYZ(mom[0],mom[1],mom[2]);
	282
	283	//Acceptance selection
ff45398a	284	Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,"CTS") ;
ff45398a	285	if(GetDebug() > 1) printf("AliAnaChargedParticles::MakeAnalysisFillAOD() - Track pt %2.2f, phi %2.2f, in fiducial cut %d\n",p3.Pt(), p3.Phi(), in);
477d6cee	286	if(p3.Pt() > GetMinPt() && in) {
	287	//Keep only particles identified with fPdg
	288	//Selection not done for the moment
	289	//Should be done here.
5025c139	290	if (GetMixedEvent()){
	291	evtIndex = GetMixedEvent()->EventIndex(track->GetID()) ;
	292	}
	293	GetVertex(vert,evtIndex);
2244659d	294	if(TMath::Abs(vert[2])> GetZvertexCut()) return;
477d6cee	295	AliAODPWG4Particle tr = AliAODPWG4Particle(mom[0],mom[1],mom[2],0);
	296	tr.SetDetector("CTS");
	297	tr.SetLabel(track->GetLabel());
	298	tr.SetTrackLabel(track->GetID(),-1);
5025c139	299	// tr.SetChargedBit(track->Charge());
7cd4e982	300	//Input from second AOD?
1e68a3f4	301	//if(GetReader()->GetAODCTSNormalInputEntries() <= i) tr.SetInputFileIndex(1);
7cd4e982	302
477d6cee	303	AddAODParticle(tr);
	304	}//selection
	305	}//loop
	306
	307	if(GetDebug() > 0)
a3aebfff	308	printf("AliAnaChargedParticles::MakeAnalysisFillAOD() - Final aod branch entries %d\n", GetOutputAODBranch()->GetEntriesFast());
477d6cee	309	}
	310
	311	//__________________________________________________________________
	312	void AliAnaChargedParticles::MakeAnalysisFillHistograms()
	313	{
	314	//Do analysis and fill histograms
	315
	316	//Loop on stored AODParticles
	317	Int_t naod = GetOutputAODBranch()->GetEntriesFast();
5025c139	318	Double_t v[3] = {0,0,0}; //vertex ;
	319	GetReader()->GetVertex(v);
	320	fhVertex->Fill(v[0],v[1],v[2]);
2244659d	321	if(TMath::Abs(v[2]) >GetZvertexCut()) return ;
2244659d	322	fhNtracks->Fill(GetReader()->GetTrackMultiplicity()) ;
a3aebfff	323	if(GetDebug() > 0) printf("AliAnaChargedParticles::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
477d6cee	324	for(Int_t iaod = 0; iaod < naod ; iaod++){
477d6cee	325	AliAODPWG4Particle* tr = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
2244659d	326
477d6cee	327	fhPt->Fill(tr->Pt());
	328	fhPhi->Fill(tr->Pt(), tr->Phi());
	329	fhEta->Fill(tr->Pt(), tr->Eta());
5025c139	330	//for charge information
	331	AliAODTrack * track = (AliAODTrack*) (GetAODCTS()->At(iaod));
	332	if(track->Charge()>0)fhPtEtaPhiPos->Fill(tr->Pt(), tr->Eta(),tr->Phi());
	333	if(track->Charge()<0)fhPtEtaPhiNeg->Fill(tr->Pt(), tr->Eta(),tr->Phi());
477d6cee	334
	335	if(IsDataMC()){
	336	//Play with the MC stack if available
5025c139	337	Int_t mompdg = -1;
f6661bbc	338	Int_t label = TMath::Abs(tr->GetLabel());
5025c139	339	if(GetReader()->ReadStack()){
	340	TParticle * mom = GetMCStack()->Particle(label);
	341	mompdg =TMath::Abs(mom->GetPdgCode());
591cc579	342	}
5025c139	343	else if(GetReader()->ReadAODMCParticles()){
	344	AliAODMCParticle * aodmom = 0;
	345	//Get the list of MC particles
	346	aodmom = (AliAODMCParticle*) (GetReader()->GetAODMCParticles(tr->GetInputFileIndex()))->At(label);
	347	mompdg =TMath::Abs(aodmom->GetPdgCode());
	348	}
	349
f6661bbc	350
5025c139	351	if(mompdg==211){
	352	fhPtPion->Fill(tr->Pt());
	353	fhPhiPion->Fill(tr->Pt(), tr->Phi());
	354	fhEtaPion->Fill(tr->Pt(), tr->Eta());
477d6cee	355	}
477d6cee	356	else if(mompdg==2212){
5025c139	357	fhPtProton->Fill(tr->Pt());
	358	fhPhiProton->Fill(tr->Pt(), tr->Phi());
	359	fhEtaProton->Fill(tr->Pt(), tr->Eta());
477d6cee	360	}
477d6cee	361	else if(mompdg==321){
5025c139	362	fhPtKaon->Fill(tr->Pt());
	363	fhPhiKaon->Fill(tr->Pt(), tr->Phi());
	364	fhEtaKaon->Fill(tr->Pt(), tr->Eta());
477d6cee	365	}
477d6cee	366	else if(mompdg==11){
5025c139	367	fhPtElectron->Fill(tr->Pt());
	368	fhPhiElectron->Fill(tr->Pt(), tr->Phi());
	369	fhEtaElectron->Fill(tr->Pt(), tr->Eta());
477d6cee	370	}
477d6cee	371	else {
5025c139	372	//printf("unknown pdg %d\n",mompdg);
	373	fhPtUnknown->Fill(tr->Pt());
	374	fhPhiUnknown->Fill(tr->Pt(), tr->Phi());
	375	fhEtaUnknown->Fill(tr->Pt(), tr->Eta());
477d6cee	376	}
477d6cee	377	}//Work with stack also
f6661bbc	378
477d6cee	379	}// aod branch loop
	380
	381	}