#include "AliAnalysisDataContainer.h"\r
#include "AliSpectraAODTrackCuts.h"\r
#include "AliSpectraAODEventCuts.h"\r
+#include "AliHelperPID.h"\r
+#include "AliPIDCombined.h"\r
#include "AliCentrality.h"\r
#include "TProof.h"\r
#include "AliVEvent.h"\r
\r
#include <iostream>\r
\r
+using namespace AliHelperPIDNameSpace;\r
using namespace std;\r
\r
ClassImp(AliAnalysisTaskSpectraAllChAOD)\r
\r
//________________________________________________________________________\r
-AliAnalysisTaskSpectraAllChAOD::AliAnalysisTaskSpectraAllChAOD(const char *name) : AliAnalysisTaskSE(name), fAOD(0), fTrackCuts(0), fEventCuts(0), fIsMC(0), fOutput(0)\r
+AliAnalysisTaskSpectraAllChAOD::AliAnalysisTaskSpectraAllChAOD(const char *name) : AliAnalysisTaskSE(name), \r
+ fAOD(0x0),\r
+ fTrackCuts(0x0),\r
+ fEventCuts(0x0),\r
+ fHelperPID(0x0),\r
+ fIsMC(0),\r
+ fCharge(0),\r
+ fVZEROside(0),\r
+ fOutput(0x0),\r
+ fnCentBins(20),\r
+ fnQvecBins(50)\r
{\r
// Default constructor\r
- \r
DefineInput(0, TChain::Class());\r
DefineOutput(1, TList::Class());\r
DefineOutput(2, AliSpectraAODEventCuts::Class());\r
DefineOutput(3, AliSpectraAODTrackCuts::Class());\r
- \r
+ DefineOutput(4, AliHelperPID::Class());\r
}\r
-//________________________________________________________________________\r
+\r
//________________________________________________________________________\r
void AliAnalysisTaskSpectraAllChAOD::UserCreateOutputObjects()\r
{\r
\r
if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");\r
if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");\r
+ if (!fHelperPID) AliFatal("HelperPID should be set in the steering macro");\r
\r
- //dimensions of standard THnSparse\r
- const Int_t nvar=4;\r
- const Double_t ptBins[] = {0.05,0.1,0.12,0.14,0.16,0.18,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.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.0,3.2,3.4,3.6,3.8,4.0,4.2,4.4,4.6,4.8,5.0,5.5,6.0,6.5,7.0,7.5,8.0,8.5,9.0,9.5,10.,11.,12.,13.,14.,15.};\r
- Int_t nptBins=67;\r
- \r
- Int_t binsHistReal[nvar] = { 100, 100, nptBins, 20};\r
- Double_t xminHistReal[nvar] = { 0., 0, 0., 0.};\r
- Double_t xmaxHistReal[nvar] = { 10., 10., 5., 100.};\r
+ // binning common to all the THn\r
+ 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.};\r
+ const Int_t nptBins=26;\r
\r
- THnSparseF* NSparseHist = new THnSparseF("NSparseHist","NSparseHist",nvar,binsHistReal,xminHistReal,xmaxHistReal);\r
- NSparseHist->GetAxis(0)->SetTitle("Q vec VZERO-C");\r
- NSparseHist->GetAxis(1)->SetTitle("Q vec VZERO-A");\r
- NSparseHist->GetAxis(2)->SetTitle("#it{p}_{T}");\r
- NSparseHist->SetBinEdges(2,ptBins);\r
- NSparseHist->GetAxis(3)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));\r
- fOutput->Add(NSparseHist);\r
+ //dimensions of THnSparse for tracks\r
+ const Int_t nvartrk=8;\r
+ // pt cent Q vec IDrec IDgen isph iswd y\r
+ Int_t binsHistRealTrk[nvartrk] = { nptBins, fnCentBins, fnQvecBins, 3, 3, 2, 2, 2};\r
+ Double_t xminHistRealTrk[nvartrk] = { 0., 0., 0., -0.5, -0.5, -0.5, -0.5, -0.5};\r
+ Double_t xmaxHistRealTrk[nvartrk] = { 10., 100., 10., 2.5, 2.5, 1.5, 1.5, 0.5}; \r
+ THnSparseF* NSparseHistTrk = new THnSparseF("NSparseHistTrk","NSparseHistTrk",nvartrk,binsHistRealTrk,xminHistRealTrk,xmaxHistRealTrk);\r
+ NSparseHistTrk->GetAxis(0)->SetTitle("#it{p}_{T,rec}");\r
+ NSparseHistTrk->GetAxis(0)->SetName("pT_rec");\r
+ NSparseHistTrk->SetBinEdges(0,ptBins);\r
+ NSparseHistTrk->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));\r
+ NSparseHistTrk->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));\r
+ NSparseHistTrk->GetAxis(2)->SetTitle("Q vec");\r
+ NSparseHistTrk->GetAxis(2)->SetName("Q_vec");\r
+ NSparseHistTrk->GetAxis(3)->SetTitle("ID rec");\r
+ NSparseHistTrk->GetAxis(3)->SetName("ID_rec");\r
+ NSparseHistTrk->GetAxis(4)->SetTitle("ID gen");\r
+ NSparseHistTrk->GetAxis(4)->SetName("ID_gen");\r
+ NSparseHistTrk->GetAxis(5)->SetTitle("isph");\r
+ NSparseHistTrk->GetAxis(5)->SetName("isph");\r
+ NSparseHistTrk->GetAxis(6)->SetTitle("iswd");\r
+ NSparseHistTrk->GetAxis(6)->SetName("iswd");\r
+ NSparseHistTrk->GetAxis(7)->SetTitle("y");\r
+ NSparseHistTrk->GetAxis(7)->SetName("y");\r
+ fOutput->Add(NSparseHistTrk);\r
\r
- TH2F* hGen = new TH2F("hGen","hGen",nptBins,ptBins,20,0.,100.);\r
- hGen->GetXaxis()->SetTitle("#it{p}_{T,Gen}");\r
- hGen->GetYaxis()->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));\r
- fOutput->Add(hGen);\r
+ //dimensions of THnSparse for stack\r
+ const Int_t nvarst=5;\r
+ // pt cent IDgen isph y\r
+ Int_t binsHistRealSt[nvarst] = { nptBins, fnCentBins, 3, 2, 2};\r
+ Double_t xminHistRealSt[nvarst] = { 0., 0., -0.5, -0.5, -0.5};\r
+ Double_t xmaxHistRealSt[nvarst] = { 10., 100., 2.5, 1.5, 0.5};\r
+ THnSparseF* NSparseHistSt = new THnSparseF("NSparseHistSt","NSparseHistSt",nvarst,binsHistRealSt,xminHistRealSt,xmaxHistRealSt);\r
+ NSparseHistSt->GetAxis(0)->SetTitle("#it{p}_{T,gen}");\r
+ NSparseHistSt->SetBinEdges(0,ptBins);\r
+ NSparseHistSt->GetAxis(0)->SetName("pT_rec");\r
+ NSparseHistSt->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));\r
+ NSparseHistSt->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));\r
+ NSparseHistSt->GetAxis(2)->SetTitle("ID gen");\r
+ NSparseHistSt->GetAxis(2)->SetName("ID_gen");\r
+ NSparseHistSt->GetAxis(3)->SetTitle("isph");\r
+ NSparseHistSt->GetAxis(3)->SetName("isph");\r
+ NSparseHistSt->GetAxis(4)->SetTitle("y");\r
+ NSparseHistSt->GetAxis(4)->SetName("y");\r
+ fOutput->Add(NSparseHistSt);\r
\r
- TH2F* hRec = new TH2F("hRec","hRec",nptBins,ptBins,20,0.,100.);\r
- hRec->GetXaxis()->SetTitle("#it{p}_{T,Rec}");\r
- hRec->GetYaxis()->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));\r
- fOutput->Add(hRec);\r
+ //dimensions of THnSparse for the normalization\r
+ const Int_t nvarev=2;\r
+ // cent Q vec \r
+ Int_t binsHistRealEv[nvarev] = { fnCentBins, fnQvecBins};\r
+ Double_t xminHistRealEv[nvarev] = { 0., 0.};\r
+ Double_t xmaxHistRealEv[nvarev] = { 100., 10.};\r
+ THnSparseF* NSparseHistEv = new THnSparseF("NSparseHistEv","NSparseHistEv",nvarev,binsHistRealEv,xminHistRealEv,xmaxHistRealEv);\r
+ NSparseHistEv->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));\r
+ NSparseHistEv->GetAxis(0)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));\r
+ NSparseHistEv->GetAxis(1)->SetTitle("Q vec");\r
+ NSparseHistEv->GetAxis(1)->SetName("Q_vec");\r
+ fOutput->Add(NSparseHistEv);\r
\r
PostData(1, fOutput );\r
PostData(2, fEventCuts);\r
PostData(3, fTrackCuts);\r
- \r
+ PostData(4, fHelperPID);\r
}\r
+\r
//________________________________________________________________________\r
+\r
void AliAnalysisTaskSpectraAllChAOD::UserExec(Option_t *)\r
{\r
+ //Printf("An event");\r
// main event loop\r
fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);\r
if (!fAOD) {\r
}\r
\r
if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection\r
+\r
+ //Default TPC priors\r
+ if(fHelperPID->GetPIDType()==kBayes)fHelperPID->GetPIDCombined()->SetDefaultTPCPriors();//FIXME maybe this can go in the UserCreateOutputObject?\r
\r
- // First do MC to fill up the MC particle array, such that we can use it later\r
+ Double_t Qvec=0.;//in case of MC we save space in the memory\r
+ if(!fIsMC){\r
+ if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();\r
+ else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();\r
+ }\r
+\r
+ Double_t Cent=fEventCuts->GetCent();\r
+ \r
+ // First do MC to fill up the MC particle array\r
TClonesArray *arrayMC = 0;\r
if (fIsMC)\r
{\r
{\r
AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);\r
if(!partMC->Charge()) continue;//Skip neutrals\r
- if(!partMC->IsPhysicalPrimary()) continue; //skip secondaries\r
+ if(fCharge != 0 && partMC->Charge() != fCharge) continue;//if fCharge != 0 only select fCharge\r
+ \r
if(partMC->Eta() < fTrackCuts->GetEtaMin() || partMC->Eta() > fTrackCuts->GetEtaMax())continue;//ETA CUT ON GENERATED!!!!!!!!!!!!!!!!!!!!!!!!!!\r
\r
- ((TH2F*)fOutput->FindObject("hGen"))->Fill(partMC->Pt(),fEventCuts->GetCent());\r
+ //pt cent IDgen isph y\r
+ Double_t varSt[5];\r
+ varSt[0]=partMC->Pt();\r
+ varSt[1]=Cent;\r
+ varSt[2]=(Double_t)fHelperPID->GetParticleSpecies(partMC);\r
+ varSt[3]=(Double_t)partMC->IsPhysicalPrimary();\r
+ varSt[4]=partMC->Y();\r
+ ((THnSparseF*)fOutput->FindObject("NSparseHistSt"))->Fill(varSt);//stack loop\r
+ \r
//Printf("a particle");\r
- \r
+ \r
}\r
}\r
\r
//main loop on tracks\r
for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {\r
AliAODTrack* track = fAOD->GetTrack(iTracks);\r
- if (!fTrackCuts->IsSelected(track,kTRUE)) continue;\r
+ if(fCharge != 0 && track->Charge() != fCharge) continue;//if fCharge != 0 only select fCharge \r
+ if (!fTrackCuts->IsSelected(track,kTRUE)) continue; //track selection (rapidity selection NOT in the standard cuts)\r
+ Int_t IDrec=fHelperPID->GetParticleSpecies(track,kTRUE);//id from detector\r
+ Double_t y= track->Y(fHelperPID->GetMass((AliHelperParticleSpecies_t)IDrec));\r
+ Int_t IDgen=kSpUndefined;//set if MC\r
+ Int_t isph=-999;\r
+ Int_t iswd=-999;\r
\r
- ((TH2F*)fOutput->FindObject("hRec"))->Fill(track->Pt(),fEventCuts->GetCent());\r
+ if (arrayMC) {\r
+ AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));\r
+ if (!partMC) { \r
+ AliError("Cannot get MC particle");\r
+ continue; \r
+ }\r
+ IDgen=fHelperPID->GetParticleSpecies(partMC);\r
+ isph=partMC->IsPhysicalPrimary();\r
+ iswd=partMC->IsSecondaryFromWeakDecay();//FIXME not working on old productions\r
+ }\r
\r
- Double_t var[4];\r
- var[0]=fEventCuts->GetqV0C();\r
- var[1]=fEventCuts->GetqV0A();\r
- var[2]=track->Pt();\r
- var[3]=fEventCuts->GetCent();\r
- ((THnSparseF*)fOutput->FindObject("NSparseHist"))->Fill(var);\r
+ //pt cent Q vec IDrec IDgen isph iswd y\r
+ Double_t varTrk[8];\r
+ varTrk[0]=track->Pt();\r
+ varTrk[1]=Cent;\r
+ varTrk[2]=Qvec;\r
+ varTrk[3]=(Double_t)IDrec;\r
+ varTrk[4]=(Double_t)IDgen;\r
+ varTrk[5]=(Double_t)isph;\r
+ varTrk[6]=(Double_t)iswd;\r
+ varTrk[7]=y;\r
+ ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop\r
\r
//Printf("a track");\r
\r
} // end loop on tracks\r
\r
+ Double_t varEv[2];\r
+ varEv[0]=Cent;\r
+ varEv[1]=Qvec;\r
+ ((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop\r
+ \r
PostData(1, fOutput );\r
PostData(2, fEventCuts);\r
PostData(3, fTrackCuts);\r
+ PostData(4, fHelperPID);\r
}\r
\r
//_________________________________________________________________\r