[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;

      //if the q vector is done using the TPC, we avoid overlap
      if (fVZEROside==2 && TMath::Abs(track->Eta())<0.5)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;
896040ab	298
	299	//if the q vector is done using the TPC, we avoid overlap
	300	if (fVZEROside==2 && TMath::Abs(track->Eta())<0.5)continue;
	301
eba59d12	302	//pt cent Q vec IDrec IDgen isph y
3dd5ed48	303	Double_t varTrk[8];
dbbc07ac	304	varTrk[0]=track->Pt();
dbbc07ac	305	varTrk[1]=Cent;
09f5725e	306	if(fIsMC && fQtrk==1) varTrk[2]=QvecMC;
388cfa2a	307	else varTrk[2]=Qvec;
dbbc07ac	308	varTrk[3]=(Double_t)IDrec;
	309	varTrk[4]=(Double_t)IDgen;
	310	varTrk[5]=(Double_t)isph;
eba59d12	311	varTrk[6]=y;
93c24428	312	varTrk[7]=dcaxy;
dbbc07ac	313	((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
	314
	315	//for nsigma PID fill double counting of ID
	316	if(fHelperPID->GetPIDType()<kBayes && fDoDoubleCounting){//only nsigma
	317	Bool_t *HasDC;
	318	HasDC=fHelperPID->GetDoubleCounting(track,kTRUE);//get the array with double counting
	319	for(Int_t ipart=0;ipart<kNSpecies;ipart++){
	320	if(HasDC[ipart]==kTRUE){
	321	varTrk[3]=(Double_t)ipart;
	322	((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
	323	}
c683985a	324	}
c683985a	325	}
dbbc07ac	326
	327	//fill all charged (3)
	328	varTrk[3]=3.;
	329	varTrk[4]=3.;
	330	((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
	331	}//end if fFillOnlyEvents
c683985a	332
c683985a	333	//Printf("a track");
dbbc07ac	334
dbbc07ac	335	Nch++;
c683985a	336	} // end loop on tracks
c683985a	337
ed33b9ef	338	Double_t varEv[4];
c683985a	339	varEv[0]=Cent;
c683985a	340	varEv[1]=Qvec;
e19bc133	341	varEv[2]=Nch;
ed33b9ef	342	varEv[3]=QvecMC;
c683985a	343	((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
	344
	345	PostData(1, fOutput );
	346	PostData(2, fEventCuts);
	347	PostData(3, fTrackCuts);
	348	PostData(4, fHelperPID);
	349	}
	350
665f3db3	351	//_________________________________________________________________
	352	Bool_t AliAnalysisTaskSpectraAllChAOD::GetDCA(const AliAODTrack* trk, Double_t * p){
	353
	354	//AliAODTrack::DCA(): for newest AOD fTrack->DCA() always gives -999. This should fix.
	355	//FIXME should update EventCuts?
	356	//FIXME add track->GetXYZ(p) method
	357
	358	double xyz[3],cov[3];
	359
	360	if (!trk->GetXYZ(xyz)) { // dca is not stored
	361	AliExternalTrackParam etp;
	362	etp.CopyFromVTrack(trk);
	363	AliVEvent* ev = (AliVEvent*)trk->GetEvent();
	364	if (!ev) {/printf("Event is not connected to the track\n");/ return kFALSE;}
	365	if (!etp.PropagateToDCA(ev->GetPrimaryVertex(), ev->GetMagneticField(),999,xyz,cov)) return kFALSE; // failed, track is too far from vertex
	366	}
	367	p[0] = xyz[0];
	368	p[1] = xyz[1];
	369	return kTRUE;
	370
	371	}
	372
c683985a	373	//_________________________________________________________________
	374	void AliAnalysisTaskSpectraAllChAOD::Terminate(Option_t *)
	375	{
	376	// Terminate
	377	}