PWGCF/EBYE/PIDFluctuation/task/AliEbyEPidTTask.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: ALICE Offline.                                                 *
   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 //=========================================================================//
  18 //             AliEbyE Analysis for Particle Ratio Fluctuation             //
  19 //                   Deepika Rathee  | Satyajit Jena                       //
  20 //                   drathee@cern.ch | sjena@cern.ch                       //
  21 //          Dealing with Wide pT Addition (Test Only - not for use)
  22 //=========================================================================//
  23
  24 #include "TChain.h"
  25 #include "TList.h"
  26 #include "TFile.h"
  27 #include "TTree.h"
  28 #include "TH1D.h"
  29 #include "TH2F.h"
  30 #include "TH3F.h"
  31 #include "TCanvas.h"
  32 #include "AliAnalysisTask.h"
  33 #include "AliAnalysisManager.h"
  34 #include "AliVEvent.h"
  35 #include "AliESD.h"
  36 #include "AliESDEvent.h"
  37 #include "AliAODEvent.h"
  38 #include "AliAODMCParticle.h"
  39 #include "AliAODMCHeader.h"
  40 #include "AliPIDResponse.h"
  41 #include "AliPIDCombined.h"
  42 #include "AliAODHeader.h"
  43 #include "AliAODpidUtil.h"
  44 #include "AliHelperPID.h"
  45 using std::endl;
  46 using std::cout;
  47 #include "AliEbyEPidTTask.h"
  48
  49 ClassImp(AliEbyEPidTTask)
  50
  51 //-----------------------------------------------------------------------
  52 AliEbyEPidTTask::AliEbyEPidTTask( const char *name ) : AliAnalysisTaskSE( name ),
  53   fThnList(NULL),
  54   isMC(kFALSE),
  55   fCentralityEstimator("V0M"),
  56   fAODtrackCutBit(128),
  57   fHelperPID(0x0),
  58   fEventCounter(NULL),
  59   fEventTree(NULL),
  60   isQA(kFALSE),
  61   fDebug(kFALSE),
  62   fRunNumber(0),
  63   fNumberOfTracks(8000),
  64   fNumberOfTracksM(8000),
  65   fFilterBit(99),
  66   fCentPercentile(-1),
  67   fVertexX(-99),
  68   fVertexY(-99),
  69   fVertexZ(-99) {
  70   /*
  71     fTrackPt[kTrack];
  72     fTrackEta[kTrack];
  73     fTrackPhi[kTrack];
  74     fTrackPtM[kTrack];
  75     fTrackEtaM[kTrack];
  76     fTrackPhiM[kTrack];
  77     fTrackCharge[kTrack];
  78     fTrackPid[kTrack];
  79     fTrackChargeM[kTrack];
  80     fTrackPidM[kTrack];
  81   */
  82
  83   DefineOutput(1, TList::Class()); //! Connect Outpput....
  84   DefineOutput(2, TTree::Class()); //! Connect Outpput....
  85 }
  86
  87 AliEbyEPidTTask::~AliEbyEPidTTask() {
  88   //!   Cleaning up
  89   if (fThnList)   delete fThnList;
  90   if (fHelperPID) delete fHelperPID;
  91   if (fEventTree) delete fEventTree;
  92 }
  93
  94 //---------------------------------------------------------------------------------
  95 void AliEbyEPidTTask::UserCreateOutputObjects() {
  96   fThnList = new TList();
  97   fThnList->SetOwner(kTRUE);
  98
  99   fEventCounter = new TH1D("fEventCounter","EventCounter", 70, 0.5,70.5);
 100   fThnList->Add(fEventCounter);
 101
 102   TList *ll = (TList*)fHelperPID->GetOutputList();
 103   for (Int_t ikey = 0; ikey < ll->GetEntries(); ikey++) {
 104     fThnList->Add(ll->At(ikey));
 105   }
 106   TDirectory *owd = gDirectory;
 107   OpenFile(1);
 108   fEventTree = new TTree("fEventTree","fEventTree");
 109   owd->cd();
 110
 111   fEventTree->Branch("fRunNumber",      &fRunNumber,     "fRunNumber/I");
 112   fEventTree->Branch("fFilterBit",      &fFilterBit,     "fFilterBit/I");
 113   fEventTree->Branch("fNumberOfTracks", &fNumberOfTracks,"fNumberOfTracks/I");
 114   fEventTree->Branch("fCentPercentile", &fCentPercentile,"fCentPercentile/F");
 115
 116   fEventTree->Branch("fVertexX",        &fVertexX,       "fVertexX/F");
 117   fEventTree->Branch("fVertexY",        &fVertexY,       "fVertexY/F");
 118   fEventTree->Branch("fVertexZ",        &fVertexZ,       "fVertexZ/F");
 119
 120   fEventTree->Branch("fTrackPt",        fTrackPt,        "fTrackPt[fNumberOfTracks]/F");
 121   fEventTree->Branch("fTrackPhi",       fTrackPhi,       "fTrackPhi[fNumberOfTracks]/F");
 122   fEventTree->Branch("fTrackEta",       fTrackEta,       "fTrackEta[fNumberOfTracks]/F");
 123   fEventTree->Branch("fTrackCharge",    fTrackCharge,    "fTrackCharge[fNumberOfTracks]/I");
 124   fEventTree->Branch("fTrackPid",       fTrackPid,       "fTrackPid[fNumberOfTracks]/I");
 125
 126   if(isMC) {
 127     fEventTree->Branch("fNumberOfTracksM", &fNumberOfTracksM, "fNumberOfTracksM/I");
 128     fEventTree->Branch("fTrackPtM",        fTrackPtM,        "fTrackPtM[fNumberOfTracksM]/F");
 129     fEventTree->Branch("fTrackPhiM",       fTrackPhiM,       "fTrackPhiM[fNumberOfTracksM]/F");
 130     fEventTree->Branch("fTrackEtaM",       fTrackEtaM,       "fTrackEtaM[fNumberOfTracksM]/F");
 131     fEventTree->Branch("fTrackChargeM",    fTrackChargeM,    "fTrackChargeM[fNumberOfTracksM]/I");
 132     fEventTree->Branch("fTrackPidM",       fTrackPidM,       "fTrackPidM[fNumberOfTracksM]/I");
 133   }
 134
 135   PostData(1, fThnList);
 136   PostData(2, fEventTree);
 137 }
 138
 139 //----------------------------------------------------------------------------------
 140 void AliEbyEPidTTask::UserExec( Option_t * ){
 141
 142    fEventCounter->Fill(1);
 143
 144   AliAODEvent* event = dynamic_cast<AliAODEvent*>(InputEvent());
 145   if (!event) {
 146     Printf("ERROR 01: AOD not found ");
 147     return;
 148   }
 149
 150   fEventCounter->Fill(2);
 151
 152   Int_t gCent   = -1;
 153  // Float_t gRefMul = -1;
 154
 155   AliAODHeader *aodHeader = dynamic_cast<AliAODHeader*>(event->GetHeader());
 156   if(!aodHeader) AliFatal("Not a standard AOD");
 157   gCent = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
 158  // gRefMul = aodHeader->GetRefMultiplicity();
 159   if (gCent < 0 || gCent > 100) return;
 160   fEventCounter->Fill(3);
 161
 162   const AliAODVertex *vertex = event->GetPrimaryVertex();
 163   if(!vertex) return;
 164   fEventCounter->Fill(4);
 165   Bool_t vtest = kFALSE;
 166   Double32_t fCov[6];
 167   vertex->GetCovarianceMatrix(fCov);
 168   if(vertex->GetNContributors() > 0) {
 169       if(fCov[5] != 0) {
 170         vtest = kTRUE;
 171       }
 172   }
 173   if(!vtest)return;
 174
 175   fEventCounter->Fill(5);
 176
 177   AliCentrality *centrality = event->GetCentrality();
 178   if (centrality->GetQuality() != 0) return;
 179
 180   fEventCounter->Fill(6);
 181
 182   fRunNumber = event->GetRunNumber();
 183   fFilterBit = fAODtrackCutBit;
 184   fCentPercentile = gCent;
 185   fVertexX = vertex->GetX();
 186   fVertexY = vertex->GetY();
 187   fVertexZ = vertex->GetZ();
 188
 189 //Default TPC priors
 190   if(fHelperPID->GetPIDType()==kBayes)fHelperPID->GetPIDCombined()->SetDefaultTPCPriors();//FIXME maybe this can go in the UserCreateOutputObject?
 191
 192
 193
 194
 195   Int_t iTracks = 0;
 196   for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); itrk++) {
 197     AliAODTrack* track = dynamic_cast<AliAODTrack *>(event->GetTrack(itrk));
 198     fEventCounter->Fill(10);
 199     if (!track) continue;
 200     fEventCounter->Fill(11);
 201     if (!AcceptTrack(track)) continue;
 202     fEventCounter->Fill(12);
 203     Int_t a = fHelperPID->GetParticleSpecies((AliVTrack*) track,kTRUE);
 204     //  if(a < 0 || a > 2) continue;
 205     fEventCounter->Fill(13);
 206
 207     Int_t b = -999;
 208     if (a == 0 ) b = 1;
 209     else if (a == 1 ) b = 2;
 210     else if (a == 2 ) b = 3;
 211     else b = 4;
 212
 213     if (track->Charge()  < 0 ) b = -1*b;
 214
 215     //    Int_t icharge = track->Charge() > 0 ? 0 : 1;
 216
 217     // cout << icharge << "  " << track->Charge() << endl;
 218
 219     fTrackPt[iTracks]     = (Float_t)track->Pt();
 220     fTrackPhi[iTracks]    = (Float_t)track->Phi();
 221     fTrackEta[iTracks]    = (Float_t)track->Eta();
 222     fTrackCharge[iTracks] = track->Charge();
 223     fTrackPid[iTracks] = b;
 224     iTracks++;
 225   }
 226   fNumberOfTracks = iTracks;
 227
 228   if(isMC) {
 229     fEventCounter->Fill(21);
 230     TClonesArray *arrayMC= 0;
 231     arrayMC = dynamic_cast<TClonesArray*> (event->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
 232     if (!arrayMC) {
 233       Printf("Error: MC particles branch not found!\n");
 234       return;
 235     }
 236     fEventCounter->Fill(22);
 237     AliAODMCHeader *mcHdr=0;
 238     mcHdr=(AliAODMCHeader*)event->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 239     if(!mcHdr) {
 240       Printf("MC header branch not found!\n");
 241       return;
 242     }
 243
 244     fEventCounter->Fill(23);
 245
 246     Int_t nMC = arrayMC->GetEntries();
 247     Int_t mTracks = 0;
 248     for (Int_t iMC = 0; iMC < nMC; iMC++) {
 249       fEventCounter->Fill(24);
 250       AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
 251       if(!AcceptMCTrack(partMC)) continue;
 252
 253       fEventCounter->Fill(25);
 254       Int_t a  = partMC->GetPdgCode();
 255
 256       // Int_t a = fHelperPID->GetMCParticleSpecie((AliVEvent*) event,(AliVTrack*)partMC,1);
 257
 258       //if(a < 0 || a > 2) continue;
 259       Int_t icharge = a > 0 ? 0 : 1;
 260
 261       //  cout << a << "  " << icharge << endl;
 262
 263       fTrackPtM[mTracks]     = (Float_t)partMC->Pt();
 264       fTrackPhiM[mTracks]    = (Float_t)partMC->Phi();
 265       fTrackEtaM[mTracks]    = (Float_t)partMC->Eta();
 266       fTrackChargeM[mTracks] = icharge;
 267       fTrackPidM[mTracks] = a;
 268
 269       mTracks++;
 270     }
 271     fEventCounter->Fill(26);
 272     fNumberOfTracksM = mTracks;
 273   }
 274
 275   fEventCounter->Fill(30);
 276   fEventTree->Fill();
 277
 278   PostData(1, fThnList);
 279   PostData(2, fEventTree);
 280 }
 281
 282 //___________________________________________________________
 283 void AliEbyEPidTTask::Terminate( Option_t * ){
 284   Info("AliEbyEPidTTask"," Task Successfully finished");
 285 }
 286
 287 //___________________________________________________________
 288 Bool_t AliEbyEPidTTask::AcceptTrack(AliAODTrack *track) const {
 289   if(!track)                                  return kFALSE;
 290   if (track->Charge() == 0 )                  return kFALSE;
 291   if (!track->TestFilterBit(fAODtrackCutBit)) return kFALSE;
 292   return kTRUE;
 293 }
 294
 295
 296 //___________________________________________________________
 297 Bool_t AliEbyEPidTTask::AcceptMCTrack(AliAODMCParticle *track) const {
 298   if(!track) return kFALSE;
 299   if(!track->IsPhysicalPrimary()) return kFALSE;
 300   if (track->Charge() == 0 ) return kFALSE;
 301   return kTRUE;
 302 }