[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / AliAnalysisTaskSpectraAllChAOD.cxx

/**************************************************************************
 * Copyright(c) 1998-2009, 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.                  *
 **************************************************************************/

//-----------------------------------------------------------------
//         AliAnalysisTaskSpectraAllChAOD class
//-----------------------------------------------------------------

#include "TChain.h"
#include "TTree.h"
#include "TLegend.h"
#include "TH1F.h"
#include "TH2F.h"
#include "THnSparse.h"
#include "TCanvas.h"
#include "AliAnalysisTask.h"
#include "AliAODTrack.h"
#include "AliAODMCParticle.h"
#include "AliVParticle.h"
#include "AliAODEvent.h"
#include "AliAODInputHandler.h"
#include "AliAnalysisTaskSpectraAllChAOD.h"
#include "AliAnalysisTaskESDfilter.h"
#include "AliAnalysisDataContainer.h"
#include "AliSpectraAODTrackCuts.h"
#include "AliSpectraAODEventCuts.h"
#include "AliHelperPID.h"
#include "AliPIDCombined.h"
#include "AliCentrality.h"
#include "TProof.h"
#include "AliVEvent.h"
#include "AliStack.h"
#include <TMCProcess.h>

#include <iostream>

using namespace AliHelperPIDNameSpace;
using namespace std;

ClassImp(AliAnalysisTaskSpectraAllChAOD)

//________________________________________________________________________
AliAnalysisTaskSpectraAllChAOD::AliAnalysisTaskSpectraAllChAOD(const char *name) : AliAnalysisTaskSE(name),  
  fAOD(0x0),
  fTrackCuts(0x0),
  fEventCuts(0x0),
  fHelperPID(0x0),
  fIsMC(0),
  fDoDoubleCounting(0),
  fFillOnlyEvents(0),
  fCharge(0),
  fVZEROside(0),
  fOutput(0x0),
  fnCentBins(20),
  fnQvecBins(100),
  fnNchBins(200),
  fIsQvecCalibMode(0),
  fQvecUpperLim(100)
{
  // Default constructor
  DefineInput(0, TChain::Class());
  DefineOutput(1, TList::Class());
  DefineOutput(2, AliSpectraAODEventCuts::Class());
  DefineOutput(3, AliSpectraAODTrackCuts::Class());
  DefineOutput(4, AliHelperPID::Class());
}

//________________________________________________________________________
void AliAnalysisTaskSpectraAllChAOD::UserCreateOutputObjects()
{
  // create output objects
  fOutput=new TList();
  fOutput->SetOwner();
  fOutput->SetName("fOutput");
  
  if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");
  if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");
  if (!fHelperPID)  AliFatal("HelperPID should be set in the steering macro");
  
  // binning common to all the THn
  const Double_t ptBins[] = {0.20,0.30,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.4,2.8,3.2,3.6,4.0,5.0,6.0,7.0,8.0,9.0,10.,12.,15.};
  const Int_t nptBins=26;
  
  //dimensions of THnSparse for tracks
  const Int_t nvartrk=7;
  //                                             pt          cent          Q vec     IDrec      IDgen       isph         y
  Int_t    binsHistRealTrk[nvartrk] = {      nptBins, fnCentBins,     fnQvecBins,        4,        3,         2,        2};
  Double_t xminHistRealTrk[nvartrk] = {         0.,          0.,              0.,       -.5,      -0.5,      -0.5,    -0.5};
  Double_t xmaxHistRealTrk[nvartrk] = {       10.,       100., fQvecUpperLim,      3.5,      2.5,       1.5,     0.5};    
  THnSparseF* NSparseHistTrk = new THnSparseF("NSparseHistTrk","NSparseHistTrk",nvartrk,binsHistRealTrk,xminHistRealTrk,xmaxHistRealTrk);
  NSparseHistTrk->GetAxis(0)->SetTitle("#it{p}_{T,rec}");
  NSparseHistTrk->GetAxis(0)->SetName("pT_rec");
  NSparseHistTrk->SetBinEdges(0,ptBins);
  NSparseHistTrk->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
  NSparseHistTrk->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
  NSparseHistTrk->GetAxis(2)->SetTitle("Q vec");
  NSparseHistTrk->GetAxis(2)->SetName("Q_vec");
  NSparseHistTrk->GetAxis(3)->SetTitle("ID rec");
  NSparseHistTrk->GetAxis(3)->SetName("ID_rec");
  NSparseHistTrk->GetAxis(4)->SetTitle("ID gen");
  NSparseHistTrk->GetAxis(4)->SetName("ID_gen");
  NSparseHistTrk->GetAxis(5)->SetTitle("isph");
  NSparseHistTrk->GetAxis(5)->SetName("isph");
  NSparseHistTrk->GetAxis(6)->SetTitle("y");
  NSparseHistTrk->GetAxis(6)->SetName("y");
  fOutput->Add(NSparseHistTrk);
  
  //dimensions of THnSparse for stack
  const Int_t nvarst=6;
  //                                             pt          cent    IDgen        isph        y    etaselected
  Int_t    binsHistRealSt[nvarst] = {      nptBins,  fnCentBins,        3,         2,        2,             1};
  Double_t xminHistRealSt[nvarst] = {         0.,           0.,      -0.5,      -0.5,    -0.5,            0.5};
  Double_t xmaxHistRealSt[nvarst] = {       10.,        100.,      2.5,       1.5,      0.5,           1.5};
  THnSparseF* NSparseHistSt = new THnSparseF("NSparseHistSt","NSparseHistSt",nvarst,binsHistRealSt,xminHistRealSt,xmaxHistRealSt);
  NSparseHistSt->GetAxis(0)->SetTitle("#it{p}_{T,gen}");
  NSparseHistSt->SetBinEdges(0,ptBins);
  NSparseHistSt->GetAxis(0)->SetName("pT_rec");
  NSparseHistSt->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
  NSparseHistSt->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
  NSparseHistSt->GetAxis(2)->SetTitle("ID gen");
  NSparseHistSt->GetAxis(2)->SetName("ID_gen");
  NSparseHistSt->GetAxis(3)->SetTitle("isph");
  NSparseHistSt->GetAxis(3)->SetName("isph");
  NSparseHistSt->GetAxis(4)->SetTitle("y");
  NSparseHistSt->GetAxis(4)->SetName("y");
  NSparseHistSt->GetAxis(5)->SetTitle("etaselected");
  NSparseHistSt->GetAxis(5)->SetName("etaselected");
  fOutput->Add(NSparseHistSt);
  
  //dimensions of THnSparse for the normalization
  const Int_t nvarev=3;
  //                                             cent             Q vec                Nch
  Int_t    binsHistRealEv[nvarev] = {    fnCentBins,      fnQvecBins,           fnNchBins};
  Double_t xminHistRealEv[nvarev] = {           0.,               0.,                   0.};
  Double_t xmaxHistRealEv[nvarev] = {       100.,  fQvecUpperLim,               2000.};
  THnSparseF* NSparseHistEv = new THnSparseF("NSparseHistEv","NSparseHistEv",nvarev,binsHistRealEv,xminHistRealEv,xmaxHistRealEv);
  NSparseHistEv->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
  NSparseHistEv->GetAxis(0)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
  NSparseHistEv->GetAxis(1)->SetTitle("Q vec");
  NSparseHistEv->GetAxis(1)->SetName("Q_vec");
  NSparseHistEv->GetAxis(2)->SetTitle("N charged");
  NSparseHistEv->GetAxis(2)->SetName("N_ch");
  fOutput->Add(NSparseHistEv);
  
  PostData(1, fOutput  );
  PostData(2, fEventCuts);
  PostData(3, fTrackCuts);
  PostData(4, fHelperPID);
}

//________________________________________________________________________

void AliAnalysisTaskSpectraAllChAOD::UserExec(Option_t *)
{
  //Printf("An event");
  // main event loop
  fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);
  if (!fAOD) {
    AliWarning("ERROR: AliAODEvent not available \n");
    return;
  }
  
  if (strcmp(fAOD->ClassName(), "AliAODEvent"))
    {
      AliFatal("Not processing AODs");
    }
  
  if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection

  //Default TPC priors
  if(fHelperPID->GetPIDType()==kBayes)fHelperPID->GetPIDCombined()->SetDefaultTPCPriors();//FIXME maybe this can go in the UserCreateOutputObject?
  
  Double_t Qvec=0.;//in case of MC we save space in the memory
  if(!fIsMC){
    if(fIsQvecCalibMode){
      if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();
      else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();
    }
    else Qvec=fEventCuts->GetQvecPercentile(fVZEROside);
  }
  
  Double_t Cent=fEventCuts->GetCent();
  
  // First do MC to fill up the MC particle array
  TClonesArray *arrayMC = 0;
  if (fIsMC)
    {
      arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
      if (!arrayMC) {
	AliFatal("Error: MC particles branch not found!\n");
      }
      Int_t nMC = arrayMC->GetEntries();
      for (Int_t iMC = 0; iMC < nMC; iMC++)
	{
	  AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
	  if(!partMC->Charge()) continue;//Skip neutrals
	  if(fCharge != 0 && partMC->Charge()*fCharge < 0.) continue;//if fCharge != 0 only select fCharge
	  
	  //flag to select particle in the same ETA acceptance of the tracks (to be used for the comparison with AllCh)
	  Double_t etaselected=-1.;
	  if(partMC->Eta()<fTrackCuts->GetEtaMin() || partMC->Eta()>fTrackCuts->GetEtaMax())etaselected=1.;
	  
	  //pt     cent    IDgen        isph        y
 	  Double_t varSt[6];
	  varSt[0]=partMC->Pt();
	  varSt[1]=Cent;
	  varSt[2]=(Double_t)fHelperPID->GetParticleSpecies(partMC);
	  varSt[3]=(Double_t)partMC->IsPhysicalPrimary();
	  varSt[4]=partMC->Y();
	  varSt[5]=etaselected;
	  ((THnSparseF*)fOutput->FindObject("NSparseHistSt"))->Fill(varSt);//stack loop
	  
	  //Printf("a particle");
	  
	}
    }
  
  //main loop on tracks
  
  Int_t Nch = 0.;
  
  for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
    AliAODTrack* track = fAOD->GetTrack(iTracks);
    if(fCharge != 0 && track->Charge() != fCharge) continue;//if fCharge != 0 only select fCharge 
    if (!fTrackCuts->IsSelected(track,kTRUE)) continue; //track selection (rapidity selection NOT in the standard cuts)
    if(!fFillOnlyEvents){
      Int_t IDrec=fHelperPID->GetParticleSpecies(track,kTRUE);//id from detector
      Double_t y= track->Y(fHelperPID->GetMass((AliHelperParticleSpecies_t)IDrec));
      Int_t IDgen=kSpUndefined;//set if MC
      Int_t isph=-999;
      //Int_t iswd=-999;
      
      if (arrayMC) {
	AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));
	if (!partMC) { 
	  AliError("Cannot get MC particle");
	  continue; 
	}
	IDgen=fHelperPID->GetParticleSpecies(partMC);
	isph=partMC->IsPhysicalPrimary();
	//iswd=partMC->IsSecondaryFromWeakDecay();//FIXME not working on old productions - removed Apr 8th 2014
      }
      
    //pt     cent    Q vec     IDrec     IDgen       isph      y
      Double_t varTrk[7];
      varTrk[0]=track->Pt();
      varTrk[1]=Cent;
      varTrk[2]=Qvec;
      varTrk[3]=(Double_t)IDrec;
      varTrk[4]=(Double_t)IDgen;
      varTrk[5]=(Double_t)isph;
      varTrk[6]=y;
      ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
      
      //for nsigma PID fill double counting of ID
      if(fHelperPID->GetPIDType()<kBayes && fDoDoubleCounting){//only nsigma
	Bool_t *HasDC;
	HasDC=fHelperPID->GetDoubleCounting(track,kTRUE);//get the array with double counting
	for(Int_t ipart=0;ipart<kNSpecies;ipart++){
	  if(HasDC[ipart]==kTRUE){
	    varTrk[3]=(Double_t)ipart;
	    ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
	  }
	}
      }
      
      //fill all charged (3)
      varTrk[3]=3.;
      varTrk[4]=3.;
      ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
    }//end if fFillOnlyEvents
    
    //Printf("a track");
      
    Nch++;
  } // end loop on tracks
  
  Double_t varEv[3];
  varEv[0]=Cent;
  varEv[1]=Qvec;
  varEv[2]=Nch;
  ((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
  
  PostData(1, fOutput  );
  PostData(2, fEventCuts);
  PostData(3, fTrackCuts);
  PostData(4, fHelperPID);
}

//_________________________________________________________________
void   AliAnalysisTaskSpectraAllChAOD::Terminate(Option_t *)
{
  // Terminate
}
Commit	Line	Data
c683985a	1	/**************************************************************************
	2	* Copyright(c) 1998-2009, 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	// AliAnalysisTaskSpectraAllChAOD class
	18	//-----------------------------------------------------------------
	19
	20	#include "TChain.h"
	21	#include "TTree.h"
	22	#include "TLegend.h"
	23	#include "TH1F.h"
	24	#include "TH2F.h"
	25	#include "THnSparse.h"
	26	#include "TCanvas.h"
	27	#include "AliAnalysisTask.h"
	28	#include "AliAODTrack.h"
	29	#include "AliAODMCParticle.h"
	30	#include "AliVParticle.h"
	31	#include "AliAODEvent.h"
	32	#include "AliAODInputHandler.h"
	33	#include "AliAnalysisTaskSpectraAllChAOD.h"
	34	#include "AliAnalysisTaskESDfilter.h"
	35	#include "AliAnalysisDataContainer.h"
	36	#include "AliSpectraAODTrackCuts.h"
	37	#include "AliSpectraAODEventCuts.h"
	38	#include "AliHelperPID.h"
	39	#include "AliPIDCombined.h"
	40	#include "AliCentrality.h"
	41	#include "TProof.h"
	42	#include "AliVEvent.h"
	43	#include "AliStack.h"
	44	#include <TMCProcess.h>
	45
	46	#include <iostream>
	47
	48	using namespace AliHelperPIDNameSpace;
	49	using namespace std;
	50
	51	ClassImp(AliAnalysisTaskSpectraAllChAOD)
	52
	53	//________________________________________________________________________
	54	AliAnalysisTaskSpectraAllChAOD::AliAnalysisTaskSpectraAllChAOD(const char *name) : AliAnalysisTaskSE(name),
	55	fAOD(0x0),
	56	fTrackCuts(0x0),
	57	fEventCuts(0x0),
	58	fHelperPID(0x0),
	59	fIsMC(0),
bc07eb4b	60	fDoDoubleCounting(0),
dbbc07ac	61	fFillOnlyEvents(0),
c683985a	62	fCharge(0),
	63	fVZEROside(0),
	64	fOutput(0x0),
	65	fnCentBins(20),
eba59d12	66	fnQvecBins(100),
	67	fnNchBins(200),
	68	fIsQvecCalibMode(0),
	69	fQvecUpperLim(100)
c683985a	70	{
	71	// Default constructor
	72	DefineInput(0, TChain::Class());
	73	DefineOutput(1, TList::Class());
	74	DefineOutput(2, AliSpectraAODEventCuts::Class());
	75	DefineOutput(3, AliSpectraAODTrackCuts::Class());
	76	DefineOutput(4, AliHelperPID::Class());
	77	}
	78
	79	//________________________________________________________________________
	80	void AliAnalysisTaskSpectraAllChAOD::UserCreateOutputObjects()
	81	{
	82	// create output objects
	83	fOutput=new TList();
	84	fOutput->SetOwner();
	85	fOutput->SetName("fOutput");
	86
	87	if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");
	88	if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");
	89	if (!fHelperPID) AliFatal("HelperPID should be set in the steering macro");
	90
	91	// binning common to all the THn
	92	const Double_t ptBins[] = {0.20,0.30,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.4,2.8,3.2,3.6,4.0,5.0,6.0,7.0,8.0,9.0,10.,12.,15.};
	93	const Int_t nptBins=26;
	94
	95	//dimensions of THnSparse for tracks
eba59d12	96	const Int_t nvartrk=7;
	97	// pt cent Q vec IDrec IDgen isph y
	98	Int_t binsHistRealTrk[nvartrk] = { nptBins, fnCentBins, fnQvecBins, 4, 3, 2, 2};
	99	Double_t xminHistRealTrk[nvartrk] = { 0., 0., 0., -.5, -0.5, -0.5, -0.5};
	100	Double_t xmaxHistRealTrk[nvartrk] = { 10., 100., fQvecUpperLim, 3.5, 2.5, 1.5, 0.5};
c683985a	101	THnSparseF* NSparseHistTrk = new THnSparseF("NSparseHistTrk","NSparseHistTrk",nvartrk,binsHistRealTrk,xminHistRealTrk,xmaxHistRealTrk);
	102	NSparseHistTrk->GetAxis(0)->SetTitle("#it{p}_{T,rec}");
	103	NSparseHistTrk->GetAxis(0)->SetName("pT_rec");
	104	NSparseHistTrk->SetBinEdges(0,ptBins);
	105	NSparseHistTrk->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
	106	NSparseHistTrk->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
	107	NSparseHistTrk->GetAxis(2)->SetTitle("Q vec");
	108	NSparseHistTrk->GetAxis(2)->SetName("Q_vec");
	109	NSparseHistTrk->GetAxis(3)->SetTitle("ID rec");
	110	NSparseHistTrk->GetAxis(3)->SetName("ID_rec");
	111	NSparseHistTrk->GetAxis(4)->SetTitle("ID gen");
	112	NSparseHistTrk->GetAxis(4)->SetName("ID_gen");
	113	NSparseHistTrk->GetAxis(5)->SetTitle("isph");
	114	NSparseHistTrk->GetAxis(5)->SetName("isph");
eba59d12	115	NSparseHistTrk->GetAxis(6)->SetTitle("y");
eba59d12	116	NSparseHistTrk->GetAxis(6)->SetName("y");
c683985a	117	fOutput->Add(NSparseHistTrk);
	118
	119	//dimensions of THnSparse for stack
eba59d12	120	const Int_t nvarst=6;
	121	// pt cent IDgen isph y etaselected
	122	Int_t binsHistRealSt[nvarst] = { nptBins, fnCentBins, 3, 2, 2, 1};
	123	Double_t xminHistRealSt[nvarst] = { 0., 0., -0.5, -0.5, -0.5, 0.5};
	124	Double_t xmaxHistRealSt[nvarst] = { 10., 100., 2.5, 1.5, 0.5, 1.5};
c683985a	125	THnSparseF* NSparseHistSt = new THnSparseF("NSparseHistSt","NSparseHistSt",nvarst,binsHistRealSt,xminHistRealSt,xmaxHistRealSt);
	126	NSparseHistSt->GetAxis(0)->SetTitle("#it{p}_{T,gen}");
	127	NSparseHistSt->SetBinEdges(0,ptBins);
	128	NSparseHistSt->GetAxis(0)->SetName("pT_rec");
	129	NSparseHistSt->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
	130	NSparseHistSt->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
	131	NSparseHistSt->GetAxis(2)->SetTitle("ID gen");
	132	NSparseHistSt->GetAxis(2)->SetName("ID_gen");
	133	NSparseHistSt->GetAxis(3)->SetTitle("isph");
	134	NSparseHistSt->GetAxis(3)->SetName("isph");
	135	NSparseHistSt->GetAxis(4)->SetTitle("y");
	136	NSparseHistSt->GetAxis(4)->SetName("y");
eba59d12	137	NSparseHistSt->GetAxis(5)->SetTitle("etaselected");
eba59d12	138	NSparseHistSt->GetAxis(5)->SetName("etaselected");
c683985a	139	fOutput->Add(NSparseHistSt);
	140
	141	//dimensions of THnSparse for the normalization
e19bc133	142	const Int_t nvarev=3;
67b8dfd2	143	// cent Q vec Nch
	144	Int_t binsHistRealEv[nvarev] = { fnCentBins, fnQvecBins, fnNchBins};
	145	Double_t xminHistRealEv[nvarev] = { 0., 0., 0.};
eba59d12	146	Double_t xmaxHistRealEv[nvarev] = { 100., fQvecUpperLim, 2000.};
c683985a	147	THnSparseF* NSparseHistEv = new THnSparseF("NSparseHistEv","NSparseHistEv",nvarev,binsHistRealEv,xminHistRealEv,xmaxHistRealEv);
	148	NSparseHistEv->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
	149	NSparseHistEv->GetAxis(0)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
	150	NSparseHistEv->GetAxis(1)->SetTitle("Q vec");
	151	NSparseHistEv->GetAxis(1)->SetName("Q_vec");
e19bc133	152	NSparseHistEv->GetAxis(2)->SetTitle("N charged");
e19bc133	153	NSparseHistEv->GetAxis(2)->SetName("N_ch");
c683985a	154	fOutput->Add(NSparseHistEv);
	155
	156	PostData(1, fOutput );
	157	PostData(2, fEventCuts);
	158	PostData(3, fTrackCuts);
	159	PostData(4, fHelperPID);
	160	}
	161
	162	//________________________________________________________________________
	163
	164	void AliAnalysisTaskSpectraAllChAOD::UserExec(Option_t *)
	165	{
	166	//Printf("An event");
	167	// main event loop
	168	fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);
	169	if (!fAOD) {
	170	AliWarning("ERROR: AliAODEvent not available \n");
	171	return;
	172	}
	173
	174	if (strcmp(fAOD->ClassName(), "AliAODEvent"))
	175	{
	176	AliFatal("Not processing AODs");
	177	}
	178
	179	if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection
	180
	181	//Default TPC priors
	182	if(fHelperPID->GetPIDType()==kBayes)fHelperPID->GetPIDCombined()->SetDefaultTPCPriors();//FIXME maybe this can go in the UserCreateOutputObject?
	183
	184	Double_t Qvec=0.;//in case of MC we save space in the memory
	185	if(!fIsMC){
eba59d12	186	if(fIsQvecCalibMode){
	187	if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();
	188	else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();
	189	}
	190	else Qvec=fEventCuts->GetQvecPercentile(fVZEROside);
c683985a	191	}
eba59d12	192
c683985a	193	Double_t Cent=fEventCuts->GetCent();
eba59d12	194
c683985a	195	// First do MC to fill up the MC particle array
	196	TClonesArray *arrayMC = 0;
	197	if (fIsMC)
	198	{
	199	arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
	200	if (!arrayMC) {
	201	AliFatal("Error: MC particles branch not found!\n");
	202	}
	203	Int_t nMC = arrayMC->GetEntries();
	204	for (Int_t iMC = 0; iMC < nMC; iMC++)
	205	{
	206	AliAODMCParticle partMC = (AliAODMCParticle) arrayMC->At(iMC);
	207	if(!partMC->Charge()) continue;//Skip neutrals
660a66d4	208	if(fCharge != 0 && partMC->Charge()*fCharge < 0.) continue;//if fCharge != 0 only select fCharge
c683985a	209
eba59d12	210	//flag to select particle in the same ETA acceptance of the tracks (to be used for the comparison with AllCh)
	211	Double_t etaselected=-1.;
	212	if(partMC->Eta()<fTrackCuts->GetEtaMin() \|\| partMC->Eta()>fTrackCuts->GetEtaMax())etaselected=1.;
c683985a	213
c683985a	214	//pt cent IDgen isph y
eba59d12	215	Double_t varSt[6];
c683985a	216	varSt[0]=partMC->Pt();
	217	varSt[1]=Cent;
	218	varSt[2]=(Double_t)fHelperPID->GetParticleSpecies(partMC);
	219	varSt[3]=(Double_t)partMC->IsPhysicalPrimary();
	220	varSt[4]=partMC->Y();
eba59d12	221	varSt[5]=etaselected;
c683985a	222	((THnSparseF*)fOutput->FindObject("NSparseHistSt"))->Fill(varSt);//stack loop
	223
	224	//Printf("a particle");
	225
	226	}
	227	}
	228
	229	//main loop on tracks
e19bc133	230
	231	Int_t Nch = 0.;
	232
c683985a	233	for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
	234	AliAODTrack* track = fAOD->GetTrack(iTracks);
	235	if(fCharge != 0 && track->Charge() != fCharge) continue;//if fCharge != 0 only select fCharge
	236	if (!fTrackCuts->IsSelected(track,kTRUE)) continue; //track selection (rapidity selection NOT in the standard cuts)
dbbc07ac	237	if(!fFillOnlyEvents){
	238	Int_t IDrec=fHelperPID->GetParticleSpecies(track,kTRUE);//id from detector
	239	Double_t y= track->Y(fHelperPID->GetMass((AliHelperParticleSpecies_t)IDrec));
	240	Int_t IDgen=kSpUndefined;//set if MC
	241	Int_t isph=-999;
eba59d12	242	//Int_t iswd=-999;
dbbc07ac	243
	244	if (arrayMC) {
	245	AliAODMCParticle partMC = (AliAODMCParticle) arrayMC->At(TMath::Abs(track->GetLabel()));
	246	if (!partMC) {
	247	AliError("Cannot get MC particle");
	248	continue;
	249	}
	250	IDgen=fHelperPID->GetParticleSpecies(partMC);
	251	isph=partMC->IsPhysicalPrimary();
eba59d12	252	//iswd=partMC->IsSecondaryFromWeakDecay();//FIXME not working on old productions - removed Apr 8th 2014
c683985a	253	}
dbbc07ac	254
eba59d12	255	//pt cent Q vec IDrec IDgen isph y
eba59d12	256	Double_t varTrk[7];
dbbc07ac	257	varTrk[0]=track->Pt();
	258	varTrk[1]=Cent;
	259	varTrk[2]=Qvec;
	260	varTrk[3]=(Double_t)IDrec;
	261	varTrk[4]=(Double_t)IDgen;
	262	varTrk[5]=(Double_t)isph;
eba59d12	263	varTrk[6]=y;
dbbc07ac	264	((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
	265
	266	//for nsigma PID fill double counting of ID
	267	if(fHelperPID->GetPIDType()<kBayes && fDoDoubleCounting){//only nsigma
	268	Bool_t *HasDC;
	269	HasDC=fHelperPID->GetDoubleCounting(track,kTRUE);//get the array with double counting
	270	for(Int_t ipart=0;ipart<kNSpecies;ipart++){
	271	if(HasDC[ipart]==kTRUE){
	272	varTrk[3]=(Double_t)ipart;
	273	((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
	274	}
c683985a	275	}
c683985a	276	}
dbbc07ac	277
	278	//fill all charged (3)
	279	varTrk[3]=3.;
	280	varTrk[4]=3.;
	281	((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
	282	}//end if fFillOnlyEvents
c683985a	283
c683985a	284	//Printf("a track");
dbbc07ac	285
dbbc07ac	286	Nch++;
c683985a	287	} // end loop on tracks
c683985a	288
e19bc133	289	Double_t varEv[3];
c683985a	290	varEv[0]=Cent;
c683985a	291	varEv[1]=Qvec;
e19bc133	292	varEv[2]=Nch;
c683985a	293	((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
	294
	295	PostData(1, fOutput );
	296	PostData(2, fEventCuts);
	297	PostData(3, fTrackCuts);
	298	PostData(4, fHelperPID);
	299	}
	300
	301	//_________________________________________________________________
	302	void AliAnalysisTaskSpectraAllChAOD::Terminate(Option_t *)
	303	{
	304	// Terminate
	305	}