PWGCF/EBYE/PIDFluctuation/task/AliEbyEPidRatioQA.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 //                  Date: Wed Jul  9 18:38:30 CEST 2014                    //
  22 //          New approch to find particle ratio to reduce memory            //
  23 //                             (Test Only)                                 //
  24 //=========================================================================//
  25
  26 #include "TMath.h"
  27 #include "TAxis.h"
  28 #include "TSystem.h"
  29 #include "TProfile.h"
  30 #include "TH2F.h"
  31 #include "TH3F.h"
  32 #include "TFile.h"
  33 #include "TPRegexp.h"
  34 #include "AliStack.h"
  35 #include "AliMCEvent.h"
  36 #include "AliMCParticle.h"
  37 #include "AliESDtrackCuts.h"
  38 #include "AliESDInputHandler.h"
  39 #include "AliESDpid.h"
  40 #include "AliCentrality.h"
  41 #include "AliTracker.h"
  42 #include "AliAODInputHandler.h"
  43 #include "AliAODEvent.h"
  44 #include "AliAODTrack.h"
  45 #include "AliAODMCParticle.h"
  46 #include "AliEbyEPidRatioQA.h"
  47
  48 using namespace std;
  49
  50 ClassImp(AliEbyEPidRatioQA)
  51 //________________________________________________________________________
  52 AliEbyEPidRatioQA::AliEbyEPidRatioQA() :
  53   AliEbyEPidRatioBase("QA", "QA"),
  54
  55   fHnQA(NULL) {
  56   // Constructor
  57
  58   AliLog::SetClassDebugLevel("AliEbyEPidRatioQA",10);
  59 }
  60
  61 //________________________________________________________________________
  62 AliEbyEPidRatioQA::~AliEbyEPidRatioQA() {
  63   // Destructor
  64
  65   return;
  66 }
  67
  68 //________________________________________________________________________
  69 void AliEbyEPidRatioQA::CreateHistograms() {
  70
  71   Int_t    binHnQA[17] = {AliEbyEPidRatioHelper::fgkfHistNBinsCent, AliEbyEPidRatioHelper::fgkfHistNBinsEta,  //       cent |       eta
  72                           AliEbyEPidRatioHelper::fgkfHistNBinsRap,  AliEbyEPidRatioHelper::fgkfHistNBinsPhi,  //          y |       phi
  73                           AliEbyEPidRatioHelper::fgkfHistNBinsPt,   AliEbyEPidRatioHelper::fgkfHistNBinsPt,   //         pt |    pInner
  74                           AliEbyEPidRatioHelper::fgkfHistNBinsSign, 500,                                      //       sign | TPCsignal |
  75                           50, 50, 50,                                                                         //  nSigmaITS | nSigmaTPC |  nSigmaTOF
  76                           50, 50, 50, 50,                                                                     //  DCAr      |      DCAz | nSigmaDCAr | nSigmaDCAz
  77                           3,4};                                                                                 //  MCisProbe
  78
  79   Double_t minHnQA[17] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[0], AliEbyEPidRatioHelper::fgkfHistRangeEta[0],
  80                           AliEbyEPidRatioHelper::fgkfHistRangeRap[0],  AliEbyEPidRatioHelper::fgkfHistRangePhi[0],
  81                           AliEbyEPidRatioHelper::fgkfHistRangePt[0],   AliEbyEPidRatioHelper::fgkfHistRangePt[0],
  82                           AliEbyEPidRatioHelper::fgkfHistRangeSign[0], 30,
  83                           -10., -10., -10.,
  84                           -10., -10., -10., -10.,
  85                           -0.5,-0.5};
  86
  87   Double_t maxHnQA[17] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[1], AliEbyEPidRatioHelper::fgkfHistRangeEta[1],
  88                           AliEbyEPidRatioHelper::fgkfHistRangeRap[1],  AliEbyEPidRatioHelper::fgkfHistRangePhi[1],
  89                           AliEbyEPidRatioHelper::fgkfHistRangePt[1],   AliEbyEPidRatioHelper::fgkfHistRangePt[1],
  90                           AliEbyEPidRatioHelper::fgkfHistRangeSign[1], 500,
  91                           10., 10., 10.,
  92                           10.,  10., 10., 10.,
  93                           2.5,3.5};
  94
  95
  96   fHnQA = new THnSparseF("hnQA", "cent:eta:y:phi:pt:pInner:sign:TPCsignal:nSigmaITS:nSigmaTPC:nSigmaTOF:DCAr:DCAz:nSigmaDCAr:nSigmaDCAz:MCisProbe",
  97                                  17, binHnQA, minHnQA, maxHnQA);
  98
  99   fHnQA->Sumw2();
 100
 101   fHnQA->GetAxis(0)->SetTitle("centrality");                   //  0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90 --> 10 bins
 102   fHnQA->GetAxis(1)->SetTitle("#eta");                         //  eta  [-0.9, 0.9]
 103   fHnQA->GetAxis(2)->SetTitle("#it{y}");                       //  rapidity [-0.5, 0.5]
 104   fHnQA->GetAxis(3)->SetTitle("#varphi");                      //  phi  [ 0. , 2Pi]
 105   fHnQA->GetAxis(4)->SetTitle("#it{p}_{T} (GeV/#it{c})");      //  pt   [ 0.46, 2.22]
 106   fHnQA->GetAxis(5)->SetTitle("#it{p}_{Inner} (GeV/#it{c})");  //  pInner [ 0.1, 3.0]
 107   fHnQA->GetAxis(6)->SetTitle("sign");                         //  -1 | 0 | +1
 108
 109   fHnQA->GetAxis(7)->SetTitle("TPC signal");                   //  TPCsignal [30, 500]
 110   fHnQA->GetAxis(8)->SetTitle("n #sigma ITS");                 //  nSigma ITS [-10, 10]
 111   fHnQA->GetAxis(9)->SetTitle("n #sigma TPC");                 //  nSigma TPC [-10, 10]
 112   fHnQA->GetAxis(10)->SetTitle("n #sigma TOF");                //  nSigma TOF [-10, 10]
 113
 114   fHnQA->GetAxis(11)->SetTitle("DCAr");                        //  DCAr [-10, 10]
 115   fHnQA->GetAxis(12)->SetTitle("DCAz");                        //  DCAz [-10, 10]
 116   fHnQA->GetAxis(13)->SetTitle("n #sigma #sqrt(Cdd)/DCAr");    //  nSigma DCAr [-10, 10]
 117   fHnQA->GetAxis(14)->SetTitle("n #sigma #sqrt(Czz)/DCAz");    //  nSigma DCAz [-10, 10]
 118   fHnQA->GetAxis(15)->SetTitle("MCisProbe");                   //  0 | 1 (isProbeParticle) | 2 (isProbeParticle wrong sign)
 119   fHnQA->GetAxis(16)->SetTitle("PID");                         //  0 | 1 | 2 | 3
 120
 121   fHelper->BinLogAxis(fHnQA, 4);
 122   fHelper->BinLogAxis(fHnQA, 5);
 123
 124   return;
 125 }
 126
 127 //________________________________________________________________________
 128 void AliEbyEPidRatioQA::Process() {
 129   Float_t etaRange[2];
 130   fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
 131
 132   Float_t ptRange[2];
 133   fESDTrackCuts->GetPtRange(ptRange[0],ptRange[1]);
 134
 135   // -- Track Loop
 136   for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
 137
 138     AliVTrack *track = (fESD) ? static_cast<AliVTrack*>(fESD->GetTrack(idxTrack)) : static_cast<AliVTrack*>(fAOD->GetTrack(idxTrack));
 139
 140     // -- Check if track is accepted for basic parameters
 141     if (!fHelper->IsTrackAcceptedBasicCharged(track))
 142       continue;
 143
 144     // -- Check if accepted - ESD
 145     if (fESD && !fESDTrackCuts->AcceptTrack(dynamic_cast<AliESDtrack*>(track)))
 146       continue;
 147
 148     // -- Check if accepted - AOD
 149     if (fAOD){
 150       AliAODTrack * trackAOD = dynamic_cast<AliAODTrack*>(track);
 151
 152       if (!trackAOD) {
 153         AliError("Pointer to dynamic_cast<AliAODTrack*>(track) = ZERO");
 154         continue;
 155       }
 156       if (!trackAOD->TestFilterBit(fAODtrackCutBit))
 157         continue;
 158
 159       // -- Check if in pT and eta range (is done in ESDTrackCuts for ESDs)
 160       if(!(track->Pt() > ptRange[0] && track->Pt() <= ptRange[1] && TMath::Abs(track->Eta()) <= etaRange[1]))
 161         continue;
 162     }
 163
 164     Int_t gPdgCode = 0;
 165
 166     Int_t iPid = 0;
 167     Double_t pid[3];
 168     if      (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kPion)))  {  iPid = 1; gPdgCode = 211;}
 169     else if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kKaon)))  {  iPid = 2; gPdgCode = 321;}
 170     else if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kProton))){  iPid = 3; gPdgCode = 2212;}
 171     else iPid = 0;
 172
 173     Double_t yP;
 174     if ((iPid != 0) && !fHelper->IsTrackAcceptedRapidity(track, yP, iPid))
 175       continue;
 176
 177     if (!fHelper->IsTrackAcceptedDCA(track))
 178       continue;
 179
 180     // -- Check if probe particle
 181     Int_t isProbeParticle = 0;
 182     if (fIsMC) {
 183       Int_t label  = TMath::Abs(track->GetLabel());
 184
 185       AliVParticle* particle = (fESD) ? fMCEvent->GetTrack(label) : static_cast<AliVParticle*>(fArrayMC->At(label));
 186       if (particle) {
 187         if (TMath::Abs(particle->PdgCode()) == gPdgCode) {
 188           ++isProbeParticle;
 189           if (particle->PdgCode() != (track->Charge()*gPdgCode))
 190             ++isProbeParticle;
 191         }
 192       }
 193     }
 194
 195     // -- Get dca r/z
 196     Float_t dca[] = {0.,0.};     // dca_xy, dca_z,
 197     Float_t cov[] = {0.,0.,0.}; // sigma_xy, sigma_xy_z, sigma_z
 198     if (fESD)
 199       (static_cast<AliESDtrack*>(track))->GetImpactParameters(dca,cov);
 200
 201     Float_t dcaRoverCdd = ( TMath::Sqrt(cov[0]) != 0. )  ? dca[0]/TMath::Sqrt(cov[0]) : -9.99;
 202     Float_t dcaZoverCzz = ( TMath::Sqrt(cov[2]) != 0. )  ? dca[1]/TMath::Sqrt(cov[2]) : -9.99;
 203
 204     if (iPid == 0)
 205       yP = track->Eta();
 206
 207
 208     // cout << pid[0] << "  " << pid[1] << " " << pid[2] << yP << "  " << iPid << endl;
 209
 210     if (iPid != 0) {
 211       Double_t aTrack[17] = {
 212         Double_t(fCentralityBin),               //  0 centrality
 213         track->Eta(),                           //  1 eta
 214         yP,                                     //  2 rapidity
 215         track->Phi(),                           //  3 phi
 216         track->Pt(),                            //  4 pt
 217         track->GetTPCmomentum(),                //  5 pInner
 218         track->Charge(),                        //  6 sign
 219         track->GetTPCsignal(),                  //  7 TPC dE/dx
 220         pid[0],                                 //  8 n Sigma ITS
 221         pid[1],                                 //  9 n Sigma TPC
 222         pid[2],                                 // 10 n Sigma TOF
 223         dca[0],                                 // 11 dca r
 224         dca[1],                                 // 12 dca z
 225         dcaRoverCdd,                            // 13 sqrt(cov[dd])
 226         dcaZoverCzz,                            // 14 sqrt(cov[zz])
 227         isProbeParticle,                        // 15 isProbe
 228         0                                    // 16 PID
 229       };
 230
 231       fHnQA->Fill(aTrack);
 232     }
 233
 234     Double_t aTrack[17] = {
 235       Double_t(fCentralityBin),               //  0 centrality
 236       track->Eta(),                           //  1 eta
 237       yP,                                     //  2 rapidity
 238       track->Phi(),                           //  3 phi
 239       track->Pt(),                            //  4 pt
 240       track->GetTPCmomentum(),                //  5 pInner
 241       track->Charge(),                        //  6 sign
 242       track->GetTPCsignal(),                  //  7 TPC dE/dx
 243       pid[0],                                 //  8 n Sigma ITS
 244       pid[1],                                 //  9 n Sigma TPC
 245       pid[2],                                 // 10 n Sigma TOF
 246       dca[0],                                 // 11 dca r
 247       dca[1],                                 // 12 dca z
 248       dcaRoverCdd,                            // 13 sqrt(cov[dd])
 249       dcaZoverCzz,                            // 14 sqrt(cov[zz])
 250       isProbeParticle,                        // 15 isProbe
 251       iPid                                    // 16 PID
 252     };
 253
 254    fHnQA->Fill(aTrack);
 255
 256   } // for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
 257
 258   return;
 259 }