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