[u/mrichter/AliRoot.git] / PWG / DevNanoAOD / AliAnalysisTaskSpectraAllChNanoAOD.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.                  *
 **************************************************************************/

//-----------------------------------------------------------------
//         AliAnalysisTaskSpectraAllChNanoAOD class
// NanoAOD edition: this is only meant to test the developement
// version of NanoAODs
//
//-----------------------------------------------------------------

#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 "AliAnalysisTaskSpectraAllChNanoAOD.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>
#include "AliNanoAODHeader.h"
#include "AliNanoAODTrack.h"

using namespace AliHelperPIDNameSpace;
using namespace std;

ClassImp(AliAnalysisTaskSpectraAllChNanoAOD)

//________________________________________________________________________
AliAnalysisTaskSpectraAllChNanoAOD::AliAnalysisTaskSpectraAllChNanoAOD(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(40),
  fnNchBins(200)
{
  // Default constructor
  DefineInput(0, TChain::Class());
  DefineOutput(1, TList::Class());
  DefineOutput(2, AliSpectraAODEventCuts::Class());
  DefineOutput(3, AliSpectraAODTrackCuts::Class());
  DefineOutput(4, AliHelperPID::Class());
}

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

//________________________________________________________________________

void AliAnalysisTaskSpectraAllChNanoAOD::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");
    } 

  AliNanoAODHeader * headNano = dynamic_cast<AliNanoAODHeader*>((TObject*)fAOD->GetHeader());
  
  Bool_t isNano = (headNano != 0);
 
  if(!isNano) {
    if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection
  }
  //Default TPC priors
  if(fHelperPID->GetPIDType()==kBayes)fHelperPID->GetPIDCombined()->SetDefaultTPCPriors();//FIXME maybe this can go in the UserCreateOutputObject?
  
  Double_t Qvec=0.;//in case of MC we save space in the memory
  if(!fIsMC){
    if(!isNano) {
      if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();
      else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();
    } else {
      
      if(fVZEROside==0)      Qvec=headNano->GetVar(0); // FIXME: we need proper getters here
      else if (fVZEROside==1)Qvec=headNano->GetVar(1);
      
    }
  }

  Double_t Cent=isNano ? headNano->GetVar(2) : 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
	  
	  if(partMC->Eta() < fTrackCuts->GetEtaMin() || partMC->Eta() > fTrackCuts->GetEtaMax())continue;//ETA CUT ON GENERATED!!!!!!!!!!!!!!!!!!!!!!!!!!
	  
	  //pt     cent    IDgen        isph        y
 	  Double_t varSt[5];
	  varSt[0]=partMC->Pt();
	  varSt[1]=Cent;
	  varSt[2]=(Double_t)fHelperPID->GetParticleSpecies(partMC);
	  varSt[3]=(Double_t)partMC->IsPhysicalPrimary();
	  varSt[4]=partMC->Y();
	  ((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++) {
    AliVTrack* track = (AliVTrack*) fAOD->GetTrack(iTracks);
    if(fCharge != 0 && track->Charge() != fCharge) continue;//if fCharge != 0 only select fCharge 
    if(!isNano) {
      if (!fTrackCuts->IsSelected((AliAODTrack*)track,kTRUE)) continue; //track selection (rapidity selection NOT in the standard cuts)
    }
    
    if(!fFillOnlyEvents){
      Int_t IDrec=isNano ? GetNanoTrackID (track) : fHelperPID->GetParticleSpecies(track,kTRUE);//id from detector      
      Double_t y= 0;
      if(isNano) y = ((AliNanoAODTrack*)track)->Y(fHelperPID->GetMass((AliHelperParticleSpecies_t)IDrec));
      else y = ((AliAODTrack*)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
      }
      
    //pt     cent    Q vec     IDrec     IDgen       isph           iswd      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]=(Double_t)iswd;
      varTrk[7]=y;
      ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
      
      //for nsigma PID fill double counting of ID
      if(fHelperPID->GetPIDType()<kBayes && fDoDoubleCounting){//only nsigma
	Bool_t *HasDC;
	HasDC=fHelperPID->GetDoubleCounting(track,kTRUE);//get the array with double counting
	for(Int_t ipart=0;ipart<kNSpecies;ipart++){
	  if(HasDC[ipart]==kTRUE){
	    varTrk[3]=(Double_t)ipart;
	    ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
	  }
	}
      }
      
      //fill all charged (3)
      varTrk[3]=3.;
      varTrk[4]=3.;
      ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
    }//end if fFillOnlyEvents
    
    //Printf("a track");
      
    Nch++;
  } // end loop on tracks
  
  Double_t varEv[3];
  varEv[0]=Cent;
  varEv[1]=Qvec;
  varEv[2]=Nch;
  ((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
  
  PostData(1, fOutput  );
  PostData(2, fEventCuts);
  PostData(3, fTrackCuts);
  PostData(4, fHelperPID);
}

//_________________________________________________________________
void   AliAnalysisTaskSpectraAllChNanoAOD::Terminate(Option_t *)
{
  // Terminate
}


Int_t AliAnalysisTaskSpectraAllChNanoAOD::GetNanoTrackID(AliVTrack * track) {
  // Applies nsigma PID to nano tracks
  AliNanoAODTrack * nanoTrack = dynamic_cast<AliNanoAODTrack*>(track);
  if(!nanoTrack) AliFatal("Not a nano AOD track");

  // Cache indexes
  static const Int_t kcstNSigmaTPCPi  = AliNanoAODTrackMapping::GetInstance()->GetVarIndex("cstNSigmaTPCPi");
  static const Int_t kcstNSigmaTPCKa  = AliNanoAODTrackMapping::GetInstance()->GetVarIndex("cstNSigmaTPCKa");
  static const Int_t kcstNSigmaTPCPr  = AliNanoAODTrackMapping::GetInstance()->GetVarIndex("cstNSigmaTPCPr");
  static const Int_t kcstNSigmaTOFPi  = AliNanoAODTrackMapping::GetInstance()->GetVarIndex("cstNSigmaTOFPi");
  static const Int_t kcstNSigmaTOFKa  = AliNanoAODTrackMapping::GetInstance()->GetVarIndex("cstNSigmaTOFKa");
  static const Int_t kcstNSigmaTOFPr  = AliNanoAODTrackMapping::GetInstance()->GetVarIndex("cstNSigmaTOFPr");

  Double_t nSigmaPID = 3.0;


  //get the identity of the particle with the minimum Nsigma
  Double_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
  if(nanoTrack->Pt() > fTrackCuts->GetPtTOFMatching()) {
    nsigmaProton =  TMath::Sqrt(nanoTrack->GetVar(kcstNSigmaTPCPr)*nanoTrack->GetVar(kcstNSigmaTPCPr)+nanoTrack->GetVar(kcstNSigmaTOFPr)*nanoTrack->GetVar(kcstNSigmaTOFPr));
    nsigmaKaon   =  TMath::Sqrt(nanoTrack->GetVar(kcstNSigmaTPCKa)*nanoTrack->GetVar(kcstNSigmaTPCKa)+nanoTrack->GetVar(kcstNSigmaTOFKa)*nanoTrack->GetVar(kcstNSigmaTOFKa));
    nsigmaPion   =  TMath::Sqrt(nanoTrack->GetVar(kcstNSigmaTPCPi)*nanoTrack->GetVar(kcstNSigmaTPCPi)+nanoTrack->GetVar(kcstNSigmaTOFPi)*nanoTrack->GetVar(kcstNSigmaTOFPi));
  }
  else {
    nsigmaProton =  TMath::Abs(nanoTrack->GetVar(kcstNSigmaTPCPr));
    nsigmaKaon   =  TMath::Abs(nanoTrack->GetVar(kcstNSigmaTPCKa));  
    nsigmaPion   =  TMath::Abs(nanoTrack->GetVar(kcstNSigmaTPCPi));  
  }

  
  // guess the particle based on the smaller nsigma (within nSigmaPID)
  if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton )) return kSpUndefined;//if is the default value for the three
  
  if( ( nsigmaKaon   < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon   < nSigmaPID)){
    return kSpKaon;
  }
  if( ( nsigmaPion   < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion   < nSigmaPID)){    
    return kSpPion;
  }
  if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < nSigmaPID)){
    return kSpProton;
  }
  // else, return undefined
  return kSpUndefined;
  

}
Commit	Line	Data
9d6172fe	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	// AliAnalysisTaskSpectraAllChNanoAOD class
	18	// NanoAOD edition: this is only meant to test the developement
	19	// version of NanoAODs
	20	//
	21	//-----------------------------------------------------------------
	22
	23	#include "TChain.h"
	24	#include "TTree.h"
	25	#include "TLegend.h"
	26	#include "TH1F.h"
	27	#include "TH2F.h"
	28	#include "THnSparse.h"
	29	#include "TCanvas.h"
	30	#include "AliAnalysisTask.h"
	31	#include "AliAODTrack.h"
	32	#include "AliAODMCParticle.h"
	33	#include "AliVParticle.h"
	34	#include "AliAODEvent.h"
	35	#include "AliAODInputHandler.h"
	36	#include "AliAnalysisTaskSpectraAllChNanoAOD.h"
	37	#include "AliAnalysisTaskESDfilter.h"
	38	#include "AliAnalysisDataContainer.h"
	39	#include "AliSpectraAODTrackCuts.h"
	40	#include "AliSpectraAODEventCuts.h"
	41	#include "AliHelperPID.h"
	42	#include "AliPIDCombined.h"
	43	#include "AliCentrality.h"
	44	#include "TProof.h"
	45	#include "AliVEvent.h"
	46	#include "AliStack.h"
	47	#include <TMCProcess.h>
	48
	49	#include <iostream>
	50	#include "AliNanoAODHeader.h"
	51	#include "AliNanoAODTrack.h"
	52
	53	using namespace AliHelperPIDNameSpace;
	54	using namespace std;
	55
	56	ClassImp(AliAnalysisTaskSpectraAllChNanoAOD)
	57
	58	//________________________________________________________________________
	59	AliAnalysisTaskSpectraAllChNanoAOD::AliAnalysisTaskSpectraAllChNanoAOD(const char *name) : AliAnalysisTaskSE(name),
	60	fAOD(0x0),
	61	fTrackCuts(0x0),
	62	fEventCuts(0x0),
	63	fHelperPID(0x0),
	64	fIsMC(0),
65	fDoDoubleCounting(0),
66	fFillOnlyEvents(0),
67	fCharge(0),
68	fVZEROside(0),
69	fOutput(0x0),
70	fnCentBins(20),
71	fnQvecBins(40),
72	fnNchBins(200)
73	{
74	// Default constructor
75	DefineInput(0, TChain::Class());
76	DefineOutput(1, TList::Class());
77	DefineOutput(2, AliSpectraAODEventCuts::Class());
78	DefineOutput(3, AliSpectraAODTrackCuts::Class());
79	DefineOutput(4, AliHelperPID::Class());
80	}
81
82	//________________________________________________________________________
83	void AliAnalysisTaskSpectraAllChNanoAOD::UserCreateOutputObjects()
84	{
85	// create output objects
86	fOutput=new TList();
87	fOutput->SetOwner();
88	fOutput->SetName("fOutput");
89
90	if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");
91	if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");
92	if (!fHelperPID) AliFatal("HelperPID should be set in the steering macro");
93
94	// binning common to all the THn
95	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.};
96	const Int_t nptBins=26;
97
98	//dimensions of THnSparse for tracks
99	const Int_t nvartrk=8;
100	// pt cent Q vec IDrec IDgen isph iswd y
101	Int_t binsHistRealTrk[nvartrk] = { nptBins, fnCentBins, fnQvecBins, 4, 3, 2, 2, 2};
102	Double_t xminHistRealTrk[nvartrk] = { 0., 0., 0., -.5, -0.5, -0.5, -0.5, -0.5};
103	Double_t xmaxHistRealTrk[nvartrk] = { 10., 100., 8., 3.5, 2.5, 1.5, 1.5, 0.5};
104	THnSparseF* NSparseHistTrk = new THnSparseF("NSparseHistTrk","NSparseHistTrk",nvartrk,binsHistRealTrk,xminHistRealTrk,xmaxHistRealTrk);
105	NSparseHistTrk->GetAxis(0)->SetTitle("#it{p}_{T,rec}");
106	NSparseHistTrk->GetAxis(0)->SetName("pT_rec");
107	NSparseHistTrk->SetBinEdges(0,ptBins);
108	NSparseHistTrk->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
109	NSparseHistTrk->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
110	NSparseHistTrk->GetAxis(2)->SetTitle("Q vec");
111	NSparseHistTrk->GetAxis(2)->SetName("Q_vec");
112	NSparseHistTrk->GetAxis(3)->SetTitle("ID rec");
113	NSparseHistTrk->GetAxis(3)->SetName("ID_rec");
114	NSparseHistTrk->GetAxis(4)->SetTitle("ID gen");
115	NSparseHistTrk->GetAxis(4)->SetName("ID_gen");
116	NSparseHistTrk->GetAxis(5)->SetTitle("isph");
117	NSparseHistTrk->GetAxis(5)->SetName("isph");
118	NSparseHistTrk->GetAxis(6)->SetTitle("iswd");
119	NSparseHistTrk->GetAxis(6)->SetName("iswd");
120	NSparseHistTrk->GetAxis(7)->SetTitle("y");
121	NSparseHistTrk->GetAxis(7)->SetName("y");
122	fOutput->Add(NSparseHistTrk);
123
124	//dimensions of THnSparse for stack
125	const Int_t nvarst=5;
126	// pt cent IDgen isph y
127	Int_t binsHistRealSt[nvarst] = { nptBins, fnCentBins, 3, 2, 2};
128	Double_t xminHistRealSt[nvarst] = { 0., 0., -0.5, -0.5, -0.5};
129	Double_t xmaxHistRealSt[nvarst] = { 10., 100., 2.5, 1.5, 0.5};
130	THnSparseF* NSparseHistSt = new THnSparseF("NSparseHistSt","NSparseHistSt",nvarst,binsHistRealSt,xminHistRealSt,xmaxHistRealSt);
131	NSparseHistSt->GetAxis(0)->SetTitle("#it{p}_{T,gen}");
132	NSparseHistSt->SetBinEdges(0,ptBins);
133	NSparseHistSt->GetAxis(0)->SetName("pT_rec");
134	NSparseHistSt->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
135	NSparseHistSt->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
136	NSparseHistSt->GetAxis(2)->SetTitle("ID gen");
137	NSparseHistSt->GetAxis(2)->SetName("ID_gen");
138	NSparseHistSt->GetAxis(3)->SetTitle("isph");
139	NSparseHistSt->GetAxis(3)->SetName("isph");
140	NSparseHistSt->GetAxis(4)->SetTitle("y");
141	NSparseHistSt->GetAxis(4)->SetName("y");
142	fOutput->Add(NSparseHistSt);
143
144	//dimensions of THnSparse for the normalization
145	const Int_t nvarev=3;
146	// cent Q vec Nch
147	Int_t binsHistRealEv[nvarev] = { fnCentBins, fnQvecBins, fnNchBins};
148	Double_t xminHistRealEv[nvarev] = { 0., 0., 0.};
149	Double_t xmaxHistRealEv[nvarev] = { 100., 8., 2000.};
150	THnSparseF* NSparseHistEv = new THnSparseF("NSparseHistEv","NSparseHistEv",nvarev,binsHistRealEv,xminHistRealEv,xmaxHistRealEv);
151	NSparseHistEv->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
152	NSparseHistEv->GetAxis(0)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
153	NSparseHistEv->GetAxis(1)->SetTitle("Q vec");
154	NSparseHistEv->GetAxis(1)->SetName("Q_vec");
155	NSparseHistEv->GetAxis(2)->SetTitle("N charged");
156	NSparseHistEv->GetAxis(2)->SetName("N_ch");
157	fOutput->Add(NSparseHistEv);
158
159	PostData(1, fOutput );
160	PostData(2, fEventCuts);
161	PostData(3, fTrackCuts);
162	PostData(4, fHelperPID);
163	}
164
165	//________________________________________________________________________
166
167	void AliAnalysisTaskSpectraAllChNanoAOD::UserExec(Option_t *)
168	{
169	//Printf("An event");
170	// main event loop
171	fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);
172	if (!fAOD) {
173	AliWarning("ERROR: AliAODEvent not available \n");
174	return;
175	}
176
177	if (strcmp(fAOD->ClassName(), "AliAODEvent"))
178	{
179	AliFatal("Not processing AODs");
180	}
181
182	AliNanoAODHeader * headNano = dynamic_cast<AliNanoAODHeader>((TObject)fAOD->GetHeader());
183
184	Bool_t isNano = (headNano != 0);
185
186	if(!isNano) {
187	if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection
188	}
189	//Default TPC priors
190	if(fHelperPID->GetPIDType()==kBayes)fHelperPID->GetPIDCombined()->SetDefaultTPCPriors();//FIXME maybe this can go in the UserCreateOutputObject?
191
192	Double_t Qvec=0.;//in case of MC we save space in the memory
193	if(!fIsMC){
194	if(!isNano) {
195	if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();
196	else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();
197	} else {
198
199	if(fVZEROside==0) Qvec=headNano->GetVar(0); // FIXME: we need proper getters here
200	else if (fVZEROside==1)Qvec=headNano->GetVar(1);
201
202	}
203	}
204
205	Double_t Cent=isNano ? headNano->GetVar(2) : fEventCuts->GetCent();
206
207	// First do MC to fill up the MC particle array
208	TClonesArray *arrayMC = 0;
209	if (fIsMC)
210	{
211	arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
212	if (!arrayMC) {
213	AliFatal("Error: MC particles branch not found!\n");
214	}
215	Int_t nMC = arrayMC->GetEntries();
216	for (Int_t iMC = 0; iMC < nMC; iMC++)
217	{
218	AliAODMCParticle partMC = (AliAODMCParticle) arrayMC->At(iMC);
219	if(!partMC->Charge()) continue;//Skip neutrals
220	if(fCharge != 0 && partMC->Charge()*fCharge < 0.) continue;//if fCharge != 0 only select fCharge
221
222	if(partMC->Eta() < fTrackCuts->GetEtaMin() \|\| partMC->Eta() > fTrackCuts->GetEtaMax())continue;//ETA CUT ON GENERATED!!!!!!!!!!!!!!!!!!!!!!!!!!
223
224	//pt cent IDgen isph y
225	Double_t varSt[5];
226	varSt[0]=partMC->Pt();
227	varSt[1]=Cent;
228	varSt[2]=(Double_t)fHelperPID->GetParticleSpecies(partMC);
229	varSt[3]=(Double_t)partMC->IsPhysicalPrimary();
230	varSt[4]=partMC->Y();
231	((THnSparseF*)fOutput->FindObject("NSparseHistSt"))->Fill(varSt);//stack loop
232
233	//Printf("a particle");
234
235	}
236	}
237
238	//main loop on tracks
239
240	Int_t Nch = 0.;
241
242	for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
243	AliVTrack* track = (AliVTrack*) fAOD->GetTrack(iTracks);
244	if(fCharge != 0 && track->Charge() != fCharge) continue;//if fCharge != 0 only select fCharge
245	if(!isNano) {
246	if (!fTrackCuts->IsSelected((AliAODTrack*)track,kTRUE)) continue; //track selection (rapidity selection NOT in the standard cuts)
247	}
248
249	if(!fFillOnlyEvents){
250	Int_t IDrec=isNano ? GetNanoTrackID (track) : fHelperPID->GetParticleSpecies(track,kTRUE);//id from detector
251	Double_t y= 0;
252	if(isNano) y = ((AliNanoAODTrack*)track)->Y(fHelperPID->GetMass((AliHelperParticleSpecies_t)IDrec));
253	else y = ((AliAODTrack*)track)->Y(fHelperPID->GetMass((AliHelperParticleSpecies_t)IDrec));
254	Int_t IDgen=kSpUndefined;//set if MC
255	Int_t isph=-999;
256	Int_t iswd=-999;
257
258	if (arrayMC) {
259	AliAODMCParticle partMC = (AliAODMCParticle) arrayMC->At(TMath::Abs(track->GetLabel()));
260	if (!partMC) {
261	AliError("Cannot get MC particle");
262	continue;
263	}
264	IDgen=fHelperPID->GetParticleSpecies(partMC);
265	isph=partMC->IsPhysicalPrimary();
266	iswd=partMC->IsSecondaryFromWeakDecay();//FIXME not working on old productions
267	}
268
269	//pt cent Q vec IDrec IDgen isph iswd y
270	Double_t varTrk[8];
271	varTrk[0]=track->Pt();
272	varTrk[1]=Cent;
273	varTrk[2]=Qvec;
274	varTrk[3]=(Double_t)IDrec;
275	varTrk[4]=(Double_t)IDgen;
276	varTrk[5]=(Double_t)isph;
277	varTrk[6]=(Double_t)iswd;
278	varTrk[7]=y;
279	((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
280
281	//for nsigma PID fill double counting of ID
282	if(fHelperPID->GetPIDType()<kBayes && fDoDoubleCounting){//only nsigma
283	Bool_t *HasDC;
284	HasDC=fHelperPID->GetDoubleCounting(track,kTRUE);//get the array with double counting
285	for(Int_t ipart=0;ipart<kNSpecies;ipart++){
286	if(HasDC[ipart]==kTRUE){
287	varTrk[3]=(Double_t)ipart;
288	((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
289	}
290	}
291	}
292
293	//fill all charged (3)
294	varTrk[3]=3.;
295	varTrk[4]=3.;
296	((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
297	}//end if fFillOnlyEvents
298
299	//Printf("a track");
300
301	Nch++;
302	} // end loop on tracks
303
304	Double_t varEv[3];
305	varEv[0]=Cent;
306	varEv[1]=Qvec;
307	varEv[2]=Nch;
308	((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
309
310	PostData(1, fOutput );
311	PostData(2, fEventCuts);
312	PostData(3, fTrackCuts);
313	PostData(4, fHelperPID);
314	}
315
316	//_________________________________________________________________
317	void AliAnalysisTaskSpectraAllChNanoAOD::Terminate(Option_t *)
318	{
319	// Terminate
320	}
321
322
323	Int_t AliAnalysisTaskSpectraAllChNanoAOD::GetNanoTrackID(AliVTrack * track) {
324	// Applies nsigma PID to nano tracks
325	AliNanoAODTrack * nanoTrack = dynamic_cast<AliNanoAODTrack*>(track);
326	if(!nanoTrack) AliFatal("Not a nano AOD track");
327
328	// Cache indexes
329	static const Int_t kcstNSigmaTPCPi = AliNanoAODTrackMapping::GetInstance()->GetVarIndex("cstNSigmaTPCPi");
330	static const Int_t kcstNSigmaTPCKa = AliNanoAODTrackMapping::GetInstance()->GetVarIndex("cstNSigmaTPCKa");
331	static const Int_t kcstNSigmaTPCPr = AliNanoAODTrackMapping::GetInstance()->GetVarIndex("cstNSigmaTPCPr");
332	static const Int_t kcstNSigmaTOFPi = AliNanoAODTrackMapping::GetInstance()->GetVarIndex("cstNSigmaTOFPi");
333	static const Int_t kcstNSigmaTOFKa = AliNanoAODTrackMapping::GetInstance()->GetVarIndex("cstNSigmaTOFKa");
334	static const Int_t kcstNSigmaTOFPr = AliNanoAODTrackMapping::GetInstance()->GetVarIndex("cstNSigmaTOFPr");
335
336	Double_t nSigmaPID = 3.0;
337
338
339
340	//get the identity of the particle with the minimum Nsigma
341	Double_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
342	if(nanoTrack->Pt() > fTrackCuts->GetPtTOFMatching()) {
343	nsigmaProton = TMath::Sqrt(nanoTrack->GetVar(kcstNSigmaTPCPr)nanoTrack->GetVar(kcstNSigmaTPCPr)+nanoTrack->GetVar(kcstNSigmaTOFPr)nanoTrack->GetVar(kcstNSigmaTOFPr));
344	nsigmaKaon = TMath::Sqrt(nanoTrack->GetVar(kcstNSigmaTPCKa)nanoTrack->GetVar(kcstNSigmaTPCKa)+nanoTrack->GetVar(kcstNSigmaTOFKa)nanoTrack->GetVar(kcstNSigmaTOFKa));
345	nsigmaPion = TMath::Sqrt(nanoTrack->GetVar(kcstNSigmaTPCPi)nanoTrack->GetVar(kcstNSigmaTPCPi)+nanoTrack->GetVar(kcstNSigmaTOFPi)nanoTrack->GetVar(kcstNSigmaTOFPi));
346	}
347	else {
348	nsigmaProton = TMath::Abs(nanoTrack->GetVar(kcstNSigmaTPCPr));
349	nsigmaKaon = TMath::Abs(nanoTrack->GetVar(kcstNSigmaTPCKa));
350	nsigmaPion = TMath::Abs(nanoTrack->GetVar(kcstNSigmaTPCPi));
351	}
352
353
354
355	// guess the particle based on the smaller nsigma (within nSigmaPID)
356	if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton )) return kSpUndefined;//if is the default value for the three
357
358	if( ( nsigmaKaon < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon < nSigmaPID)){
359	return kSpKaon;
360	}
361	if( ( nsigmaPion < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion < nSigmaPID)){
362	return kSpPion;
363	}
364	if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < nSigmaPID)){
365	return kSpProton;
366	}
367	// else, return undefined
368	return kSpUndefined;
369
370
371	}