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