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   fHnQAa(NULL), fHnQAb(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    binHnQAa[10] = {AliEbyEPidRatioHelper::fgkfHistNBinsCent, 4,
  72                           AliEbyEPidRatioHelper::fgkfHistNBinsSign,
  73                           AliEbyEPidRatioHelper::fgkfHistNBinsPt,
  74                            AliEbyEPidRatioHelper::fgkfHistNBinsPt,
  75                           500, 50, 50, 50, 3};
  76
  77   Double_t minHnQAa[10] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[0],-0.5,
  78                           AliEbyEPidRatioHelper::fgkfHistRangeSign[0],
  79                           AliEbyEPidRatioHelper::fgkfHistRangePt[0],
  80                            AliEbyEPidRatioHelper::fgkfHistRangePt[0],
  81                           30, -10,-10, -10, -0.5};
  82
  83   Double_t maxHnQAa[10] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[1], 3.5,
  84                           AliEbyEPidRatioHelper::fgkfHistRangeSign[1],
  85                           AliEbyEPidRatioHelper::fgkfHistRangePt[1],
  86                            AliEbyEPidRatioHelper::fgkfHistRangePt[1],
  87                           500, 10., 10., 10., 2.5};
  88
  89   Int_t    binHnQAb[9] = {AliEbyEPidRatioHelper::fgkfHistNBinsCent, 4,
  90                           AliEbyEPidRatioHelper::fgkfHistNBinsSign,
  91                           AliEbyEPidRatioHelper::fgkfHistNBinsPt,
  92                           AliEbyEPidRatioHelper::fgkfHistNBinsEta,
  93                           AliEbyEPidRatioHelper::fgkfHistNBinsRap,
  94                           AliEbyEPidRatioHelper::fgkfHistNBinsPhi,
  95                           50, 50};
  96
  97   Double_t minHnQAb[9] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[0],-0.5,
  98                            AliEbyEPidRatioHelper::fgkfHistRangeSign[0],
  99                           AliEbyEPidRatioHelper::fgkfHistRangePt[0],
 100                           AliEbyEPidRatioHelper::fgkfHistRangeEta[0],
 101                           AliEbyEPidRatioHelper::fgkfHistRangeRap[0],
 102                           AliEbyEPidRatioHelper::fgkfHistRangePhi[0],
 103                           -5,-5};
 104
 105   Double_t maxHnQAb[9] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[1], 3.5,
 106                           AliEbyEPidRatioHelper::fgkfHistRangeSign[1],
 107                           AliEbyEPidRatioHelper::fgkfHistRangePt[1],
 108                           AliEbyEPidRatioHelper::fgkfHistRangeEta[1],
 109                           AliEbyEPidRatioHelper::fgkfHistRangeRap[1],
 110                           AliEbyEPidRatioHelper::fgkfHistRangePhi[1],
 111                           5., 5.};
 112
 113   fHnQAa = new THnSparseF("hnQAPid", "cent:pid:sign:pt:pInner:TPCsignal:nSigmaITS:nSigmaTPC:nSigmaTOF:MCisProbe", 10, binHnQAa, minHnQAa, maxHnQAa);
 114   fHnQAb = new THnSparseF("hnQADca", "cent:pid:sign:pt:eta:y:phi:DCAr:DCAz", 9, binHnQAb, minHnQAb, maxHnQAb);
 115   fHnQAa->Sumw2();
 116   fHnQAb->Sumw2();
 117
 118   fHnQAa->GetAxis(0)->SetTitle("centrality");                   //  0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90|90-100 --> 11 bins
 119   fHnQAa->GetAxis(1)->SetTitle("N_{ch}|N_{#pi}|N_{K}|N_{p}");   //  0 | 1 | 2 | 3
 120   fHnQAa->GetAxis(2)->SetTitle("sign");                         //  -1 | 0 | +1
 121   fHnQAa->GetAxis(3)->SetTitle("#it{p}_{T} (GeV/#it{c})");      //  pt   [ 0.46, 2.22]
 122   fHnQAa->GetAxis(4)->SetTitle("#it{p}_{Inner} (GeV/#it{c})");  //  pInner [ 0.1, 3.0]
 123   fHnQAa->GetAxis(5)->SetTitle("TPC signal");                   //  TPCsignal [30, 500]
 124   fHnQAa->GetAxis(6)->SetTitle("n #sigma ITS");                 //  nSigma ITS [-10, 10]
 125   fHnQAa->GetAxis(7)->SetTitle("n #sigma TPC");                 //  nSigma TPC [-10, 10]
 126   fHnQAa->GetAxis(8)->SetTitle("n #sigma TOF");                 //  nSigma TOF [-10, 10]
 127   fHnQAa->GetAxis(9)->SetTitle("MCisProbe");                    //  0 | 1 (isProbeParticle) | 2 (isProbeParticle wrong sign)
 128
 129   fHnQAb->GetAxis(0)->SetTitle("centrality");                   //  0-5|5-10|10-20|20-30|30-40|40-50|50-60|60-70|70-80|80-90!90-100 --> 11 bins
 130   fHnQAb->GetAxis(1)->SetTitle("N_{ch}|N_{#pi}|N_{K}|N_{p}");   //  0 | 1 | 2 | 3
 131   fHnQAb->GetAxis(2)->SetTitle("sign");                         //  -1 | 0 | +1
 132   fHnQAb->GetAxis(3)->SetTitle("#it{p}_{T} (GeV/#it{c})");      //  pt   [ 0.46, 2.22]
 133   fHnQAb->GetAxis(4)->SetTitle("#eta");                         //  eta  [-0.9, 0.9]
 134   fHnQAb->GetAxis(5)->SetTitle("#it{y}");                       //  rapidity [-0.5, 0.5]
 135   fHnQAb->GetAxis(6)->SetTitle("#varphi");                      //  phi  [ 0. , 2Pi]
 136   fHnQAb->GetAxis(7)->SetTitle("DCAr");                         //  DCAr [-10, 10]
 137   fHnQAb->GetAxis(8)->SetTitle("DCAz");                         //  DCAz [-10, 10]
 138 // needed to be add later
 139   // fHnQA->GetAxis(9)->SetTitle("n #sigma #sqrt(Cdd)/DCAr");   //  nSigma DCAr [-10, 10]
 140   // fHnQA->GetAxis(10)->SetTitle("n #sigma #sqrt(Czz)/DCAz");  //  nSigma DCAz [-10, 10]
 141
 142   fHelper->BinLogAxis(fHnQAa, 3);
 143   fHelper->BinLogAxis(fHnQAb, 3);
 144
 145   return;
 146 }
 147
 148 //________________________________________________________________________
 149 void AliEbyEPidRatioQA::Process() {
 150   Float_t etaRange[2];
 151   fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
 152
 153   Float_t ptRange[2];
 154   fESDTrackCuts->GetPtRange(ptRange[0],ptRange[1]);
 155
 156   // -- Track Loop
 157   for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
 158
 159     AliVTrack *track = (fESD) ? static_cast<AliVTrack*>(fESD->GetTrack(idxTrack)) : static_cast<AliVTrack*>(fAOD->GetTrack(idxTrack));
 160
 161     // -- Check if track is accepted for basic parameters
 162     if (!fHelper->IsTrackAcceptedBasicCharged(track))
 163       continue;
 164
 165     // -- Check if accepted - ESD
 166     if (fESD){
 167       if (!fESDTrackCuts->AcceptTrack(dynamic_cast<AliESDtrack*>(track)))
 168         continue;
 169     }
 170
 171
 172     // -- Check if accepted - AOD
 173     if (fAOD){
 174       AliAODTrack * trackAOD = dynamic_cast<AliAODTrack*>(track);
 175
 176       if (!trackAOD) {
 177         AliError("Pointer to dynamic_cast<AliAODTrack*>(track) = ZERO");
 178         continue;
 179       }
 180       if (!trackAOD->TestFilterBit(fAODtrackCutBit))
 181         continue;
 182
 183       // -- Check if in pT and eta range (is done in ESDTrackCuts for ESDs)
 184       if(!(track->Pt() > ptRange[0] && track->Pt() <= ptRange[1] && TMath::Abs(track->Eta()) <= etaRange[1]))
 185         continue;
 186     }
 187
 188     Int_t gPdgCode = 0;
 189
 190     Int_t iPid = 0;
 191     Double_t pid[3];
 192     if      (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kPion)))  {  iPid = 1; gPdgCode = 211;}
 193     else if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kKaon)))  {  iPid = 2; gPdgCode = 321;}
 194     else if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kProton))){  iPid = 3; gPdgCode = 2212;}
 195     else iPid = 0;
 196
 197
 198     //cout << idxTrack << " --- QA ---- " << iPid << "  " << gPdgCode << endl;
 199
 200     Double_t yP;
 201     if ((iPid != 0) && !fHelper->IsTrackAcceptedRapidity(track, yP, iPid))
 202       continue;
 203
 204     if (!fHelper->IsTrackAcceptedDCA(track))
 205       continue;
 206
 207     // -- Check if probe particle
 208     Int_t isProbeParticle = 0;
 209     if (fIsMC) {
 210       Int_t label  = TMath::Abs(track->GetLabel());
 211
 212       AliVParticle* particle = (fESD) ? fMCEvent->GetTrack(label) : static_cast<AliVParticle*>(fArrayMC->At(label));
 213       if (particle) {
 214         if (TMath::Abs(particle->PdgCode()) == gPdgCode) {
 215           ++isProbeParticle;
 216           if (particle->PdgCode() != (track->Charge()*gPdgCode))
 217             ++isProbeParticle;
 218         }
 219       }
 220     }
 221
 222     // -- Get dca r/z
 223     Float_t dca[] = {0.,0.};     // dca_xy, dca_z,
 224     Float_t cov[] = {0.,0.,0.}; // sigma_xy, sigma_xy_z, sigma_z
 225     if (fESD)
 226       (static_cast<AliESDtrack*>(track))->GetImpactParameters(dca,cov);
 227
 228     //  Float_t dcaRoverCdd = ( TMath::Sqrt(cov[0]) != 0. )  ? dca[0]/TMath::Sqrt(cov[0]) : -9.99;
 229     //  Float_t dcaZoverCzz = ( TMath::Sqrt(cov[2]) != 0. )  ? dca[1]/TMath::Sqrt(cov[2]) : -9.99;
 230
 231     if (iPid == 0)
 232       yP = track->Eta();
 233
 234
 235     // cout << pid[0] << "  " << pid[1] << " " << pid[2] << yP << "  " << iPid << endl;
 236
 237     if (iPid != 0) {
 238
 239       Double_t aTracka[10] = {fCentralityBin,0,
 240                               static_cast<Double_t>(track->Charge()),
 241                               track->Pt(),
 242                               track->GetTPCmomentum(),track->GetTPCsignal(),pid[0],pid[1],pid[2],
 243                               static_cast<Double_t>(isProbeParticle)};
 244       Double_t aTrackb[9] = {fCentralityBin,0,
 245                              static_cast<Double_t>(track->Charge()),
 246                              track->Pt(),track->Eta(),yP,
 247                              track->Phi(),dca[0],dca[1]};
 248       fHnQAa->Fill(aTracka);
 249       fHnQAb->Fill(aTrackb);
 250     }
 251
 252     Double_t aTracka[10] = {fCentralityBin,
 253                             static_cast<Double_t>(iPid),
 254                             static_cast<Double_t>(track->Charge()),
 255                             track->Pt(),
 256                             track->GetTPCmomentum(),
 257                             track->GetTPCsignal(),pid[0],pid[1],pid[2],
 258                             static_cast<Double_t>(isProbeParticle)};
 259     Double_t aTrackb[9] = {fCentralityBin,
 260                            static_cast<Double_t>(iPid),
 261                            static_cast<Double_t>(track->Charge()),
 262                            track->Pt(),track->Eta(),yP, track->Phi(),dca[0],dca[1]};
 263       fHnQAa->Fill(aTracka);
 264       fHnQAb->Fill(aTrackb);
 265
 266   } // for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
 267
 268   return;
 269 }