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