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 fFillOnlyEvents(0),
  62 fCharge(0),
  63 fVZEROside(0),
  64 fOutput(0x0),
  65 fnCentBins(20),
  66 fnQvecBins(100),
  67 fnNchBins(200),
  68 fIsQvecCalibMode(0),
  69 fQvecUpperLim(100),
  70 fIsAOD160(1),
  71 fnDCABins(60),
  72 fDCAmin(-3),
  73 fDCAmax(3),
  74 fDCAzCut(999999.),
  75 fQst(1),
  76 fQtrk(0),
  77 fQgenType(0),
  78 fDoCentrSystCentrality(0)
  79 {
  80   // Default constructor
  81   DefineInput(0, TChain::Class());
  82   DefineOutput(1, TList::Class());
  83   DefineOutput(2, AliSpectraAODEventCuts::Class());
  84   DefineOutput(3, AliSpectraAODTrackCuts::Class());
  85   DefineOutput(4, AliHelperPID::Class());
  86 }
  87
  88 //________________________________________________________________________
  89 void AliAnalysisTaskSpectraAllChAOD::UserCreateOutputObjects()
  90 {
  91   // create output objects
  92   fOutput=new TList();
  93   fOutput->SetOwner();
  94   fOutput->SetName("fOutput");
  95
  96   if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");
  97   if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");
  98   if (!fHelperPID)  AliFatal("HelperPID should be set in the steering macro");
  99
 100   // binning common to all the THn
 101   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.,20.,25.,30.,35.,40.,50.,75.,100.};
 102   const Int_t nptBins=34;
 103   //const Double_t ptBins[] = {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.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.,20.,25.,30.,35.,40.,50.,75.,100.};
 104   //const Int_t nptBins=44;
 105
 106   //dimensions of THnSparse for tracks
 107   //const Int_t nvartrk=8;
 108   //                                             pt          cent          Q vec       IDrec          IDgen        isph         y             DCA
 109   //  Int_t    binsHistRealTrk[nvartrk] = {      nptBins, fnCentBins,     fnQvecBins,        4,             3,          3,        2,        fnDCABins};
 110   //  Double_t xminHistRealTrk[nvartrk] = {         0.,          0.,              0.,       -.5,          -0.5,        0.5,      -0.5,          fDCAmin};
 111   //  Double_t xmaxHistRealTrk[nvartrk] = {       10.,       100.,     fQvecUpperLim,       3.5,      2.5,        3.5,      0.5,         fDCAmax};
 112   //  THnSparseF* NSparseHistTrk = new THnSparseF("NSparseHistTrk","NSparseHistTrk",nvartrk,binsHistRealTrk,xminHistRealTrk,xmaxHistRealTrk);
 113   //  NSparseHistTrk->GetAxis(0)->SetTitle("#it{p}_{T,rec}");
 114   //  NSparseHistTrk->GetAxis(0)->SetName("pT_rec");
 115   //  NSparseHistTrk->SetBinEdges(0,ptBins);
 116   //  NSparseHistTrk->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
 117   //  NSparseHistTrk->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
 118   //  NSparseHistTrk->GetAxis(2)->SetTitle("Q vec");
 119   //  NSparseHistTrk->GetAxis(2)->SetName("Q_vec");
 120   //  NSparseHistTrk->GetAxis(3)->SetTitle("ID rec");
 121   //  NSparseHistTrk->GetAxis(3)->SetName("ID_rec");
 122   //  NSparseHistTrk->GetAxis(4)->SetTitle("ID gen");
 123   //  NSparseHistTrk->GetAxis(4)->SetName("ID_gen");
 124   //  NSparseHistTrk->GetAxis(5)->SetTitle("isph");
 125   //  NSparseHistTrk->GetAxis(5)->SetName("isph");
 126   //  NSparseHistTrk->GetAxis(6)->SetTitle("y");
 127   //  NSparseHistTrk->GetAxis(6)->SetName("y");
 128   //  NSparseHistTrk->GetAxis(7)->SetTitle("dca");
 129   //  NSparseHistTrk->GetAxis(7)->SetName("dca");
 130   const Int_t nvartrk=7;
 131   Int_t    binsHistRealTrk[nvartrk] = {      nptBins, fnCentBins,     fnQvecBins,        4,             3,          3,        2};
 132   Double_t xminHistRealTrk[nvartrk] = {         0.,          0.,              0.,       -.5,          -0.5,        0.5,      -0.5};
 133   Double_t xmaxHistRealTrk[nvartrk] = {       10.,       100.,     fQvecUpperLim,       3.5,      2.5,        3.5,      0.5};
 134   THnSparseF* NSparseHistTrk = new THnSparseF("NSparseHistTrk","NSparseHistTrk",nvartrk,binsHistRealTrk,xminHistRealTrk,xmaxHistRealTrk);
 135   NSparseHistTrk->GetAxis(0)->SetTitle("#it{p}_{T,rec}");
 136   NSparseHistTrk->GetAxis(0)->SetName("pT_rec");
 137   NSparseHistTrk->SetBinEdges(0,ptBins);
 138   NSparseHistTrk->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
 139   NSparseHistTrk->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
 140   NSparseHistTrk->GetAxis(2)->SetTitle("Q vec");
 141   NSparseHistTrk->GetAxis(2)->SetName("Q_vec");
 142   NSparseHistTrk->GetAxis(3)->SetTitle("ID rec");
 143   NSparseHistTrk->GetAxis(3)->SetName("ID_rec");
 144   NSparseHistTrk->GetAxis(4)->SetTitle("ID gen");
 145   NSparseHistTrk->GetAxis(4)->SetName("ID_gen");
 146   NSparseHistTrk->GetAxis(5)->SetTitle("isph");
 147   NSparseHistTrk->GetAxis(5)->SetName("isph");
 148   NSparseHistTrk->GetAxis(6)->SetTitle("y");
 149   NSparseHistTrk->GetAxis(6)->SetName("y");
 150
 151
 152   fOutput->Add(NSparseHistTrk);
 153
 154   //dimensions of THnSparse for stack
 155   const Int_t nvarst=7;
 156   //                                             pt          cent    IDgen        isph        y    etaselected       Qvec gen
 157   Int_t    binsHistRealSt[nvarst] = {      nptBins,  fnCentBins,        3,         2,        2,             1,      fnQvecBins};
 158   Double_t xminHistRealSt[nvarst] = {         0.,           0.,      -0.5,      -0.5,    -0.5,            0.5,              0.};
 159   Double_t xmaxHistRealSt[nvarst] = {       10.,        100.,      2.5,       1.5,      0.5,           1.5,      fQvecUpperLim};
 160   THnSparseF* NSparseHistSt = new THnSparseF("NSparseHistSt","NSparseHistSt",nvarst,binsHistRealSt,xminHistRealSt,xmaxHistRealSt);
 161   NSparseHistSt->GetAxis(0)->SetTitle("#it{p}_{T,gen}");
 162   NSparseHistSt->SetBinEdges(0,ptBins);
 163   NSparseHistSt->GetAxis(0)->SetName("pT_rec");
 164   NSparseHistSt->GetAxis(1)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
 165   NSparseHistSt->GetAxis(1)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
 166   NSparseHistSt->GetAxis(2)->SetTitle("ID gen");
 167   NSparseHistSt->GetAxis(2)->SetName("ID_gen");
 168   NSparseHistSt->GetAxis(3)->SetTitle("isph");
 169   NSparseHistSt->GetAxis(3)->SetName("isph");
 170   NSparseHistSt->GetAxis(4)->SetTitle("y");
 171   NSparseHistSt->GetAxis(4)->SetName("y");
 172   NSparseHistSt->GetAxis(5)->SetTitle("etaselected");
 173   NSparseHistSt->GetAxis(5)->SetName("etaselected");
 174   NSparseHistSt->GetAxis(6)->SetTitle("Q vec gen");
 175   NSparseHistSt->GetAxis(6)->SetName("Q_vec_gen");
 176   fOutput->Add(NSparseHistSt);
 177
 178   //dimensions of THnSparse for the normalization
 179   const Int_t nvarev=3;
 180   //                                             cent             Q vec                Nch
 181   Int_t    binsHistRealEv[nvarev] = {    fnCentBins,      fnQvecBins,           fnNchBins};
 182   Double_t xminHistRealEv[nvarev] = {           0.,               0.,                   0.};
 183   Double_t xmaxHistRealEv[nvarev] = {       100.,  fQvecUpperLim,               2000.};
 184   THnSparseF* NSparseHistEv = new THnSparseF("NSparseHistEv","NSparseHistEv",nvarev,binsHistRealEv,xminHistRealEv,xmaxHistRealEv);
 185   NSparseHistEv->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
 186   NSparseHistEv->GetAxis(0)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
 187   NSparseHistEv->GetAxis(1)->SetTitle("Q vec");
 188   NSparseHistEv->GetAxis(1)->SetName("Q_vec");
 189   NSparseHistEv->GetAxis(2)->SetTitle("N charged");
 190   NSparseHistEv->GetAxis(2)->SetName("N_ch");
 191   fOutput->Add(NSparseHistEv);
 192
 193   PostData(1, fOutput  );
 194   PostData(2, fEventCuts);
 195   PostData(3, fTrackCuts);
 196   PostData(4, fHelperPID);
 197 }
 198
 199 //________________________________________________________________________
 200
 201 void AliAnalysisTaskSpectraAllChAOD::UserExec(Option_t *)
 202 {
 203   //Printf("An event");
 204   // main event loop
 205   fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);
 206   if (!fAOD) {
 207     AliWarning("ERROR: AliAODEvent not available \n");
 208     return;
 209   }
 210
 211   if (strcmp(fAOD->ClassName(), "AliAODEvent"))
 212   {
 213     AliFatal("Not processing AODs");
 214   }
 215
 216   if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection
 217
 218   //Default TPC priors
 219   if(fHelperPID->GetPIDType()==kBayes)fHelperPID->GetPIDCombined()->SetDefaultTPCPriors();//FIXME we should modify the task to change priors
 220
 221   Double_t Qvec=0.;
 222   if(fIsQvecCalibMode){
 223     if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();
 224     else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();
 225     else if (fVZEROside==2)Qvec=fEventCuts->GetqTPC();
 226   }
 227   else Qvec=fEventCuts->GetQvecPercentile(fVZEROside);
 228
 229   Double_t QvecMC = 0.;
 230   if(fIsMC){
 231     if(fIsQvecCalibMode){
 232       QvecMC = fEventCuts->CalculateQVectorMC(fVZEROside, fQgenType);
 233     }
 234     else QvecMC = fEventCuts->GetQvecPercentileMC(fVZEROside, fQgenType);
 235   }
 236
 237   Double_t Cent=(fDoCentrSystCentrality)?1.01*fEventCuts->GetCent():fEventCuts->GetCent();
 238
 239   // First do MC to fill up the MC particle array
 240   TClonesArray *arrayMC = 0;
 241   if (fIsMC)
 242   {
 243     arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 244     if (!arrayMC) {
 245       AliFatal("Error: MC particles branch not found!\n");
 246     }
 247     Int_t nMC = arrayMC->GetEntries();
 248     for (Int_t iMC = 0; iMC < nMC; iMC++)
 249     {
 250       AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
 251       if(!partMC->Charge()) continue;//Skip neutrals
 252       if(fCharge != 0 && partMC->Charge()*fCharge < 0.) continue;//if fCharge != 0 only select fCharge
 253
 254       //flag to select particle in the same ETA acceptance of the tracks (to be used for the comparison with AllCh)
 255       Double_t etaselected=-1.;
 256       if(partMC->Eta()>=fTrackCuts->GetEtaMin() && partMC->Eta()<=fTrackCuts->GetEtaMax())etaselected=1.;
 257
 258       //pt     cent    IDgen        isph        y
 259       Double_t varSt[7];
 260       varSt[0]=partMC->Pt();
 261       varSt[1]=Cent;
 262       varSt[2]=(Double_t)fHelperPID->GetParticleSpecies(partMC);
 263       varSt[3]=(Double_t)partMC->IsPhysicalPrimary();
 264       varSt[4]=partMC->Y();
 265       varSt[5]=etaselected;
 266
 267       if(fQst==0) varSt[6]=Qvec;
 268       else varSt[6]=QvecMC;
 269
 270       ((THnSparseF*)fOutput->FindObject("NSparseHistSt"))->Fill(varSt);//stack loop
 271
 272       //Printf("a particle");
 273
 274     }
 275   }
 276
 277   //main loop on tracks
 278
 279   Int_t Nch = 0.;
 280
 281   for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
 282     AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
 283     if(!track) AliFatal("Not a standard AOD");
 284     if(fCharge != 0 && track->Charge() != fCharge) continue;//if fCharge != 0 only select fCharge
 285     if (!fTrackCuts->IsSelected(track,kTRUE)) continue; //track selection (rapidity selection NOT in the standard cuts)
 286     if(!fFillOnlyEvents){
 287       Int_t IDrec=fHelperPID->GetParticleSpecies(track,kTRUE);//id from detector
 288       Double_t y= track->Y(fHelperPID->GetMass((AliHelperParticleSpecies_t)IDrec));
 289       Int_t IDgen=kSpUndefined;//set if MC
 290       Int_t isph=-999;
 291       //       Int_t iswd=-999;
 292
 293       if (arrayMC) {
 294         AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));
 295         if (!partMC) {
 296           AliError("Cannot get MC particle");
 297           continue;
 298         }
 299         IDgen=fHelperPID->GetParticleSpecies(partMC);
 300         isph=partMC->IsPhysicalPrimary();
 301         //iswd=partMC->IsSecondaryFromWeakDecay();//FIXME not working on old productions - removed Apr 8th 2014
 302
 303         if(fIsAOD160){// enabled for new ADO160 only
 304           if(partMC->IsSecondaryFromWeakDecay()) isph=2.;
 305           if(partMC->IsSecondaryFromMaterial()) isph=3.;
 306
 307         }
 308       }
 309
 310       /*** DCA ***/
 311       Double_t dcaxy = -999.;
 312       Double_t dcaz = -999.;
 313
 314       Double_t p[2];
 315       if(GetDCA(track,p)){ dcaxy=p[0]; dcaz=p[1];}
 316
 317       if(dcaz >= fDCAzCut) continue;
 318
 319       //if the q vector is done using the TPC, we avoid overlap
 320       if (fVZEROside==2 && TMath::Abs(track->Eta())<0.5)continue;
 321
 322       //pt     cent    Q vec     IDrec     IDgen       isph      y
 323       //      Double_t varTrk[8];
 324       Double_t varTrk[7];
 325       varTrk[0]=track->Pt();
 326       varTrk[1]=Cent;
 327       if(fIsMC && fQtrk==1) varTrk[2]=QvecMC;
 328       else varTrk[2]=Qvec;
 329       varTrk[3]=(Double_t)IDrec;
 330       varTrk[4]=(Double_t)IDgen;
 331       varTrk[5]=(Double_t)isph;
 332       varTrk[6]=y;
 333       //varTrk[7]=dcaxy;
 334       ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
 335
 336       //for nsigma PID fill double counting of ID
 337       if(fHelperPID->GetPIDType()<kBayes && fDoDoubleCounting){//only nsigma
 338         Bool_t *HasDC;
 339         HasDC=fHelperPID->GetDoubleCounting(track,kTRUE);//get the array with double counting
 340         for(Int_t ipart=0;ipart<kNSpecies;ipart++){
 341           if(HasDC[ipart]==kTRUE){
 342             varTrk[3]=(Double_t)ipart;
 343             ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
 344           }
 345         }
 346       }
 347
 348       //fill all charged (3)
 349       varTrk[3]=3.;
 350       varTrk[4]=3.;
 351       ((THnSparseF*)fOutput->FindObject("NSparseHistTrk"))->Fill(varTrk);//track loop
 352     }//end if fFillOnlyEvents
 353
 354     //Printf("a track");
 355
 356     Nch++;
 357   } // end loop on tracks
 358
 359   Double_t varEv[4];
 360   varEv[0]=Cent;
 361   varEv[1]=Qvec;
 362   varEv[2]=Nch;
 363   ((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
 364
 365   PostData(1, fOutput  );
 366   PostData(2, fEventCuts);
 367   PostData(3, fTrackCuts);
 368   PostData(4, fHelperPID);
 369 }
 370
 371 //_________________________________________________________________
 372 Bool_t  AliAnalysisTaskSpectraAllChAOD::GetDCA(const AliAODTrack* trk, Double_t * p){
 373
 374   //AliAODTrack::DCA(): for newest AOD fTrack->DCA() always gives -999. This should fix.
 375   //FIXME should update EventCuts?
 376   //FIXME add track->GetXYZ(p) method
 377
 378   double xyz[3],cov[3];
 379
 380   if (!trk->GetXYZ(xyz)) { // dca is not stored
 381     AliExternalTrackParam etp;
 382     etp.CopyFromVTrack(trk);
 383     AliVEvent* ev = (AliVEvent*)trk->GetEvent();
 384     if (!ev) {/*printf("Event is not connected to the track\n");*/ return kFALSE;}
 385     if (!etp.PropagateToDCA(ev->GetPrimaryVertex(), ev->GetMagneticField(),999,xyz,cov)) return kFALSE; // failed, track is too far from vertex
 386   }
 387   p[0] = xyz[0];
 388   p[1] = xyz[1];
 389   return kTRUE;
 390
 391 }
 392
 393 //_________________________________________________________________
 394 void   AliAnalysisTaskSpectraAllChAOD::Terminate(Option_t *)
 395 {
 396   // Terminate
 397 }