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 = event->GetHeader();
 156   gCent = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
 157  // gRefMul = aodHeader->GetRefMultiplicity();
 158   if (gCent < 0 || gCent > 100) return;
 159   fEventCounter->Fill(3);
 160
 161   const AliAODVertex *vertex = event->GetPrimaryVertex();
 162   if(!vertex) return;
 163   fEventCounter->Fill(4);
 164   Bool_t vtest = kFALSE;
 165   Double32_t fCov[6];
 166   vertex->GetCovarianceMatrix(fCov);
 167   if(vertex->GetNContributors() > 0) {
 168       if(fCov[5] != 0) {
 169         vtest = kTRUE;
 170       }
 171   }
 172   if(!vtest)return;
 173
 174   fEventCounter->Fill(5);
 175
 176   AliCentrality *centrality = event->GetCentrality();
 177   if (centrality->GetQuality() != 0) return;
 178
 179   fEventCounter->Fill(6);
 180
 181   fRunNumber = event->GetRunNumber();
 182   fFilterBit = fAODtrackCutBit;
 183   fCentPercentile = gCent;
 184   fVertexX = vertex->GetX();
 185   fVertexY = vertex->GetY();
 186   fVertexZ = vertex->GetZ();
 187
 188 //Default TPC priors
 189   if(fHelperPID->GetPIDType()==kBayes)fHelperPID->GetPIDCombined()->SetDefaultTPCPriors();//FIXME maybe this can go in the UserCreateOutputObject?
 190
 191
 192
 193
 194   Int_t iTracks = 0;
 195   for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); itrk++) {
 196     AliAODTrack* track = dynamic_cast<AliAODTrack *>(event->GetTrack(itrk));
 197     fEventCounter->Fill(10);
 198     if (!track) continue;
 199     fEventCounter->Fill(11);
 200     if (!AcceptTrack(track)) continue;
 201     fEventCounter->Fill(12);
 202     Int_t a = fHelperPID->GetParticleSpecies((AliVTrack*) track,kTRUE);
 203     //  if(a < 0 || a > 2) continue;
 204     fEventCounter->Fill(13);
 205
 206     Int_t b = -999;
 207     if (a == 0 ) b = 1;
 208     else if (a == 1 ) b = 2;
 209     else if (a == 2 ) b = 3;
 210     else b = 4;
 211
 212     if (track->Charge()  < 0 ) b = -1*b;
 213
 214     //    Int_t icharge = track->Charge() > 0 ? 0 : 1;
 215
 216     // cout << icharge << "  " << track->Charge() << endl;
 217
 218     fTrackPt[iTracks]     = (Float_t)track->Pt();
 219     fTrackPhi[iTracks]    = (Float_t)track->Phi();
 220     fTrackEta[iTracks]    = (Float_t)track->Eta();
 221     fTrackCharge[iTracks] = track->Charge();
 222     fTrackPid[iTracks] = b;
 223     iTracks++;
 224   }
 225   fNumberOfTracks = iTracks;
 226
 227   if(isMC) {
 228     fEventCounter->Fill(21);
 229     TClonesArray *arrayMC= 0;
 230     arrayMC = dynamic_cast<TClonesArray*> (event->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
 231     if (!arrayMC) {
 232       Printf("Error: MC particles branch not found!\n");
 233       return;
 234     }
 235     fEventCounter->Fill(22);
 236     AliAODMCHeader *mcHdr=0;
 237     mcHdr=(AliAODMCHeader*)event->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 238     if(!mcHdr) {
 239       Printf("MC header branch not found!\n");
 240       return;
 241     }
 242
 243     fEventCounter->Fill(23);
 244
 245     Int_t nMC = arrayMC->GetEntries();
 246     Int_t mTracks = 0;
 247     for (Int_t iMC = 0; iMC < nMC; iMC++) {
 248       fEventCounter->Fill(24);
 249       AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
 250       if(!AcceptMCTrack(partMC)) continue;
 251
 252       fEventCounter->Fill(25);
 253       Int_t a  = partMC->GetPdgCode();
 254
 255       // Int_t a = fHelperPID->GetMCParticleSpecie((AliVEvent*) event,(AliVTrack*)partMC,1);
 256
 257       //if(a < 0 || a > 2) continue;
 258       Int_t icharge = a > 0 ? 0 : 1;
 259
 260       //  cout << a << "  " << icharge << endl;
 261
 262       fTrackPtM[mTracks]     = (Float_t)partMC->Pt();
 263       fTrackPhiM[mTracks]    = (Float_t)partMC->Phi();
 264       fTrackEtaM[mTracks]    = (Float_t)partMC->Eta();
 265       fTrackChargeM[mTracks] = icharge;
 266       fTrackPidM[mTracks] = a;
 267
 268       mTracks++;
 269     }
 270     fEventCounter->Fill(26);
 271     fNumberOfTracksM = mTracks;
 272   }
 273
 274   fEventCounter->Fill(30);
 275   fEventTree->Fill();
 276
 277   PostData(1, fThnList);
 278   PostData(2, fEventTree);
 279 }
 280
 281 //___________________________________________________________
 282 void AliEbyEPidTTask::Terminate( Option_t * ){
 283   Info("AliEbyEPidTTask"," Task Successfully finished");
 284 }
 285
 286 //___________________________________________________________
 287 Bool_t AliEbyEPidTTask::AcceptTrack(AliAODTrack *track) const {
 288   if(!track)                                  return kFALSE;
 289   if (track->Charge() == 0 )                  return kFALSE;
 290   if (!track->TestFilterBit(fAODtrackCutBit)) return kFALSE;
 291   return kTRUE;
 292 }
 293
 294
 295 //___________________________________________________________
 296 Bool_t AliEbyEPidTTask::AcceptMCTrack(AliAODMCParticle *track) const {
 297   if(!track) return kFALSE;
 298   if(!track->IsPhysicalPrimary()) return kFALSE;
 299   if (track->Charge() == 0 ) return kFALSE;
 300   return kTRUE;
 301 }