PWG/DevNanoAOD/AliAnalysisTaskSpectraAllChNanoAOD.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 //         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 }