]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGLF/SPECTRA/PiKaPr/TestAOD/AliAnalysisTaskSpectraAllChAOD.cxx
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / AliAnalysisTaskSpectraAllChAOD.cxx
index 82f24045f60862e8bb6c7a36b93db8e9268e583c..f8898b2acf0abea997fa22731ce62e3b1db996e1 100644 (file)
@@ -35,6 +35,8 @@
 #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
@@ -69,42 +81,82 @@ void AliAnalysisTaskSpectraAllChAOD::UserCreateOutputObjects()
   \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
@@ -118,8 +170,19 @@ void AliAnalysisTaskSpectraAllChAOD::UserExec(Option_t *)
     }\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
@@ -132,36 +195,71 @@ void AliAnalysisTaskSpectraAllChAOD::UserExec(Option_t *)
        {\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