PWGCF/EBYE/NetParticle/AliAnalysisNetParticleDCA.cxx

   1 //-*- Mode: C++ -*-
   2
   3 #include "TMath.h"
   4 #include "TAxis.h"
   5
   6 #include "AliESDEvent.h"
   7 #include "AliESDInputHandler.h"
   8 #include "AliStack.h"
   9 #include "AliMCEvent.h"
  10
  11 #include "AliESDtrackCuts.h"
  12
  13 #include "AliAnalysisNetParticleDCA.h"
  14
  15 using namespace std;
  16
  17 // Task for NetParticle checks
  18 // Author: Jochen Thaeder <jochen@thaeder.de>
  19
  20 ClassImp(AliAnalysisNetParticleDCA)
  21
  22 /*
  23  * ---------------------------------------------------------------------------------
  24  *                            Constructor / Destructor
  25  * ---------------------------------------------------------------------------------
  26  */
  27
  28 //________________________________________________________________________
  29 AliAnalysisNetParticleDCA::AliAnalysisNetParticleDCA() :
  30   fHelper(NULL),
  31
  32   fESD(NULL),
  33   fESDTrackCuts(NULL),
  34   fESDTrackCutsBkg(NULL),
  35
  36   fStack(NULL),
  37   fMCEvent(NULL),
  38
  39   fHnDCA(NULL) {
  40   // Constructor
  41
  42   AliLog::SetClassDebugLevel("AliAnalysisNetParticleDCA",10);
  43 }
  44
  45 //________________________________________________________________________
  46 AliAnalysisNetParticleDCA::~AliAnalysisNetParticleDCA() {
  47   // Destructor
  48
  49 }
  50
  51
  52
  53 /*
  54  * ---------------------------------------------------------------------------------
  55  *                                 Public Methods
  56  * ---------------------------------------------------------------------------------
  57  */
  58
  59 //________________________________________________________________________
  60 void AliAnalysisNetParticleDCA::Initialize(AliESDtrackCuts *cuts, AliESDtrackCuts *cutsBkg, AliAnalysisNetParticleHelper* helper) {
  61
  62   // -- Get Helper
  63   // ---------------
  64   fHelper          = helper;
  65
  66   // -- ESD track cuts
  67   // -------------------
  68   fESDTrackCuts    = cuts;
  69   fESDTrackCutsBkg = cutsBkg;
  70
  71 #if 0
  72   // -- Get particle species / pdgCode
  73   // -------------------------
  74   fPdgCode          = AliPID::ParticleCode(fHelper->GetParticleSpecies());
  75
  76   // -- MC Labels for efficiency
  77   // -----------------------------
  78   fLabelsRec        = new Int_t*[2];
  79   for (Int_t ii = 0; ii < 2 ; ++ii)
  80     fLabelsRec[ii] = NULL;
  81
  82   // -- Create THnSparse Histograms
  83   // --------------------------------
  84   CreateHistograms();
  85 #endif
  86   return;
  87 }
  88
  89 /*
  90  * ---------------------------------------------------------------------------------
  91  *                            Setup/Reset Methods - private
  92  * ---------------------------------------------------------------------------------
  93  */
  94
  95 //________________________________________________________________________
  96 Int_t AliAnalysisNetParticleDCA::SetupEvent(AliESDInputHandler *esdHandler, AliMCEvent *mcEvent) {
  97   // -- Setup event
  98
  99   // -- Reset of event
 100   // -------------------
 101   ResetEvent();
 102
 103   // -- Setup of event
 104   // -------------------
 105   fESD           = esdHandler->GetEvent();
 106
 107   fMCEvent       = mcEvent;
 108   fStack         = fMCEvent->Stack();
 109 #if 0
 110   fCentralityBin = centBin;
 111
 112   // -- Create label arrays
 113   // ------------------------
 114   fLabelsRec[0] = new Int_t[fNTracks];
 115   if(!fLabelsRec[0]) {
 116     AliError("Cannot create fLabelsRec[0]");
 117     return -1;
 118   }
 119
 120   fLabelsRec[1] = new Int_t[fNTracks];
 121   if(!fLabelsRec[1]) {
 122     AliError("Cannot create fLabelsRec[1] for TPC");
 123     return -1;
 124   }
 125
 126   for(Int_t ii = 0; ii < fNTracks; ++ii) {
 127     fLabelsRec[0][ii] = 0;
 128     fLabelsRec[1][ii] = 0;
 129   }
 130 #endif
 131   return 0;
 132 }
 133
 134 //________________________________________________________________________
 135 void AliAnalysisNetParticleDCA::ResetEvent() {
 136   // -- Reset event
 137   /*
 138   for (Int_t ii = 0; ii < 2 ; ++ii) {
 139     if (fLabelsRec[ii])
 140       delete[] fLabelsRec[ii];
 141     fLabelsRec[ii] = NULL;
 142   }
 143   */
 144   return;
 145 }
 146
 147 //________________________________________________________________________
 148 void AliAnalysisNetParticleDCA::Process() {
 149   // -- Process event
 150
 151   // -- Setup (clean, create and fill) MC labels
 152   // ---------------------------------------------
 153   //  FillMCLabels();
 154
 155   // -- Fill  MC histograms for efficiency studies
 156   // -----------------------------------------------
 157   //  FillMCEffHist();
 158
 159   return;
 160 }
 161
 162 /*
 163  * ---------------------------------------------------------------------------------
 164  *                                 Private Methods
 165  * ---------------------------------------------------------------------------------
 166  */
 167 #if 0
 168 //________________________________________________________________________
 169 void AliAnalysisNetParticleDCA::CreateHistograms() {
 170   // -- Create histograms
 171
 172   Double_t dCent[2] = {-0.5, 8.5};
 173   Int_t iCent       = 9;
 174
 175   Double_t dEta[2]  = {-0.9, 0.9};
 176   Int_t iEta        = Int_t((dEta[1]-dEta[0]) / 0.075) +1 ;
 177
 178   Double_t dRap[2]  = {-0.5, 0.5};
 179   Int_t iRap        = Int_t((dRap[1]-dRap[0]) / 0.075) +1 ;
 180
 181   Double_t dPhi[2]  = {0.0, TMath::TwoPi()};
 182   Int_t iPhi        = 42;
 183
 184   Double_t dPt[2]   = {0.1, 3.0};
 185   Int_t iPt         = Int_t((dPt[1]-dPt[0]) / 0.075);
 186
 187   Double_t dSign[2] = {-1.5, 1.5};
 188   Int_t iSign       = 3;
 189
 190   // ------------------------------------------------------------------
 191   // -- Create THnSparseF - DCA
 192   // ------------------------------------------------------------------
 193
 194   //                              cent:   etaMC:     yMC:   phiMC:   ptMC:     sign:contPart:DCArAccepted:DCAr
 195   Int_t    binHnDCA[9] = {   iCent,    iEta,    iRap,    iPhi,    iPt,    iSign,       4,           2,  50 };
 196   Double_t minHnDCA[9] = {dCent[0], dEta[0], dRap[0], dPhi[0], dPt[0], dSign[0],     0.5,        -0.5, -10.};
 197   Double_t maxHnDCA[9] = {dCent[1], dEta[1], dRap[1], dPhi[1], dPt[1], dSign[1],     4.5,         1.5,  10.};
 198   fHnDCA = new THnSparseF("fHnDCA", "cent:etaMC:yMC:phiMC:ptMC:sign:contPart:contSign:DCArAccepted:DCAr", 9, binHnDCA, minHnDCA, maxHnDCA);
 199
 200   fHnDCA->Sumw2();
 201   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
 202   fHnDCA->GetAxis(1)->SetTitle("#eta_{MC}");                    //  eta  [-0.9,0.9]
 203   fHnDCA->GetAxis(2)->SetTitle("#it{y}_{MC}");                  //  rapidity  [-0.5, 0.5]
 204   fHnDCA->GetAxis(3)->SetTitle("#varphi_{MC} (rad)");           //  phi  [ 0. ,2Pi]
 205   fHnDCA->GetAxis(4)->SetTitle("#it{p}_{T,MC} (GeV/#it{c})");   //  pT   [ 0.1,1.3]
 206   fHnDCA->GetAxis(5)->SetTitle("sign");                         //  -1 | 0 | +1
 207   fHnDCA->GetAxis(6)->SetTitle("contPart");                     //  1  primary | 2 missId | 3 from WeakDecay | 4 p from Material
 208   fHnDCA->GetAxis(7)->SetTitle("DCArAccepted");                 //  0 not accepted | 1 accepted
 209   fHnDCA->GetAxis(8)->SetTitle("DCAr");                         //  DCAr [-10,10]
 210
 211   // ------------------------------------------------------------------
 212
 213   return;
 214 }
 215
 216 //________________________________________________________________________
 217 void AliAnalysisNetParticleDCA::FillMCDCA() {
 218   // Fill MC labels
 219   // Loop over ESD tracks and fill arrays with MC lables
 220   //  fLabelsRec[0] : all Tracks
 221   //  fLabelsRec[1] : all Tracks accepted by PID of TPC
 222   // Check every accepted track if correctly identified
 223   //  otherwise check for contamination
 224
 225   for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
 226     AliESDtrack *track = fESD->GetTrack(idxTrack);
 227
 228     // -- Check if track is accepted for basic parameters
 229     if (!fHelper->IsTrackAcceptedBasicCharged(track))
 230       continue;
 231
 232     // -- Check if accepted
 233     if (!fESDTrackCuts->AcceptTrack(track))
 234       continue;
 235
 236     // -- Check if accepted in rapidity window
 237     Double_t yP;
 238     if (!fHelper->IsTrackAcceptedRapidity(track, yP))
 239       continue;
 240
 241     // -- Check if accepted with thighter DCA cuts
 242     if (!fHelper->IsTrackAcceptedDCA(track))
 243       continue;
 244
 245     Int_t label  = TMath::Abs(track->GetLabel());
 246
 247     // -- Fill Label of all reconstructed
 248     fLabelsRec[0][idxTrack] = label;
 249
 250     // -- Check if accepted by PID from TPC or TPC+TOF
 251     Double_t pid[2];
 252     if (!fHelper->IsTrackAcceptedPID(track, pid))
 253       continue;
 254
 255     // -- Fill Label of all reconstructed && recPid_TPC+TOF
 256     fLabelsRec[1][idxTrack] = label;
 257
 258     // -- Check for contamination and fill contamination THnSparse
 259     CheckContTrack(label, track->GetSign(), idxTrack);
 260
 261   } // for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
 262
 263   return;
 264 }
 265
 266 //________________________________________________________________________
 267 void AliAnalysisNetParticleDCA::CheckDCATrack(Int_t label, Float_t sign, Float_t isDCArAccepted, Float_t dcar) {
 268   // Check if DCA of contamination or correctly identified
 269   // Fill contamination DCA THnSparse
 270
 271   TParticle* particle = fStack->Particle(label);
 272   if (!particle)
 273     return;
 274
 275   Bool_t isSecondaryFromWeakDecay = kFALSE;
 276   Bool_t isSecondaryFromMaterial  = kFALSE;
 277
 278   Int_t contPart = 0;
 279
 280   // -- Check if correctly identified
 281   if (particle->GetPdgCode() == (sign*fPdgCode)) {
 282
 283     // Check if is physical primary -> all ok
 284     if (fStack->IsPhysicalPrimary(label))
 285       contPart = 1;
 286     // -- Check if secondaries from weak decay
 287     else if(isSecondaryFromWeakDecay = fStack->IsSecondaryFromWeakDecay(label))
 288       contPart = 3;
 289     // -- Check if secondaries from material decay
 290     else if (isSecondaryFromMaterial  = fStack->IsSecondaryFromMaterial(label))
 291       contPart = 4;
 292   }
 293   // -- MissIdentification
 294   else
 295     contPart = 2;
 296
 297   Double_t hnDCA[9] = {fCentralityBin, particle->Eta(), particle->Y(), particle->Phi(), particle->Pt(), sign, contPart, isDCArAccepted, dcar};
 298   fHnDCA->Fill(hnContDCA);
 299
 300   return;
 301 }
 302
 303
 304 #endif