PWGLF/SPECTRA/PiKaPr/TestAOD/AliAnalysisTaskSpectraAllChAOD.cxx

   1 /**************************************************************************\r
   2  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
   3  *                                                                        *\r
   4  * Author: The ALICE Off-line Project.                                    *\r
   5  * Contributors are mentioned in the code where appropriate.              *\r
   6  *                                                                        *\r
   7  * Permission to use, copy, modify and distribute this software and its   *\r
   8  * documentation strictly for non-commercial purposes is hereby granted   *\r
   9  * without fee, provided that the above copyright notice appears in all   *\r
  10  * copies and that both the copyright notice and this permission notice   *\r
  11  * appear in the supporting documentation. The authors make no claims     *\r
  12  * about the suitability of this software for any purpose. It is          *\r
  13  * provided "as is" without express or implied warranty.                  *\r
  14  **************************************************************************/\r
  15 \r
  16 //-----------------------------------------------------------------\r
  17 //         AliAnalysisTaskSpectraAllChAOD class\r
  18 //-----------------------------------------------------------------\r
  19 \r
  20 #include "TChain.h"\r
  21 #include "TTree.h"\r
  22 #include "TLegend.h"\r
  23 #include "TH1F.h"\r
  24 #include "TH2F.h"\r
  25 #include "THnSparse.h"\r
  26 #include "TCanvas.h"\r
  27 #include "AliAnalysisTask.h"\r
  28 #include "AliAODTrack.h"\r
  29 #include "AliAODMCParticle.h"\r
  30 #include "AliVParticle.h"\r
  31 #include "AliAODEvent.h"\r
  32 #include "AliAODInputHandler.h"\r
  33 #include "AliAnalysisTaskSpectraAllChAOD.h"\r
  34 #include "AliAnalysisTaskESDfilter.h"\r
  35 #include "AliAnalysisDataContainer.h"\r
  36 #include "AliSpectraAODTrackCuts.h"\r
  37 #include "AliSpectraAODEventCuts.h"\r
  38 #include "AliHelperPID.h"\r
  39 #include "AliCentrality.h"\r
  40 #include "TProof.h"\r
  41 #include "AliVEvent.h"\r
  42 #include "AliStack.h"\r
  43 #include <TMCProcess.h>\r
  44 \r
  45 #include <iostream>\r
  46 \r
  47 using namespace AliHelperPIDNameSpace;\r
  48 using namespace std;\r
  49 \r
  50 ClassImp(AliAnalysisTaskSpectraAllChAOD)\r
  51 \r
  52 //________________________________________________________________________\r
  53 AliAnalysisTaskSpectraAllChAOD::AliAnalysisTaskSpectraAllChAOD(const char *name) : AliAnalysisTaskSE(name),  \r
  54   fAOD(0x0),\r
  55   fTrackCuts(0x0),\r
  56   fEventCuts(0x0),\r
  57   fHelperPID(0x0),\r
  58   fIsMC(0),\r
  59   fOutput(0x0),\r
  60   fnCentBins(20),\r
  61   fnQvecBins(50)\r
  62 {\r
  63   // Default constructor\r
  64   DefineInput(0, TChain::Class());\r
  65   DefineOutput(1, TList::Class());\r
  66   DefineOutput(2, AliSpectraAODEventCuts::Class());\r
  67   DefineOutput(3, AliSpectraAODTrackCuts::Class());\r
  68   DefineOutput(4, AliHelperPID::Class());\r
  69 }\r
  70 \r
  71 //________________________________________________________________________\r
  72 void AliAnalysisTaskSpectraAllChAOD::UserCreateOutputObjects()\r
  73 {\r
  74   // create output objects\r
  75   fOutput=new TList();\r
  76   fOutput->SetOwner();\r
  77   fOutput->SetName("fOutput");\r
  78   \r
  79   if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");\r
  80   if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");\r
  81   if (!fHelperPID)  AliFatal("HelperPID should be set in the steering macro");\r
  82   \r
  83   // binning common to all the THn\r
  84   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
  85   const Int_t nptBins=26;\r
  86   \r
  87   //dimensions of THnSparse for tracks\r
  88   const Int_t nvartrk=8;\r
  89   //                                             pt          cent        Q vec     IDrec     IDgen       isph           iswd      y\r
  90   Int_t    binsHistRealTrk[nvartrk] = {      nptBins, fnCentBins,   fnQvecBins,        3,        3,         2,          2,       2};\r
  91   Double_t xminHistRealTrk[nvartrk] = {         0.,          0.,            0.,     -0.5,      -0.5,      -0.5,        -0.5,   -0.5};\r
  92   Double_t xmaxHistRealTrk[nvartrk] = {       10.,       100.,          10.,      2.5,      2.5,       1.5,         1.5,     0.5};    \r
  93   THnSparseF* NSparseHistTrk = new THnSparseF("NSparseHistTrk","NSparseHistTrk",nvartrk,binsHistRealTrk,xminHistRealTrk,xmaxHistRealTrk);\r
  94   NSparseHistTrk->GetAxis(0)->SetTitle("#it{p}_{T,rec}");\r
  95   NSparseHistTrk->GetAxis(0)->SetName("pT_rec");\r
  96   NSparseHistTrk->SetBinEdges(0,ptBins);\r
  97   NSparseHistTrk->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));\r
  98   NSparseHistTrk->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));\r
  99   NSparseHistTrk->GetAxis(2)->SetTitle("Q vec");\r
 100   NSparseHistTrk->GetAxis(2)->SetName("Q_vec");\r
 101   NSparseHistTrk->GetAxis(3)->SetTitle("ID rec");\r
 102   NSparseHistTrk->GetAxis(3)->SetName("ID_rec");\r
 103   NSparseHistTrk->GetAxis(4)->SetTitle("ID gen");\r
 104   NSparseHistTrk->GetAxis(4)->SetName("ID_gen");\r
 105   NSparseHistTrk->GetAxis(5)->SetTitle("isph");\r
 106   NSparseHistTrk->GetAxis(5)->SetName("isph");\r
 107   NSparseHistTrk->GetAxis(6)->SetTitle("iswd");\r
 108   NSparseHistTrk->GetAxis(6)->SetName("iswd");\r
 109   NSparseHistTrk->GetAxis(7)->SetTitle("y");\r
 110   NSparseHistTrk->GetAxis(7)->SetName("y");\r
 111   fOutput->Add(NSparseHistTrk);\r
 112   \r
 113   //dimensions of THnSparse for stack\r
 114   const Int_t nvarst=5;\r
 115   //                                             pt          cent    IDgen        isph        y\r
 116   Int_t    binsHistRealSt[nvarst] = {      nptBins,  fnCentBins,        3,         2,        2};\r
 117   Double_t xminHistRealSt[nvarst] = {         0.,           0.,      -0.5,      -0.5,    -0.5};\r
 118   Double_t xmaxHistRealSt[nvarst] = {       10.,        100.,      2.5,       1.5,      0.5};\r
 119   THnSparseF* NSparseHistSt = new THnSparseF("NSparseHistSt","NSparseHistSt",nvarst,binsHistRealSt,xminHistRealSt,xmaxHistRealSt);\r
 120   NSparseHistSt->GetAxis(0)->SetTitle("#it{p}_{T,gen}");\r
 121   NSparseHistSt->SetBinEdges(0,ptBins);\r
 122   NSparseHistSt->GetAxis(0)->SetName("pT_rec");\r
 123   NSparseHistSt->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));\r
 124   NSparseHistSt->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));\r
 125   NSparseHistSt->GetAxis(2)->SetTitle("ID gen");\r
 126   NSparseHistSt->GetAxis(2)->SetName("ID_gen");\r
 127   NSparseHistSt->GetAxis(3)->SetTitle("isph");\r
 128   NSparseHistSt->GetAxis(3)->SetName("isph");\r
 129   NSparseHistSt->GetAxis(4)->SetTitle("y");\r
 130   NSparseHistSt->GetAxis(4)->SetName("y");\r
 131   fOutput->Add(NSparseHistSt);\r
 132   \r
 133   //dimensions of THnSparse for the normalization\r
 134   const Int_t nvarev=2;\r
 135   //                                             cent             Q vec   \r
 136   Int_t    binsHistRealEv[nvarev] = {    fnCentBins,      fnQvecBins};\r
 137   Double_t xminHistRealEv[nvarev] = {           0.,               0.};\r
 138   Double_t xmaxHistRealEv[nvarev] = {       100.,              10.};\r
 139   THnSparseF* NSparseHistEv = new THnSparseF("NSparseHistEv","NSparseHistEv",nvarev,binsHistRealEv,xminHistRealEv,xmaxHistRealEv);\r
 140   NSparseHistEv->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));\r
 141   NSparseHistEv->GetAxis(0)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));\r
 142   NSparseHistEv->GetAxis(1)->SetTitle("Q vec");\r
 143   NSparseHistEv->GetAxis(1)->SetName("Q_vec");\r
 144   fOutput->Add(NSparseHistEv);\r
 145   \r
 146   PostData(1, fOutput  );\r
 147   PostData(2, fEventCuts);\r
 148   PostData(3, fTrackCuts);\r
 149   PostData(4, fHelperPID);\r
 150 }\r
 151 \r
 152 //________________________________________________________________________\r
 153 \r
 154 void AliAnalysisTaskSpectraAllChAOD::UserExec(Option_t *)\r
 155 {\r
 156   //Printf("An event");\r
 157   // main event loop\r
 158   fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);\r
 159   if (!fAOD) {\r
 160     AliWarning("ERROR: AliAODEvent not available \n");\r
 161     return;\r
 162   }\r
 163   \r
 164   if (strcmp(fAOD->ClassName(), "AliAODEvent"))\r
 165     {\r
 166       AliFatal("Not processing AODs");\r
 167     }\r
 168   \r
 169   if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection\r
 170   \r
 171   Double_t Qvec=0.;//in case of MC we save space in the memory\r
 172   if(!fIsMC)Qvec=fEventCuts->GetqV0A();//FIXME we'll have to combine A and C\r
 173   Double_t Cent=fEventCuts->GetCent();\r
 174     \r
 175   // First do MC to fill up the MC particle array\r
 176   TClonesArray *arrayMC = 0;\r
 177   if (fIsMC)\r
 178     {\r
 179       arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());\r
 180       if (!arrayMC) {\r
 181         AliFatal("Error: MC particles branch not found!\n");\r
 182       }\r
 183       Int_t nMC = arrayMC->GetEntries();\r
 184       for (Int_t iMC = 0; iMC < nMC; iMC++)\r
 185         {\r
 186           AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);\r
 187           if(!partMC->Charge()) continue;//Skip neutrals\r
 188           if(partMC->Eta() < fTrackCuts->GetEtaMin() || partMC->Eta() > fTrackCuts->GetEtaMax())continue;//ETA CUT ON GENERATED!!!!!!!!!!!!!!!!!!!!!!!!!!\r
 189           \r
 190           //pt     cent    IDgen        isph        y\r
 191           Double_t varSt[5];\r
 192           varSt[0]=partMC->Pt();\r
 193           varSt[1]=Cent;\r
 194           varSt[2]=(Double_t)fHelperPID->GetParticleSpecies(partMC);\r
 195           varSt[3]=(Double_t)partMC->IsPhysicalPrimary();\r
 196           varSt[4]=partMC->Y();\r
 197           ((THnSparseF*)fOutput->FindObject("NSparseHistSt"))->Fill(varSt);//stack loop\r
 198           \r
 199           //Printf("a particle");\r
 200           \r
 201         }\r
 202     }\r
 203   \r
 204   //main loop on tracks\r
 205   for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {\r
 206     AliAODTrack* track = fAOD->GetTrack(iTracks);\r
 207     if (!fTrackCuts->IsSelected(track,kTRUE)) continue; //track selection (rapidity selection NOT in the standard cuts)\r
 208     Int_t IDrec=fHelperPID->GetParticleSpecies(track,kTRUE);//id from detector\r
 209     Double_t y= track->Y(fHelperPID->GetMass((AliHelperParticleSpecies_t)IDrec));\r
 210     Int_t IDgen=kSpUndefined;//set if MC\r
 211     Int_t isph=-999;\r
 212     Int_t iswd=-999;\r
 213     \r
 214     if (arrayMC) {\r
 215       AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));\r
 216       if (!partMC) { \r
 217         AliError("Cannot get MC particle");\r
 218         continue; \r
 219       }\r
 220       IDgen=fHelperPID->GetParticleSpecies(partMC);\r
 221       isph=partMC->IsPhysicalPrimary();\r
 222       iswd=partMC->IsSecondaryFromWeakDecay();//FIXME not working on old productions\r
 223     }\r
 224     \r
 225     //pt     cent    Q vec     IDrec     IDgen       isph           iswd      y\r
 226     Double_t varTrk[8];\r
 227     varTrk[0]=track->Pt();\r
 228     varTrk[1]=Cent;\r
 229     varTrk[2]=Qvec;\r
 230     varTrk[3]=(Double_t)IDrec;\r
 231     varTrk[4]=(Double_t)IDgen;\r
 232     varTrk[5]=(Double_t)isph;\r
 233     varTrk[6]=(Double_t)iswd;\r
 234     varTrk[7]=y;\r
 235     ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop\r
 236     \r
 237     //Printf("a track");\r
 238     \r
 239   } // end loop on tracks\r
 240   \r
 241   Double_t varEv[2];\r
 242   varEv[0]=Cent;\r
 243   varEv[1]=Qvec;\r
 244   ((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop\r
 245   \r
 246   PostData(1, fOutput  );\r
 247   PostData(2, fEventCuts);\r
 248   PostData(3, fTrackCuts);\r
 249   PostData(4, fHelperPID);\r
 250 }\r
 251 \r
 252 //_________________________________________________________________\r
 253 void   AliAnalysisTaskSpectraAllChAOD::Terminate(Option_t *)\r
 254 {\r
 255   // Terminate\r
 256 }\r