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