[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),
  fIsAOD160(1),
  fnDCABins(60),
  fDCAmin(-3),
  fDCAmax(3),
  fDCAzCut(999999.)
{
  // 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.,20.,25.,30.,35.,40.,50.,75.,100.};
  const Int_t nptBins=34;
  //const Double_t ptBins[] = {0.20,0.25,0.30,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,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.,20.,25.,30.,35.,40.,50.,75.,100.};
  //const Int_t nptBins=44;
  
  //dimensions of THnSparse for tracks
  const Int_t nvartrk=8;
  //                                             pt          cent          Q vec       IDrec          IDgen        isph         y             DCA 
  Int_t    binsHistRealTrk[nvartrk] = {      nptBins, fnCentBins,     fnQvecBins,        4,             3,          3,        2,        fnDCABins};
  Double_t xminHistRealTrk[nvartrk] = {         0.,          0.,              0.,       -.5,          -0.5,        0.5,      -0.5,          fDCAmin};
  Double_t xmaxHistRealTrk[nvartrk] = {       10.,       100.,     fQvecUpperLim,       3.5,      2.5,        3.5,      0.5,         fDCAmax};    
  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");
  NSparseHistTrk->GetAxis(7)->SetTitle("dca");
  NSparseHistTrk->GetAxis(7)->SetName("dca");
  fOutput->Add(NSparseHistTrk);
  
  //dimensions of THnSparse for stack
  const Int_t nvarst=7;
  //                                             pt          cent    IDgen        isph        y    etaselected       Qvec gen
  Int_t    binsHistRealSt[nvarst] = {      nptBins,  fnCentBins,        3,         2,        2,             1,      fnQvecBins};
  Double_t xminHistRealSt[nvarst] = {         0.,           0.,      -0.5,      -0.5,    -0.5,            0.5,              0.};
  Double_t xmaxHistRealSt[nvarst] = {       10.,        100.,      2.5,       1.5,      0.5,           1.5,      fQvecUpperLim};
  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");
  NSparseHistSt->GetAxis(6)->SetTitle("Q vec gen");
  NSparseHistSt->GetAxis(6)->SetName("Q_vec_gen");
  fOutput->Add(NSparseHistSt);
  
  //dimensions of THnSparse for the normalization
  const Int_t nvarev=4;
  //                                             cent             Q vec                Nch    Qvec gen
  Int_t    binsHistRealEv[nvarev] = {    fnCentBins,      fnQvecBins,           fnNchBins,      fnQvecBins};
  Double_t xminHistRealEv[nvarev] = {           0.,               0.,                   0.,             0.};
  Double_t xmaxHistRealEv[nvarev] = {       100.,  fQvecUpperLim,               2000.,       fQvecUpperLim};
  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");
  NSparseHistEv->GetAxis(3)->SetTitle("Q vec gen");
  NSparseHistEv->GetAxis(3)->SetName("Q_vec_gen");
  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.;
  if(fIsQvecCalibMode){
    if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();
    else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();
  }
  else Qvec=fEventCuts->GetQvecPercentile(fVZEROside);
  
  Double_t QvecMC = 0.;
  if(fIsMC){
    if(fIsQvecCalibMode){
      QvecMC = fEventCuts->CalculateQVectorMC(fVZEROside);
    }
  else QvecMC = fEventCuts->GetQvecPercentileMC(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[7];
	  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;
	  varSt[6]=QvecMC;
	  ((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
	
	if(fIsAOD160){// enabled for new ADO160 only
	  if(partMC->IsSecondaryFromWeakDecay()) isph=2.;
	  if(partMC->IsSecondaryFromMaterial()) isph=3.;

	}
      }
      
      /*** DCA ***/
      Double_t dcaxy = -999.;
      Double_t dcaz = -999.;
      
      Double_t p[2]; 
      if(GetDCA(track,p)){ dcaxy=p[0]; dcaz=p[1];}
      
      if(dcaz >= fDCAzCut) continue;
      
    //pt     cent    Q vec     IDrec     IDgen       isph      y
      Double_t varTrk[8];
      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;
      varTrk[7]=dcaxy;
      ((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[4];
  varEv[0]=Cent;
  varEv[1]=Qvec;
  varEv[2]=Nch;
  varEv[3]=QvecMC;
  ((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
  
  PostData(1, fOutput  );
  PostData(2, fEventCuts);
  PostData(3, fTrackCuts);
  PostData(4, fHelperPID);
}

//_________________________________________________________________
Bool_t  AliAnalysisTaskSpectraAllChAOD::GetDCA(const AliAODTrack* trk, Double_t * p){
  
  //AliAODTrack::DCA(): for newest AOD fTrack->DCA() always gives -999. This should fix.
  //FIXME should update EventCuts?
  //FIXME add track->GetXYZ(p) method
  
  double xyz[3],cov[3];
  
  if (!trk->GetXYZ(xyz)) { // dca is not stored
    AliExternalTrackParam etp;
    etp.CopyFromVTrack(trk);
    AliVEvent* ev = (AliVEvent*)trk->GetEvent();
    if (!ev) {/*printf("Event is not connected to the track\n");*/ return kFALSE;}
    if (!etp.PropagateToDCA(ev->GetPrimaryVertex(), ev->GetMagneticField(),999,xyz,cov)) return kFALSE; // failed, track is too far from vertex
  }
  p[0] = xyz[0];
  p[1] = xyz[1];
  return kTRUE;

}

//_________________________________________________________________
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),
93c24428	69	fQvecUpperLim(100),
	70	fIsAOD160(1),
	71	fnDCABins(60),
	72	fDCAmin(-3),
3dd5ed48	73	fDCAmax(3),
3dd5ed48	74	fDCAzCut(999999.)
c683985a	75	{
	76	// Default constructor
	77	DefineInput(0, TChain::Class());
	78	DefineOutput(1, TList::Class());
	79	DefineOutput(2, AliSpectraAODEventCuts::Class());
	80	DefineOutput(3, AliSpectraAODTrackCuts::Class());
	81	DefineOutput(4, AliHelperPID::Class());
	82	}
	83
	84	//________________________________________________________________________
	85	void AliAnalysisTaskSpectraAllChAOD::UserCreateOutputObjects()
	86	{
	87	// create output objects
	88	fOutput=new TList();
	89	fOutput->SetOwner();
	90	fOutput->SetName("fOutput");
	91
	92	if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");
	93	if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");
	94	if (!fHelperPID) AliFatal("HelperPID should be set in the steering macro");
	95
	96	// binning common to all the THn
c12a5f69	97	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.,20.,25.,30.,35.,40.,50.,75.,100.};
	98	const Int_t nptBins=34;
	99	//const Double_t ptBins[] = {0.20,0.25,0.30,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,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.,20.,25.,30.,35.,40.,50.,75.,100.};
	100	//const Int_t nptBins=44;
c683985a	101
c683985a	102	//dimensions of THnSparse for tracks
3dd5ed48	103	const Int_t nvartrk=8;
1fdecdaf	104	// pt cent Q vec IDrec IDgen isph y DCA
	105	Int_t binsHistRealTrk[nvartrk] = { nptBins, fnCentBins, fnQvecBins, 4, 3, 3, 2, fnDCABins};
	106	Double_t xminHistRealTrk[nvartrk] = { 0., 0., 0., -.5, -0.5, 0.5, -0.5, fDCAmin};
	107	Double_t xmaxHistRealTrk[nvartrk] = { 10., 100., fQvecUpperLim, 3.5, 2.5, 3.5, 0.5, fDCAmax};
c683985a	108	THnSparseF* NSparseHistTrk = new THnSparseF("NSparseHistTrk","NSparseHistTrk",nvartrk,binsHistRealTrk,xminHistRealTrk,xmaxHistRealTrk);
	109	NSparseHistTrk->GetAxis(0)->SetTitle("#it{p}_{T,rec}");
	110	NSparseHistTrk->GetAxis(0)->SetName("pT_rec");
	111	NSparseHistTrk->SetBinEdges(0,ptBins);
	112	NSparseHistTrk->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
	113	NSparseHistTrk->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
	114	NSparseHistTrk->GetAxis(2)->SetTitle("Q vec");
	115	NSparseHistTrk->GetAxis(2)->SetName("Q_vec");
	116	NSparseHistTrk->GetAxis(3)->SetTitle("ID rec");
	117	NSparseHistTrk->GetAxis(3)->SetName("ID_rec");
	118	NSparseHistTrk->GetAxis(4)->SetTitle("ID gen");
	119	NSparseHistTrk->GetAxis(4)->SetName("ID_gen");
	120	NSparseHistTrk->GetAxis(5)->SetTitle("isph");
	121	NSparseHistTrk->GetAxis(5)->SetName("isph");
eba59d12	122	NSparseHistTrk->GetAxis(6)->SetTitle("y");
eba59d12	123	NSparseHistTrk->GetAxis(6)->SetName("y");
93c24428	124	NSparseHistTrk->GetAxis(7)->SetTitle("dca");
93c24428	125	NSparseHistTrk->GetAxis(7)->SetName("dca");
c683985a	126	fOutput->Add(NSparseHistTrk);
	127
	128	//dimensions of THnSparse for stack
ed33b9ef	129	const Int_t nvarst=7;
	130	// pt cent IDgen isph y etaselected Qvec gen
	131	Int_t binsHistRealSt[nvarst] = { nptBins, fnCentBins, 3, 2, 2, 1, fnQvecBins};
	132	Double_t xminHistRealSt[nvarst] = { 0., 0., -0.5, -0.5, -0.5, 0.5, 0.};
	133	Double_t xmaxHistRealSt[nvarst] = { 10., 100., 2.5, 1.5, 0.5, 1.5, fQvecUpperLim};
c683985a	134	THnSparseF* NSparseHistSt = new THnSparseF("NSparseHistSt","NSparseHistSt",nvarst,binsHistRealSt,xminHistRealSt,xmaxHistRealSt);
	135	NSparseHistSt->GetAxis(0)->SetTitle("#it{p}_{T,gen}");
	136	NSparseHistSt->SetBinEdges(0,ptBins);
	137	NSparseHistSt->GetAxis(0)->SetName("pT_rec");
	138	NSparseHistSt->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
	139	NSparseHistSt->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
	140	NSparseHistSt->GetAxis(2)->SetTitle("ID gen");
	141	NSparseHistSt->GetAxis(2)->SetName("ID_gen");
	142	NSparseHistSt->GetAxis(3)->SetTitle("isph");
	143	NSparseHistSt->GetAxis(3)->SetName("isph");
	144	NSparseHistSt->GetAxis(4)->SetTitle("y");
	145	NSparseHistSt->GetAxis(4)->SetName("y");
eba59d12	146	NSparseHistSt->GetAxis(5)->SetTitle("etaselected");
eba59d12	147	NSparseHistSt->GetAxis(5)->SetName("etaselected");
ed33b9ef	148	NSparseHistSt->GetAxis(6)->SetTitle("Q vec gen");
ed33b9ef	149	NSparseHistSt->GetAxis(6)->SetName("Q_vec_gen");
c683985a	150	fOutput->Add(NSparseHistSt);
	151
	152	//dimensions of THnSparse for the normalization
ed33b9ef	153	const Int_t nvarev=4;
	154	// cent Q vec Nch Qvec gen
	155	Int_t binsHistRealEv[nvarev] = { fnCentBins, fnQvecBins, fnNchBins, fnQvecBins};
	156	Double_t xminHistRealEv[nvarev] = { 0., 0., 0., 0.};
	157	Double_t xmaxHistRealEv[nvarev] = { 100., fQvecUpperLim, 2000., fQvecUpperLim};
c683985a	158	THnSparseF* NSparseHistEv = new THnSparseF("NSparseHistEv","NSparseHistEv",nvarev,binsHistRealEv,xminHistRealEv,xmaxHistRealEv);
	159	NSparseHistEv->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
	160	NSparseHistEv->GetAxis(0)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
	161	NSparseHistEv->GetAxis(1)->SetTitle("Q vec");
	162	NSparseHistEv->GetAxis(1)->SetName("Q_vec");
e19bc133	163	NSparseHistEv->GetAxis(2)->SetTitle("N charged");
e19bc133	164	NSparseHistEv->GetAxis(2)->SetName("N_ch");
ed33b9ef	165	NSparseHistEv->GetAxis(3)->SetTitle("Q vec gen");
ed33b9ef	166	NSparseHistEv->GetAxis(3)->SetName("Q_vec_gen");
c683985a	167	fOutput->Add(NSparseHistEv);
	168
	169	PostData(1, fOutput );
	170	PostData(2, fEventCuts);
	171	PostData(3, fTrackCuts);
	172	PostData(4, fHelperPID);
	173	}
	174
	175	//________________________________________________________________________
	176
	177	void AliAnalysisTaskSpectraAllChAOD::UserExec(Option_t *)
	178	{
	179	//Printf("An event");
	180	// main event loop
	181	fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);
	182	if (!fAOD) {
	183	AliWarning("ERROR: AliAODEvent not available \n");
	184	return;
	185	}
	186
	187	if (strcmp(fAOD->ClassName(), "AliAODEvent"))
	188	{
	189	AliFatal("Not processing AODs");
	190	}
	191
	192	if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection
	193
	194	//Default TPC priors
1fdecdaf	195	//if(fHelperPID->GetPIDType()==kBayes)fHelperPID->GetPIDCombined()->SetDefaultTPCPriors();//FIXME maybe this can go in the UserCreateOutputObject?
c683985a	196
ed33b9ef	197	Double_t Qvec=0.;
	198	if(fIsQvecCalibMode){
	199	if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();
	200	else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();
	201	}
	202	else Qvec=fEventCuts->GetQvecPercentile(fVZEROside);
	203
	204	Double_t QvecMC = 0.;
	205	if(fIsMC){
eba59d12	206	if(fIsQvecCalibMode){
ed33b9ef	207	QvecMC = fEventCuts->CalculateQVectorMC(fVZEROside);
eba59d12	208	}
cafe9ad1	209	else QvecMC = fEventCuts->GetQvecPercentileMC(fVZEROside);
c683985a	210	}
eba59d12	211
c683985a	212	Double_t Cent=fEventCuts->GetCent();
eba59d12	213
c683985a	214	// First do MC to fill up the MC particle array
	215	TClonesArray *arrayMC = 0;
	216	if (fIsMC)
	217	{
	218	arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
	219	if (!arrayMC) {
	220	AliFatal("Error: MC particles branch not found!\n");
	221	}
	222	Int_t nMC = arrayMC->GetEntries();
	223	for (Int_t iMC = 0; iMC < nMC; iMC++)
	224	{
	225	AliAODMCParticle partMC = (AliAODMCParticle) arrayMC->At(iMC);
	226	if(!partMC->Charge()) continue;//Skip neutrals
660a66d4	227	if(fCharge != 0 && partMC->Charge()*fCharge < 0.) continue;//if fCharge != 0 only select fCharge
c683985a	228
eba59d12	229	//flag to select particle in the same ETA acceptance of the tracks (to be used for the comparison with AllCh)
eba59d12	230	Double_t etaselected=-1.;
caa224ff	231	if(partMC->Eta()>=fTrackCuts->GetEtaMin() && partMC->Eta()<=fTrackCuts->GetEtaMax())etaselected=1.;
c683985a	232
c683985a	233	//pt cent IDgen isph y
ed33b9ef	234	Double_t varSt[7];
c683985a	235	varSt[0]=partMC->Pt();
	236	varSt[1]=Cent;
	237	varSt[2]=(Double_t)fHelperPID->GetParticleSpecies(partMC);
	238	varSt[3]=(Double_t)partMC->IsPhysicalPrimary();
	239	varSt[4]=partMC->Y();
eba59d12	240	varSt[5]=etaselected;
ed33b9ef	241	varSt[6]=QvecMC;
c683985a	242	((THnSparseF*)fOutput->FindObject("NSparseHistSt"))->Fill(varSt);//stack loop
	243
	244	//Printf("a particle");
	245
	246	}
	247	}
	248
	249	//main loop on tracks
e19bc133	250
	251	Int_t Nch = 0.;
	252
c683985a	253	for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
	254	AliAODTrack* track = fAOD->GetTrack(iTracks);
	255	if(fCharge != 0 && track->Charge() != fCharge) continue;//if fCharge != 0 only select fCharge
	256	if (!fTrackCuts->IsSelected(track,kTRUE)) continue; //track selection (rapidity selection NOT in the standard cuts)
dbbc07ac	257	if(!fFillOnlyEvents){
	258	Int_t IDrec=fHelperPID->GetParticleSpecies(track,kTRUE);//id from detector
	259	Double_t y= track->Y(fHelperPID->GetMass((AliHelperParticleSpecies_t)IDrec));
	260	Int_t IDgen=kSpUndefined;//set if MC
	261	Int_t isph=-999;
93c24428	262	// Int_t iswd=-999;
dbbc07ac	263
	264	if (arrayMC) {
	265	AliAODMCParticle partMC = (AliAODMCParticle) arrayMC->At(TMath::Abs(track->GetLabel()));
	266	if (!partMC) {
	267	AliError("Cannot get MC particle");
	268	continue;
	269	}
	270	IDgen=fHelperPID->GetParticleSpecies(partMC);
	271	isph=partMC->IsPhysicalPrimary();
eba59d12	272	//iswd=partMC->IsSecondaryFromWeakDecay();//FIXME not working on old productions - removed Apr 8th 2014
93c24428	273
93c24428	274	if(fIsAOD160){// enabled for new ADO160 only
3dd5ed48	275	if(partMC->IsSecondaryFromWeakDecay()) isph=2.;
3dd5ed48	276	if(partMC->IsSecondaryFromMaterial()) isph=3.;
93c24428	277
93c24428	278	}
c683985a	279	}
dbbc07ac	280
93c24428	281	/* DCA */
93c24428	282	Double_t dcaxy = -999.;
3dd5ed48	283	Double_t dcaz = -999.;
93c24428	284
665f3db3	285	Double_t p[2];
3dd5ed48	286	if(GetDCA(track,p)){ dcaxy=p[0]; dcaz=p[1];}
	287
	288	if(dcaz >= fDCAzCut) continue;
93c24428	289
eba59d12	290	//pt cent Q vec IDrec IDgen isph y
3dd5ed48	291	Double_t varTrk[8];
dbbc07ac	292	varTrk[0]=track->Pt();
	293	varTrk[1]=Cent;
	294	varTrk[2]=Qvec;
	295	varTrk[3]=(Double_t)IDrec;
	296	varTrk[4]=(Double_t)IDgen;
	297	varTrk[5]=(Double_t)isph;
eba59d12	298	varTrk[6]=y;
93c24428	299	varTrk[7]=dcaxy;
dbbc07ac	300	((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
	301
	302	//for nsigma PID fill double counting of ID
	303	if(fHelperPID->GetPIDType()<kBayes && fDoDoubleCounting){//only nsigma
	304	Bool_t *HasDC;
	305	HasDC=fHelperPID->GetDoubleCounting(track,kTRUE);//get the array with double counting
	306	for(Int_t ipart=0;ipart<kNSpecies;ipart++){
	307	if(HasDC[ipart]==kTRUE){
	308	varTrk[3]=(Double_t)ipart;
	309	((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
	310	}
c683985a	311	}
c683985a	312	}
dbbc07ac	313
	314	//fill all charged (3)
	315	varTrk[3]=3.;
	316	varTrk[4]=3.;
	317	((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
	318	}//end if fFillOnlyEvents
c683985a	319
c683985a	320	//Printf("a track");
dbbc07ac	321
dbbc07ac	322	Nch++;
c683985a	323	} // end loop on tracks
c683985a	324
ed33b9ef	325	Double_t varEv[4];
c683985a	326	varEv[0]=Cent;
c683985a	327	varEv[1]=Qvec;
e19bc133	328	varEv[2]=Nch;
ed33b9ef	329	varEv[3]=QvecMC;
c683985a	330	((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
	331
	332	PostData(1, fOutput );
	333	PostData(2, fEventCuts);
	334	PostData(3, fTrackCuts);
	335	PostData(4, fHelperPID);
	336	}
	337
665f3db3	338	//_________________________________________________________________
	339	Bool_t AliAnalysisTaskSpectraAllChAOD::GetDCA(const AliAODTrack* trk, Double_t * p){
	340
	341	//AliAODTrack::DCA(): for newest AOD fTrack->DCA() always gives -999. This should fix.
	342	//FIXME should update EventCuts?
	343	//FIXME add track->GetXYZ(p) method
	344
	345	double xyz[3],cov[3];
	346
	347	if (!trk->GetXYZ(xyz)) { // dca is not stored
	348	AliExternalTrackParam etp;
	349	etp.CopyFromVTrack(trk);
	350	AliVEvent* ev = (AliVEvent*)trk->GetEvent();
	351	if (!ev) {/printf("Event is not connected to the track\n");/ return kFALSE;}
	352	if (!etp.PropagateToDCA(ev->GetPrimaryVertex(), ev->GetMagneticField(),999,xyz,cov)) return kFALSE; // failed, track is too far from vertex
	353	}
	354	p[0] = xyz[0];
	355	p[1] = xyz[1];
	356	return kTRUE;
	357
	358	}
	359
c683985a	360	//_________________________________________________________________
	361	void AliAnalysisTaskSpectraAllChAOD::Terminate(Option_t *)
	362	{
	363	// Terminate
	364	}