PWGLF/SPECTRA/PiKaPr/TestAOD/AliAnalysisTaskSpectraAllChAOD.cxx

   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),
  60   fDoDoubleCounting(0),
  61   fCharge(0),
  62   fVZEROside(0),
  63   fOutput(0x0),
  64   fnCentBins(20),
  65   fnQvecBins(40)
  66 {
  67   // Default constructor
  68   DefineInput(0, TChain::Class());
  69   DefineOutput(1, TList::Class());
  70   DefineOutput(2, AliSpectraAODEventCuts::Class());
  71   DefineOutput(3, AliSpectraAODTrackCuts::Class());
  72   DefineOutput(4, AliHelperPID::Class());
  73 }
  74
  75 //________________________________________________________________________
  76 void AliAnalysisTaskSpectraAllChAOD::UserCreateOutputObjects()
  77 {
  78   // create output objects
  79   fOutput=new TList();
  80   fOutput->SetOwner();
  81   fOutput->SetName("fOutput");
  82
  83   if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");
  84   if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");
  85   if (!fHelperPID)  AliFatal("HelperPID should be set in the steering macro");
  86
  87   // binning common to all the THn
  88   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.};
  89   const Int_t nptBins=26;
  90
  91   //dimensions of THnSparse for tracks
  92   const Int_t nvartrk=8;
  93   //                                             pt          cent        Q vec     IDrec     IDgen       isph           iswd      y
  94   Int_t    binsHistRealTrk[nvartrk] = {      nptBins, fnCentBins,   fnQvecBins,        4,        3,         2,          2,       2};
  95   Double_t xminHistRealTrk[nvartrk] = {         0.,          0.,            0.,      -.5,      -0.5,      -0.5,        -0.5,   -0.5};
  96   Double_t xmaxHistRealTrk[nvartrk] = {       10.,       100.,            8.,      3.5,      2.5,       1.5,         1.5,     0.5};
  97   THnSparseF* NSparseHistTrk = new THnSparseF("NSparseHistTrk","NSparseHistTrk",nvartrk,binsHistRealTrk,xminHistRealTrk,xmaxHistRealTrk);
  98   NSparseHistTrk->GetAxis(0)->SetTitle("#it{p}_{T,rec}");
  99   NSparseHistTrk->GetAxis(0)->SetName("pT_rec");
 100   NSparseHistTrk->SetBinEdges(0,ptBins);
 101   NSparseHistTrk->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
 102   NSparseHistTrk->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
 103   NSparseHistTrk->GetAxis(2)->SetTitle("Q vec");
 104   NSparseHistTrk->GetAxis(2)->SetName("Q_vec");
 105   NSparseHistTrk->GetAxis(3)->SetTitle("ID rec");
 106   NSparseHistTrk->GetAxis(3)->SetName("ID_rec");
 107   NSparseHistTrk->GetAxis(4)->SetTitle("ID gen");
 108   NSparseHistTrk->GetAxis(4)->SetName("ID_gen");
 109   NSparseHistTrk->GetAxis(5)->SetTitle("isph");
 110   NSparseHistTrk->GetAxis(5)->SetName("isph");
 111   NSparseHistTrk->GetAxis(6)->SetTitle("iswd");
 112   NSparseHistTrk->GetAxis(6)->SetName("iswd");
 113   NSparseHistTrk->GetAxis(7)->SetTitle("y");
 114   NSparseHistTrk->GetAxis(7)->SetName("y");
 115   fOutput->Add(NSparseHistTrk);
 116
 117   //dimensions of THnSparse for stack
 118   const Int_t nvarst=5;
 119   //                                             pt          cent    IDgen        isph        y
 120   Int_t    binsHistRealSt[nvarst] = {      nptBins,  fnCentBins,        3,         2,        2};
 121   Double_t xminHistRealSt[nvarst] = {         0.,           0.,      -0.5,      -0.5,    -0.5};
 122   Double_t xmaxHistRealSt[nvarst] = {       10.,        100.,      2.5,       1.5,      0.5};
 123   THnSparseF* NSparseHistSt = new THnSparseF("NSparseHistSt","NSparseHistSt",nvarst,binsHistRealSt,xminHistRealSt,xmaxHistRealSt);
 124   NSparseHistSt->GetAxis(0)->SetTitle("#it{p}_{T,gen}");
 125   NSparseHistSt->SetBinEdges(0,ptBins);
 126   NSparseHistSt->GetAxis(0)->SetName("pT_rec");
 127   NSparseHistSt->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
 128   NSparseHistSt->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
 129   NSparseHistSt->GetAxis(2)->SetTitle("ID gen");
 130   NSparseHistSt->GetAxis(2)->SetName("ID_gen");
 131   NSparseHistSt->GetAxis(3)->SetTitle("isph");
 132   NSparseHistSt->GetAxis(3)->SetName("isph");
 133   NSparseHistSt->GetAxis(4)->SetTitle("y");
 134   NSparseHistSt->GetAxis(4)->SetName("y");
 135   fOutput->Add(NSparseHistSt);
 136
 137   //dimensions of THnSparse for the normalization
 138   const Int_t nvarev=3;
 139   //                                             cent             Q vec          Nch
 140   Int_t    binsHistRealEv[nvarev] = {    fnCentBins,      fnQvecBins,           360.};
 141   Double_t xminHistRealEv[nvarev] = {           0.,               0.,             0.};
 142   Double_t xmaxHistRealEv[nvarev] = {       100.,               8.,            1800.};
 143   THnSparseF* NSparseHistEv = new THnSparseF("NSparseHistEv","NSparseHistEv",nvarev,binsHistRealEv,xminHistRealEv,xmaxHistRealEv);
 144   NSparseHistEv->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
 145   NSparseHistEv->GetAxis(0)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
 146   NSparseHistEv->GetAxis(1)->SetTitle("Q vec");
 147   NSparseHistEv->GetAxis(1)->SetName("Q_vec");
 148   NSparseHistEv->GetAxis(2)->SetTitle("N charged");
 149   NSparseHistEv->GetAxis(2)->SetName("N_ch");
 150   fOutput->Add(NSparseHistEv);
 151
 152   PostData(1, fOutput  );
 153   PostData(2, fEventCuts);
 154   PostData(3, fTrackCuts);
 155   PostData(4, fHelperPID);
 156 }
 157
 158 //________________________________________________________________________
 159
 160 void AliAnalysisTaskSpectraAllChAOD::UserExec(Option_t *)
 161 {
 162   //Printf("An event");
 163   // main event loop
 164   fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);
 165   if (!fAOD) {
 166     AliWarning("ERROR: AliAODEvent not available \n");
 167     return;
 168   }
 169
 170   if (strcmp(fAOD->ClassName(), "AliAODEvent"))
 171     {
 172       AliFatal("Not processing AODs");
 173     }
 174
 175   if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection
 176
 177   //Default TPC priors
 178   if(fHelperPID->GetPIDType()==kBayes)fHelperPID->GetPIDCombined()->SetDefaultTPCPriors();//FIXME maybe this can go in the UserCreateOutputObject?
 179
 180   Double_t Qvec=0.;//in case of MC we save space in the memory
 181   if(!fIsMC){
 182     if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();
 183     else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();
 184   }
 185
 186   Double_t Cent=fEventCuts->GetCent();
 187
 188   // First do MC to fill up the MC particle array
 189   TClonesArray *arrayMC = 0;
 190   if (fIsMC)
 191     {
 192       arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 193       if (!arrayMC) {
 194         AliFatal("Error: MC particles branch not found!\n");
 195       }
 196       Int_t nMC = arrayMC->GetEntries();
 197       for (Int_t iMC = 0; iMC < nMC; iMC++)
 198         {
 199           AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
 200           if(!partMC->Charge()) continue;//Skip neutrals
 201           if(fCharge != 0 && partMC->Charge()*fCharge < 0.) continue;//if fCharge != 0 only select fCharge
 202
 203           if(partMC->Eta() < fTrackCuts->GetEtaMin() || partMC->Eta() > fTrackCuts->GetEtaMax())continue;//ETA CUT ON GENERATED!!!!!!!!!!!!!!!!!!!!!!!!!!
 204
 205           //pt     cent    IDgen        isph        y
 206           Double_t varSt[5];
 207           varSt[0]=partMC->Pt();
 208           varSt[1]=Cent;
 209           varSt[2]=(Double_t)fHelperPID->GetParticleSpecies(partMC);
 210           varSt[3]=(Double_t)partMC->IsPhysicalPrimary();
 211           varSt[4]=partMC->Y();
 212           ((THnSparseF*)fOutput->FindObject("NSparseHistSt"))->Fill(varSt);//stack loop
 213
 214           //Printf("a particle");
 215
 216         }
 217     }
 218
 219   //main loop on tracks
 220
 221   Int_t Nch = 0.;
 222
 223   for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
 224     AliAODTrack* track = fAOD->GetTrack(iTracks);
 225     if(fCharge != 0 && track->Charge() != fCharge) continue;//if fCharge != 0 only select fCharge
 226     if (!fTrackCuts->IsSelected(track,kTRUE)) continue; //track selection (rapidity selection NOT in the standard cuts)
 227     Int_t IDrec=fHelperPID->GetParticleSpecies(track,kTRUE);//id from detector
 228     Double_t y= track->Y(fHelperPID->GetMass((AliHelperParticleSpecies_t)IDrec));
 229     Int_t IDgen=kSpUndefined;//set if MC
 230     Int_t isph=-999;
 231     Int_t iswd=-999;
 232
 233     if (arrayMC) {
 234       AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));
 235       if (!partMC) {
 236         AliError("Cannot get MC particle");
 237         continue;
 238       }
 239       IDgen=fHelperPID->GetParticleSpecies(partMC);
 240       isph=partMC->IsPhysicalPrimary();
 241       iswd=partMC->IsSecondaryFromWeakDecay();//FIXME not working on old productions
 242     }
 243
 244     //pt     cent    Q vec     IDrec     IDgen       isph           iswd      y
 245     Double_t varTrk[8];
 246     varTrk[0]=track->Pt();
 247     varTrk[1]=Cent;
 248     varTrk[2]=Qvec;
 249     varTrk[3]=(Double_t)IDrec;
 250     varTrk[4]=(Double_t)IDgen;
 251     varTrk[5]=(Double_t)isph;
 252     varTrk[6]=(Double_t)iswd;
 253     varTrk[7]=y;
 254     ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
 255
 256     //for nsigma PID fill double counting of ID
 257     if(fHelperPID->GetPIDType()<kBayes && fDoDoubleCounting){//only nsigma
 258       Bool_t *HasDC;
 259       HasDC=fHelperPID->GetDoubleCounting(track,kTRUE);//get the array with double counting
 260       for(Int_t ipart=0;ipart<kNSpecies;ipart++){
 261         if(HasDC[ipart]==kTRUE){
 262           varTrk[3]=(Double_t)ipart;
 263           ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
 264         }
 265       }
 266     }
 267
 268     //fill all charged (3)
 269     varTrk[3]=3.;
 270     varTrk[4]=3.;
 271     Nch++;
 272     ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
 273
 274     //Printf("a track");
 275
 276   } // end loop on tracks
 277
 278   Double_t varEv[3];
 279   varEv[0]=Cent;
 280   varEv[1]=Qvec;
 281   varEv[2]=Nch;
 282   ((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
 283
 284   PostData(1, fOutput  );
 285   PostData(2, fEventCuts);
 286   PostData(3, fTrackCuts);
 287   PostData(4, fHelperPID);
 288 }
 289
 290 //_________________________________________________________________
 291 void   AliAnalysisTaskSpectraAllChAOD::Terminate(Option_t *)
 292 {
 293   // Terminate
 294 }