PWGCF/EBYE/PIDFluctuation/task/AliEbyEPidRatioDCA.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 //        Copied from NetParticle Classes
  25 //        Origin: Authors: Jochen Thaeder <jochen@thaeder.de>
  26 //                         Michael Weber <m.weber@cern.ch>
  27 //=========================================================================//
  28
  29 #include "TMath.h"
  30 #include "TAxis.h"
  31
  32 #include "AliESDEvent.h"
  33 #include "AliStack.h"
  34 #include "AliMCEvent.h"
  35 #include "AliESDtrackCuts.h"
  36
  37 #include "AliAODEvent.h"
  38 #include "AliAODMCParticle.h"
  39
  40 #include "AliEbyEPidRatioDCA.h"
  41
  42 using namespace std;
  43
  44 ClassImp(AliEbyEPidRatioDCA)
  45
  46 //________________________________________________________________________
  47 AliEbyEPidRatioDCA::AliEbyEPidRatioDCA() :
  48   AliEbyEPidRatioBase("DCA", "DCA"),
  49   fESDTrackCutsBkg(NULL),
  50   fHnDCA(NULL) {
  51   // Constructor
  52
  53
  54   AliLog::SetClassDebugLevel("AliEbyEPidRatioDCA",10);
  55 }
  56
  57 //________________________________________________________________________
  58 AliEbyEPidRatioDCA::~AliEbyEPidRatioDCA() {
  59   // Destructor
  60 }
  61
  62 //________________________________________________________________________
  63 void AliEbyEPidRatioDCA::Process() {
  64    Float_t etaRange[2];
  65   fESDTrackCutsBkg->GetEtaRange(etaRange[0],etaRange[1]);
  66
  67   Float_t ptRange[2];
  68   fESDTrackCutsBkg->GetPtRange(ptRange[0],ptRange[1]);
  69
  70   // -- Track Loop
  71   for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
  72
  73     AliVTrack *track = (fESD) ? static_cast<AliVTrack*>(fESD->GetTrack(idxTrack)) : static_cast<AliVTrack*>(fAOD->GetTrack(idxTrack));
  74
  75     // -- Check if track is accepted for basic parameters
  76     if (!fHelper->IsTrackAcceptedBasicCharged(track))
  77       continue;
  78
  79     // -- Check if accepted - ESD
  80     if (fESD && !fESDTrackCutsBkg->AcceptTrack(dynamic_cast<AliESDtrack*>(track)))
  81       continue;
  82
  83     // -- Check if accepted - AOD
  84     if (fAOD){
  85       AliAODTrack * trackAOD = dynamic_cast<AliAODTrack*>(track);
  86
  87       if (!trackAOD) {
  88         AliError("Pointer to dynamic_cast<AliAODTrack*>(track) = ZERO");
  89         continue;
  90       }
  91       if (!trackAOD->TestFilterBit(fAODtrackCutBit))
  92         continue;
  93
  94       // -- Check if in pT and eta range (is done in ESDTrackCuts for ESDs)
  95       if(!(track->Pt() > ptRange[0] && track->Pt() <= ptRange[1] && TMath::Abs(track->Eta()) <= etaRange[1]))
  96         continue;
  97     }
  98
  99     Int_t gPdgCode = 0;
 100     Int_t iPid = 0;
 101     Double_t pid[3];
 102     if      (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kPion)))  {  iPid = 1; gPdgCode = 211;}
 103     else if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kKaon)))  {  iPid = 2; gPdgCode = 321;}
 104     else if (fHelper->IsTrackAcceptedPID(track, pid, (AliPID::kProton))){  iPid = 3; gPdgCode = 2212;}
 105     else { iPid = 0; gPdgCode = 0;}
 106
 107     //  cout << " --- DCA ---- " << iPid << "  " << gPdgCode << endl;
 108     Double_t yP;
 109     if (!fHelper->IsTrackAcceptedRapidity(track, yP, iPid))
 110       continue;
 111
 112     Bool_t isDCArAccepted = fHelper->IsTrackAcceptedDCA(track);
 113
 114     // -- Check if accepted with thighter DCA cuts
 115     // ?!?!? How to mimic this in AODs?
 116     if (fESD && !fESDTrackCuts->AcceptTrack(dynamic_cast<AliESDtrack*>(track)))
 117       isDCArAccepted = kFALSE;
 118
 119     // -- Check for contamination
 120     Int_t contIdx = (fIsMC) ? GetContIdxTrack(TMath::Abs(track->GetLabel()), track->Charge(), gPdgCode) : 1;
 121
 122     // -- Get DCAs (dca_r, dca_z, sigma_xy, sigma_xy_z, sigma_z)
 123     Float_t dca[2], cov[3]; //
 124     if (fESD)
 125       (dynamic_cast<AliESDtrack*>(track))->GetImpactParameters(dca, cov);
 126     else
 127       dca[0] = 1.;
 128
 129     // -- Fill THnSparse
 130
 131     if(iPid != 0) {
 132       Double_t hnDCA[10] = {fCentralityBin,0.,
 133                             static_cast<Double_t>(track->Charge()),
 134                             track->Eta(),
 135                             yP,
 136                             track->Phi(),
 137                             track->Pt(),
 138                             static_cast<Double_t>(contIdx),
 139                             static_cast<Double_t>(isDCArAccepted),
 140                             dca[0]};
 141       fHnDCA->Fill(hnDCA);
 142     }
 143
 144     Double_t hnDCA[10] = {fCentralityBin,
 145                           static_cast<Double_t>(iPid),
 146                           static_cast<Double_t>(track->Charge()),
 147                           track->Eta(),
 148                           yP,
 149                           track->Phi(),
 150                           track->Pt(),
 151                           static_cast<Double_t>(contIdx),
 152                           static_cast<Double_t>(isDCArAccepted),
 153                           dca[0]};
 154       fHnDCA->Fill(hnDCA);
 155
 156
 157     } // for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
 158
 159   return;
 160 }
 161
 162 //________________________________________________________________________
 163 void AliEbyEPidRatioDCA::CreateHistograms() {
 164   Int_t    binHnDCA[10] = {AliEbyEPidRatioHelper::fgkfHistNBinsCent,4,
 165                            AliEbyEPidRatioHelper::fgkfHistNBinsSign,
 166                            AliEbyEPidRatioHelper::fgkfHistNBinsEta,
 167                            AliEbyEPidRatioHelper::fgkfHistNBinsRap,
 168                            AliEbyEPidRatioHelper::fgkfHistNBinsPhi,
 169                            AliEbyEPidRatioHelper::fgkfHistNBinsPt,
 170                            4,  2,  77};
 171
 172   Double_t minHnDCA[10] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[0],-0.5,
 173                            AliEbyEPidRatioHelper::fgkfHistRangeSign[0],
 174                            AliEbyEPidRatioHelper::fgkfHistRangeEta[0],
 175                            AliEbyEPidRatioHelper::fgkfHistRangeRap[0],
 176                            AliEbyEPidRatioHelper::fgkfHistRangePhi[0],
 177                            AliEbyEPidRatioHelper::fgkfHistRangePt[0],
 178                            0.5, -0.5, -3.};
 179
 180   Double_t maxHnDCA[10] = {AliEbyEPidRatioHelper::fgkfHistRangeCent[1],3.5,
 181                            AliEbyEPidRatioHelper::fgkfHistRangeSign[1],
 182                            AliEbyEPidRatioHelper::fgkfHistRangeEta[1],
 183                            AliEbyEPidRatioHelper::fgkfHistRangeRap[1],
 184                            AliEbyEPidRatioHelper::fgkfHistRangePhi[1],
 185                            AliEbyEPidRatioHelper::fgkfHistRangePt[1],
 186                            4.5, 1.5, 3.};
 187
 188
 189   fHnDCA = new THnSparseD("hnDCA", "cent:pid:etaRec:yRec:phiRec:ptRec:sign:contPart:contSign:DCArAccepted:DCAr", 10, binHnDCA, minHnDCA, maxHnDCA);
 190
 191   fHnDCA->Sumw2();
 192
 193   fHnDCA->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
 194   fHnDCA->GetAxis(1)->SetTitle("N_{ch}|N_{#pi}|N_{K}|N_{p}");   //  0 | 1 | 2 | 3
 195   fHnDCA->GetAxis(2)->SetTitle("sign");                         //  -1 | 0 | +1
 196   fHnDCA->GetAxis(3)->SetTitle("#eta_{Rec}");                   //  eta  [-0.9, 0.9]
 197   fHnDCA->GetAxis(4)->SetTitle("#it{y}_{Rec}");                 //  rapidity  [-0.5, 0.5]
 198   fHnDCA->GetAxis(5)->SetTitle("#varphi_{Rec} (rad)");          //  phi  [ 0. , 2Pi]
 199   fHnDCA->GetAxis(6)->SetTitle("#it{p}_{T,Rec} (GeV/#it{c})");  //  pT   [ 0.2, 2.6]
 200   fHnDCA->GetAxis(7)->SetTitle("contPart");                     //  1  primary | 2 missId | 3 from WeakDecay | 4 p from Material
 201   fHnDCA->GetAxis(8)->SetTitle("DCArAccepted");                 //  0 not accepted | 1 accepted
 202   fHnDCA->GetAxis(9)->SetTitle("DCAr");                         //  DCAr [-3, 3]
 203
 204
 205   fHelper->BinLogAxis(fHnDCA,  6, fESDTrackCuts);
 206   fHelper->BinLogAxis(fHnDCA,  6, fESDTrackCutsBkg);
 207
 208   // -- Set binning for DCAr
 209   Double_t binsDCAr[77] = {-3.,-2.85,-2.7,-2.55,-2.4,-2.25,-2.1,-1.95,-1.8,-1.65,-1.5,-1.35,-1.2,-1.05,-0.9,-0.75,-0.6,-0.45,-0.3,-0.285,-0.27,-0.255,-0.24,-0.225,-0.21,-0.195,-0.18,-0.165,-0.15,-0.135,-0.12,-0.105,-0.09,-0.075,-0.06,-0.045,-0.03,-0.015,0.,0.015,0.03,0.045,0.06,0.075,0.09,0.105,0.12,0.135,0.15,0.165,0.18,0.195,0.21,0.225,0.24,0.255,0.27,0.285,0.3,0.45,0.6,0.75,0.9,1.05,1.2,1.35,1.5,1.65,1.8,1.95,2.1,2.25,2.4,2.55,2.7,2.85,3.};
 210   fHnDCA->GetAxis(9)->Set(76, binsDCAr);
 211
 212   // ------------------------------------------------------------------
 213
 214   return;
 215 }
 216
 217 //________________________________________________________________________
 218 Int_t AliEbyEPidRatioDCA::GetContIdxTrack(Int_t label, Int_t sign, Int_t gPdgCode) {
 219   Int_t contIdx = -1;
 220
 221   AliVParticle* particle = (fESD) ? fMCEvent->GetTrack(label) : static_cast<AliVParticle*>(fArrayMC->At(label));
 222   if (!particle)
 223     return contIdx;
 224
 225   Bool_t isPhysicalPrimary        = (fESD) ? fStack->IsPhysicalPrimary(label): (static_cast<AliAODMCParticle*>(particle))->IsPhysicalPrimary();
 226   Bool_t isSecondaryFromWeakDecay = (fESD) ? fStack->IsSecondaryFromWeakDecay(label) : (static_cast<AliAODMCParticle*>(particle))->IsSecondaryFromWeakDecay();
 227   Bool_t isSecondaryFromMaterial  = (fESD) ? fStack->IsSecondaryFromMaterial(label)  : (static_cast<AliAODMCParticle*>(particle))->IsSecondaryFromMaterial();
 228
 229   if (isPhysicalPrimary) {
 230     if (gPdgCode == 0) {
 231       // -- Check if correctly identified
 232       if (particle->PdgCode() == (sign*gPdgCode))
 233         contIdx = 1;
 234       // -- MissIdentification
 235       else
 236         contIdx = 2;
 237     }
 238     else
 239       contIdx = 1;
 240   }
 241
 242   // -- Check if secondaries from material or weak decay
 243   else if(isSecondaryFromWeakDecay)
 244     contIdx = 3;
 245   else if (isSecondaryFromMaterial)
 246     contIdx = 4;
 247   else
 248     contIdx = -1;
 249
 250   return contIdx;
 251 }
 252
 253
 254