| 5ea3b761 |
1 | /************************************************************************** |
| 2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * |
| 3 | * * |
| 4 | * Author: Satyajit Jena. * |
| 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 | // // |
| 19 | // Analysis Task for Net-Charge Higher Moment Analysis // |
| 20 | // Author: Satyajit Jena || Nirbhay K. Behera // |
| 21 | // sjena@cern.ch || nbehera@cern.ch // |
| 22 | // // |
| 23 | //=========================================================================// |
| 24 | |
| 25 | #include "TChain.h" |
| 26 | #include "TList.h" |
| 27 | #include "TFile.h" |
| 28 | #include "TTree.h" |
| 29 | #include "TH1D.h" |
| 30 | #include "TH2F.h" |
| 31 | #include "THnSparse.h" |
| 32 | #include "TCanvas.h" |
| 33 | #include "TParticle.h" |
| 34 | #include <TDatabasePDG.h> |
| 35 | #include <TParticlePDG.h> |
| 36 | |
| 37 | #include "AliAnalysisTaskSE.h" |
| 38 | #include "AliAnalysisManager.h" |
| 39 | |
| 40 | #include "AliESDVertex.h" |
| 41 | #include "AliESDEvent.h" |
| 42 | #include "AliESDInputHandler.h" |
| 43 | #include "AliAODEvent.h" |
| 44 | #include "AliAODTrack.h" |
| 45 | #include "AliAODInputHandler.h" |
| 46 | #include "AliESD.h" |
| 47 | #include "AliESDEvent.h" |
| 48 | #include "AliAODEvent.h" |
| 49 | #include "AliStack.h" |
| 50 | #include "AliESDtrackCuts.h" |
| 51 | #include "AliAODMCHeader.h" |
| 52 | #include "AliAODMCParticle.h" |
| 53 | #include "AliMCEventHandler.h" |
| 54 | #include "AliMCEvent.h" |
| 55 | #include "AliMCParticle.h" |
| 56 | #include "AliStack.h" |
| 57 | #include "AliGenHijingEventHeader.h" |
| 58 | #include "AliGenEventHeader.h" |
| 59 | #include "AliPID.h" |
| 60 | #include "AliAODPid.h" |
| 61 | #include "AliPIDResponse.h" |
| 62 | #include "AliAODpidUtil.h" |
| 63 | #include "AliPIDCombined.h" |
| 64 | |
| 65 | #include "AliEbyEHigherMomentsEffContTask.h" |
| 66 | |
| 67 | using std::cout; |
| 68 | using std::endl; |
| 69 | |
| 70 | |
| 71 | ClassImp(AliEbyEHigherMomentsEffContTask) |
| 72 | |
| 73 | |
| 74 | //----------------------------------------------------------------------- |
| 75 | AliEbyEHigherMomentsEffContTask::AliEbyEHigherMomentsEffContTask( const char *name ) |
| 76 | : AliAnalysisTaskSE( name ), |
| 77 | fListOfHistosQA(0), |
| 78 | fListOfHistos(0), |
| 79 | fAOD(0), |
| 80 | fArrayMC(0), |
| 81 | fPIDResponse(0), |
| 82 | fParticleSpecies(AliPID::kProton), |
| 83 | fAnalysisType("AOD"), |
| 84 | fCentralityEstimator("V0M"), |
| 85 | fCentrality(0), |
| 86 | fVxMax(3.), |
| 87 | fVyMax(3.), |
| 88 | fVzMax(10.), |
| 89 | fDCAxy(2.4), |
| 90 | fDCAz(3.), |
| 91 | fPtLowerLimit(0.3), |
| 92 | fPtHigherLimit(1.5), |
| 93 | fEtaLowerLimit(-0.8), |
| 94 | fEtaHigherLimit(0.8), |
| 95 | fRapidityCut(0.5), |
| 96 | fNSigmaCut(3.), |
| 97 | fTPCNClus(80), |
| 98 | fChi2perNDF(4.), |
| 99 | fAODtrackCutBit(1024), |
| 100 | fLabel(NULL), |
| 101 | fUsePid(kFALSE), |
| 102 | fCheckCont(kFALSE), |
| 103 | fEff(kFALSE), |
| 104 | fEventCounter(0), |
| 105 | fHistDCA(0), |
| 106 | fTPCSig(0), |
| 107 | fTPCSigA(0), |
| 108 | fTHnCentNplusNminusCh(0), |
| 109 | fTHnCentNplusNminus(0), |
| 110 | fTHnEff(0), |
| 111 | fTHnCont(0) |
| 112 | { |
| 113 | |
| 114 | for ( Int_t i = 0; i < 13; i++) { |
| 115 | fHistQA[i] = NULL; |
| 116 | } |
| 117 | |
| 118 | for ( Int_t i = 0; i < 5; i++ ){ |
| 119 | fTHnCentNplusNminusPid[i] = NULL; |
| 120 | } |
| 121 | |
| 122 | DefineOutput(1, TList::Class()); // Outpput.... |
| 123 | DefineOutput(2, TList::Class()); |
| 124 | |
| 125 | } |
| 126 | |
| 127 | AliEbyEHigherMomentsEffContTask::~AliEbyEHigherMomentsEffContTask() { |
| 128 | |
| 129 | if(fListOfHistosQA) delete fListOfHistosQA; |
| 130 | if(fListOfHistos) delete fListOfHistos; |
| 131 | if(fLabel[0] ) delete [] (fLabel[0]); |
| 132 | if(fLabel[1] ) delete [] (fLabel[1]); |
| 133 | |
| 134 | } |
| 135 | |
| 136 | //--------------------------------------------------------------------------------- |
| 137 | void AliEbyEHigherMomentsEffContTask::UserCreateOutputObjects() { |
| 138 | //For QA-Histograms |
| 139 | fListOfHistosQA = new TList(); |
| 140 | fListOfHistosQA->SetOwner(kTRUE); |
| 141 | |
| 142 | fListOfHistos = new TList(); |
| 143 | fListOfHistos->SetOwner(kTRUE); |
| 144 | |
| 145 | fEventCounter = new TH1D("fEventCounter","EventCounter", 10, 0.5,10.5); |
| 146 | fEventCounter->GetXaxis()->SetBinLabel(1,"Event Accesed"); |
| 147 | fEventCounter->GetXaxis()->SetBinLabel(2,"Within 0-90% centrality"); |
| 148 | fEventCounter->GetXaxis()->SetBinLabel(5,"Have a vertex"); |
| 149 | fEventCounter->GetXaxis()->SetBinLabel(6,"After vertex Cut"); |
| 150 | fEventCounter->GetXaxis()->SetBinLabel(7,"Event Analyzed"); |
| 151 | fEventCounter->GetXaxis()->SetBinLabel(8,"Event Analysis finished"); |
| 152 | fListOfHistosQA->Add(fEventCounter); |
| 153 | |
| 154 | fHistQA[0] = new TH1D("fHistQAvx", "Histo Vx After Cut", 400, -4., 4.); |
| 155 | fHistQA[1] = new TH1D("fHistQAvy", "Histo Vy After Cut", 400, -4., 4.); |
| 156 | fHistQA[2] = new TH1D("fHistQAvz", "Histo Vz After Cut", 500, -25., 25.); |
| 157 | fHistQA[3] = new TH1D("fHistQAvxA", "Histo Vx All", 500, -5., 5.); |
| 158 | fHistQA[4] = new TH1D("fHistQAvyA", "Histo Vy All", 500, -5., 5.); |
| 159 | fHistQA[5] = new TH1D("fHistQAvzA", "Histo Vz All", 500, -25., 25.); |
| 160 | fHistQA[6] = new TH1D("fHistQADcaXyC", "Histo DCAxy after cut", 500, -5., 5.); |
| 161 | fHistQA[7] = new TH1D("fHistQADcaZC", "Histo DCAz after cut", 500, -5., 5.); |
| 162 | fHistQA[8] = new TH1D("fHistQAPt","p_{T} distribution",600,0,6); |
| 163 | fHistQA[9] = new TH1D("fHistQAEta","#eta distribution",240,-1.2,1.2); |
| 164 | fHistQA[10] = new TH1D("fHistQAPhi","#phi distribution",340,0,6.8); |
| 165 | fHistQA[11] = new TH1D("fHistQANCls","Number of TPC cluster",200,0,200); |
| 166 | fHistQA[12] = new TH1D("fHistQAChi2","Chi2 per NDF",100,0,10); |
| 167 | |
| 168 | for(Int_t i = 0; i < 13; i++) |
| 169 | { |
| 170 | fListOfHistosQA->Add(fHistQA[i]); |
| 171 | } |
| 172 | |
| 173 | fHistDCA = new TH2D("fHistDCA","DCAxy Vs DCAz", 500, -5., 5., 500, -5., 5.); |
| 174 | fTPCSig = new TH2D("fTPCSig","TPC signal",200, 0.0, 10. ,1000,0.,1000); |
| 175 | fTPCSig->SetMarkerColor(kRed); |
| 176 | fTPCSigA = new TH2D("fTPCSigA","TPC signal all ",200, 0.0, 10. ,1000,0.,1000); |
| 177 | fListOfHistosQA->Add(fHistDCA); |
| 178 | fListOfHistosQA->Add(fTPCSig); |
| 179 | fListOfHistosQA->Add(fTPCSigA); |
| 180 | |
| 181 | const Int_t nDim = 3; |
| 182 | const Int_t nPid = 5; |
| 183 | |
| 184 | TString Species[nPid] = {"Electron","Muon","Pion","Kaon","Proton"}; |
| 185 | |
| 186 | Int_t fBins[nPid][nDim]; |
| 187 | Double_t fMin[nPid][nDim]; |
| 188 | Double_t fMax[nPid][nDim]; |
| 189 | |
| 190 | for( Int_t i = 0; i < nPid; i++ ){ |
| 191 | for( Int_t j = 0; j < nDim; j++ ){ |
| 192 | fBins[i][j] = 0; |
| 193 | fMin[i][j] = 0.; |
| 194 | fMax[i][j] = 0.; |
| 195 | } |
| 196 | } |
| 197 | |
| 198 | |
| 199 | // if(fAnalysisType == "AOD" || fAnalysisType == "MCAOD") |
| 200 | |
| 201 | //{ |
| 202 | TString hname1; |
| 203 | TString htitle1, axisTitle1, axisTitle2; |
| 204 | |
| 205 | |
| 206 | |
| 207 | if( fUsePid ){ |
| 208 | |
| 209 | Int_t gPid = (Int_t)fParticleSpecies; |
| 210 | |
| 211 | if( gPid > 1 && gPid < 5 ){ |
| 212 | //Pion---- |
| 213 | fBins[2][0] = 100; fBins[2][1] = 1000; fBins[2][2] = 1000; |
| 214 | fMin[2][0] = -0.5; fMin[2][1] = -0.5; fMin[2][2] = -0.5; |
| 215 | fMax[2][0] = 99.5; fMax[2][1] = 999.5; fMax[2][2] = 999.5; |
| 216 | //Kaon------ |
| 217 | fBins[3][0] = 100; fBins[3][1] = 600; fBins[3][2] = 600; |
| 218 | fMin[3][0] = -0.5; fMin[3][1] = -0.5; fMin[3][2] = -0.5; |
| 219 | fMax[3][0] = 99.5; fMax[3][1] = 599.5; fMax[3][2] = 599.5; |
| 220 | //Proton----- |
| 221 | fBins[4][0] = 100; fBins[4][1] = 400; fBins[4][2] = 400; |
| 222 | fMin[4][0] = -0.5; fMin[4][1] = -0.5; fMin[4][2] = -0.5; |
| 223 | fMax[4][0] = 99.5; fMax[4][1] = 399.5; fMax[4][2] = 399.5; |
| 224 | |
| 225 | hname1 = "fCentNplusNminusPid"; hname1 +=gPid; |
| 226 | htitle1 = Species[gPid]; htitle1 +=" And Neg-"; htitle1 +=Species[gPid]; |
| 227 | axisTitle1 = Species[gPid]; |
| 228 | axisTitle2 ="Neg-"; axisTitle2 += Species[gPid]; |
| 229 | |
| 230 | fTHnCentNplusNminusPid[gPid] = new THnSparseD(hname1.Data(),htitle1.Data(), nDim, fBins[gPid], fMin[gPid],fMax[gPid]); |
| 231 | |
| 232 | fTHnCentNplusNminusPid[gPid]->GetAxis(0)->SetTitle("Centrality"); |
| 233 | fTHnCentNplusNminusPid[gPid]->GetAxis(1)->SetTitle(axisTitle1.Data()); |
| 234 | fTHnCentNplusNminusPid[gPid]->GetAxis(1)->SetTitle(axisTitle2.Data()); |
| 235 | |
| 236 | fListOfHistos->Add(fTHnCentNplusNminusPid[gPid]); |
| 237 | |
| 238 | }//Pion, Koan and Proton------- |
| 239 | |
| 240 | else{ |
| 241 | |
| 242 | Int_t fBinsX[nDim] = {100, 1500, 1500}; |
| 243 | Double_t fMinX[nDim] = { -0.5, -0.5, -0.5 }; |
| 244 | Double_t fMaxX[nDim] = { 99.5, 1499.5, 1499.5}; |
| 245 | |
| 246 | fTHnCentNplusNminus = new THnSparseD("fTHnCentNplusNminus","Cent-NplusChrg-NminusChrg", nDim, fBinsX, fMinX, fMaxX); |
| 247 | fTHnCentNplusNminus->GetAxis(0)->SetTitle("Centrality"); |
| 248 | fTHnCentNplusNminus->GetAxis(1)->SetTitle("UnwantedPlus"); |
| 249 | fTHnCentNplusNminus->GetAxis(2)->SetTitle("UnwantedMinus"); |
| 250 | fListOfHistos->Add(fTHnCentNplusNminus); |
| 251 | |
| 252 | }//Unwanted particle ------- |
| 253 | |
| 254 | }//fUsePid------- |
| 255 | |
| 256 | else{ |
| 257 | |
| 258 | Int_t fBinsCh[nDim] = {100, 1500, 1500}; |
| 259 | Double_t fMinCh[nDim] = { -0.5, -0.5, -0.5 }; |
| 260 | Double_t fMaxCh[nDim] = { 99.5, 1499.5, 1499.5}; |
| 261 | |
| 262 | fTHnCentNplusNminusCh = new THnSparseD("fTHnCentNplusNminusCh","Cent-NplusChrg-NminusChrg", nDim, fBinsCh, fMinCh, fMaxCh); |
| 263 | fTHnCentNplusNminusCh->GetAxis(0)->SetTitle("Centrality"); |
| 264 | fTHnCentNplusNminusCh->GetAxis(1)->SetTitle("Nplus"); |
| 265 | fTHnCentNplusNminusCh->GetAxis(2)->SetTitle("Nminus"); |
| 266 | fListOfHistos->Add(fTHnCentNplusNminusCh); |
| 267 | }//All charge hadrons--------- |
| 268 | |
| 269 | |
| 270 | |
| 271 | if( fAnalysisType == "MCAOD" ){ |
| 272 | |
| 273 | Double_t dCent[2] = {-0.5, 9.5}; |
| 274 | Int_t iCent = 10; |
| 275 | |
| 276 | Double_t dEta[2] = {-0.9, 0.9}; |
| 277 | Int_t iEta = Int_t((dEta[1]-dEta[0]) / 0.075) +1 ; |
| 278 | |
| 279 | Double_t dRap[2] = {-0.5, 0.5}; |
| 280 | Int_t iRap = Int_t((dRap[1]-dRap[0]) / 0.075) +1 ; |
| 281 | |
| 282 | Double_t dPhi[2] = {0.0, TMath::TwoPi()}; |
| 283 | Int_t iPhi = 42; |
| 284 | |
| 285 | Double_t dPt[2] = {0.1, 3.0}; |
| 286 | Int_t iPt = Int_t((dPt[1]-dPt[0]) / 0.075); |
| 287 | |
| 288 | Double_t dSign[2] = {-1.5, 1.5}; |
| 289 | Int_t iSign = 3; |
| 290 | |
| 291 | if(fEff){ |
| 292 | |
| 293 | Int_t nBinEff[15] = { iCent, iEta, iRap, iPhi, iPt, iSign, 2, 2, 2, iEta, iPhi, iPt, iEta, 2*iPhi+1, 2*iPt+1}; |
| 294 | Double_t nMinEff[15] = { dCent[0], dEta[0], dRap[0], dPhi[0], dPt[0], dSign[0], -0.5, -0.5, -0.5, dEta[0], dPhi[0], dPt[0], dEta[0], -dPhi[1], -dPt[1] }; |
| 295 | Double_t nMaxEff[15] = { dCent[1], dEta[1], dRap[1], dPhi[1], dPt[1], dSign[1], 1.5, 1.5, 1.5, dEta[1], dPhi[1], dPt[1], dEta[1], dPhi[1], dPt[1] }; |
| 296 | |
| 297 | fTHnEff = new THnSparseF("fHnEff", "cent:etaMC:yMC:phiMC:ptMC:sign:findable:recStatus:pidStatus:etaRec:phiRec:ptRec:deltaEta:deltaPhi:deltaPt", 15, nBinEff, nMinEff, nMaxEff); |
| 298 | |
| 299 | fTHnEff->Sumw2(); |
| 300 | |
| 301 | fTHnEff->GetAxis(0)->SetTitle("centrality"); |
| 302 | fTHnEff->GetAxis(1)->SetTitle("#eta_{MC}"); |
| 303 | fTHnEff->GetAxis(2)->SetTitle("#it{y}_{MC}"); |
| 304 | fTHnEff->GetAxis(3)->SetTitle("#varphi_{MC} (rad)"); |
| 305 | fTHnEff->GetAxis(4)->SetTitle("p_{T,MC} (GeV/#it{c})"); |
| 306 | fTHnEff->GetAxis(5)->SetTitle("sign"); |
| 307 | fTHnEff->GetAxis(6)->SetTitle("findable"); |
| 308 | fTHnEff->GetAxis(7)->SetTitle("recStatus"); |
| 309 | fTHnEff->GetAxis(8)->SetTitle("recPid"); |
| 310 | fTHnEff->GetAxis(9)->SetTitle("#eta_{Rec}"); |
| 311 | fTHnEff->GetAxis(10)->SetTitle("#varphi_{Rec} (rad)"); |
| 312 | fTHnEff->GetAxis(11)->SetTitle("p_{T,Rec} (GeV/#it{c})"); |
| 313 | fTHnEff->GetAxis(12)->SetTitle("#eta_{MC}-#eta_{Rec}"); |
| 314 | fTHnEff->GetAxis(13)->SetTitle("#varphi_{MC}-#varphi_{Rec} (rad)"); |
| 315 | fTHnEff->GetAxis(14)->SetTitle("p_{T,MC}-p_{T,Rec} (GeV/#it{c})"); |
| 316 | |
| 317 | fListOfHistos->Add(fTHnEff); |
| 318 | } |
| 319 | |
| 320 | if(fCheckCont){ |
| 321 | |
| 322 | Int_t nBinCont[14] = { iCent, iEta, iRap, iPhi, iPt, iSign, 8, iSign, iEta, iPhi, iPt, iEta, 2*iPhi+1, 2*iPt+1 }; |
| 323 | Double_t nMinCont[14] = {dCent[0], dEta[0], dRap[0], dPhi[0], dPt[0], dSign[0], 0.5, dSign[0], dEta[0], dPhi[0], dPt[0], dEta[0], -dPhi[1], -dPt[1]}; |
| 324 | Double_t nMaxCont[14] = { dCent[1], dEta[1], dRap[1], dPhi[1], dPt[1], dSign[1], 8.5, dSign[1], dEta[1], dPhi[1], dPt[1], dEta[1], dPhi[1], dPt[1] }; |
| 325 | |
| 326 | fTHnCont = new THnSparseF("fHnCont", "cent:etaMC:yMC:phiMC:ptMC:sign:contPart:contSign:etaRec:phiRec:ptRec:deltaEta:deltaPhi:deltaPt",14, nBinCont, nMinCont, nMaxCont); |
| 327 | fTHnCont->Sumw2(); |
| 328 | |
| 329 | fTHnCont->GetAxis(0)->SetTitle("centrality"); |
| 330 | fTHnCont->GetAxis(1)->SetTitle("#eta_{MC}"); |
| 331 | fTHnCont->GetAxis(2)->SetTitle("#it{y}_{MC}"); |
| 332 | fTHnCont->GetAxis(3)->SetTitle("#varphi_{MC} (rad)"); |
| 333 | fTHnCont->GetAxis(4)->SetTitle("p_{T,MC} (GeV/#it{c})"); |
| 334 | fTHnCont->GetAxis(5)->SetTitle("sign"); |
| 335 | fTHnCont->GetAxis(6)->SetTitle("contPart"); |
| 336 | fTHnCont->GetAxis(7)->SetTitle("contSign"); |
| 337 | fTHnCont->GetAxis(8)->SetTitle("#eta_{Rec}"); |
| 338 | fTHnCont->GetAxis(9)->SetTitle("#varphi_{Rec} (rad)"); |
| 339 | fTHnCont->GetAxis(10)->SetTitle("p_{T,Rec} (GeV/#it{c})"); |
| 340 | fTHnCont->GetAxis(11)->SetTitle("#eta_{MC}-#eta_{Rec}"); |
| 341 | fTHnCont->GetAxis(12)->SetTitle("#varphi_{MC}-#varphi_{Rec} (rad)"); |
| 342 | fTHnCont->GetAxis(13)->SetTitle("p_{T,MC}-p_{T,Rec} (GeV/#it{c})"); |
| 343 | |
| 344 | fListOfHistos->Add(fTHnCont); |
| 345 | } |
| 346 | |
| 347 | }//MCAOD------ |
| 348 | |
| 349 | PostData(1, fListOfHistosQA); |
| 350 | PostData(2, fListOfHistos); |
| 351 | |
| 352 | } |
| 353 | |
| 354 | |
| 355 | //---------------------------------------------------------------------------------- |
| 356 | void AliEbyEHigherMomentsEffContTask::UserExec( Option_t * ){ |
| 357 | |
| 358 | fEventCounter->Fill(1); |
| 359 | |
| 360 | |
| 361 | if(fAnalysisType == "AOD") { |
| 362 | |
| 363 | doAODEvent(); |
| 364 | |
| 365 | }//AOD--analysis----- |
| 366 | |
| 367 | else if(fAnalysisType == "MCAOD") { |
| 368 | |
| 369 | doMCAODEvent(); |
| 370 | |
| 371 | } |
| 372 | |
| 373 | else return; |
| 374 | |
| 375 | |
| 376 | |
| 377 | |
| 378 | fEventCounter->Fill(8); |
| 379 | |
| 380 | PostData(1, fListOfHistosQA); |
| 381 | PostData(2, fListOfHistos); |
| 382 | |
| 383 | } |
| 384 | |
| 385 | //-------------------------------------------------------------------------------------- |
| 386 | void AliEbyEHigherMomentsEffContTask::doAODEvent(){ |
| 387 | |
| 388 | //------------------- |
| 389 | Double_t nPlusCharge = 0.; |
| 390 | Double_t nMinusCharge = 0.; |
| 391 | Double_t nPartile = 0.; |
| 392 | Double_t nAntiParticle = 0.; |
| 393 | Int_t gPid = 0.; |
| 394 | |
| 395 | |
| 396 | //connect to the analysis mannager----- |
| 397 | |
| 398 | AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager(); |
| 399 | if (!manager) { |
| 400 | cout<<"ERROR: Analysis manager not found."<<endl; |
| 401 | return; |
| 402 | } |
| 403 | //coneect to the inputHandler------------ |
| 404 | AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler()); |
| 405 | if (!inputHandler) { |
| 406 | cout<<"ERROR: Input handler not found."<<endl; |
| 407 | return; |
| 408 | } |
| 409 | |
| 410 | fAOD = dynamic_cast<AliAODEvent*>(InputEvent()); |
| 411 | if (!fAOD) |
| 412 | { |
| 413 | cout<< "ERROR 01: AOD not found " <<endl; |
| 414 | return; |
| 415 | } |
| 416 | |
| 417 | fPIDResponse =inputHandler->GetPIDResponse(); |
| 418 | |
| 419 | |
| 420 | if (!fPIDResponse){ |
| 421 | AliFatal("This Task needs the PID response attached to the inputHandler"); |
| 422 | return; |
| 423 | } |
| 424 | |
| 425 | |
| 426 | AliAODHeader *aodHeader = fAOD->GetHeader(); |
| 427 | |
| 428 | Int_t cent = -1; |
| 429 | cent = aodHeader->GetCentralityP()->GetCentralityClass10(fCentralityEstimator.Data()); |
| 430 | |
| 431 | if (cent == 0) |
| 432 | fCentrality = aodHeader->GetCentralityP()->GetCentralityClass5(fCentralityEstimator.Data()); |
| 433 | else if (cent == 10 || cent == -1.) |
| 434 | fCentrality = -1; |
| 435 | else if (cent > 0 && cent < 9) |
| 436 | fCentrality = cent + 1; |
| 437 | |
| 438 | if(fCentrality < 0 || fCentrality >= 10) return; |
| 439 | |
| 440 | fEventCounter->Fill(2); |
| 441 | |
| 442 | if(!ProperVertex()) return; |
| 443 | |
| 444 | Int_t nTracks = fAOD->GetNumberOfTracks(); |
| 445 | |
| 446 | TExMap *trackMap = new TExMap();//Mapping matrix---- |
| 447 | |
| 448 | for(Int_t i = 0; i < nTracks; i++) { |
| 449 | |
| 450 | AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(i)); |
| 451 | |
| 452 | if(!aodTrack) { |
| 453 | AliError(Form("ERROR: Could not retrieve AODtrack %d",i)); |
| 454 | continue; |
| 455 | } |
| 456 | |
| 457 | Double_t tpcSignalAll = aodTrack->GetTPCsignal(); |
| 458 | fTPCSigA->Fill(aodTrack->GetTPCmomentum(),tpcSignalAll); |
| 459 | |
| 460 | Int_t gID = aodTrack->GetID(); |
| 461 | |
| 462 | if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, i);//Global tracks----- |
| 463 | }//1st track loop---- |
| 464 | |
| 465 | AliAODTrack* newAodTrack; |
| 466 | |
| 467 | for( Int_t j = 0; j < nTracks; j++ ) |
| 468 | { |
| 469 | |
| 470 | AliAODTrack* aodTrack1 = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(j)); |
| 471 | |
| 472 | if(!aodTrack1) { |
| 473 | AliError(Form("ERROR: Could not retrieve AODtrack %d",j)); |
| 474 | continue; |
| 475 | } |
| 476 | |
| 477 | |
| 478 | if(!aodTrack1->TestFilterBit(fAODtrackCutBit)) continue; |
| 479 | |
| 480 | Int_t gID = aodTrack1->GetID(); |
| 481 | |
| 482 | //if( aodTrack->GetID() != (-aodTrack1->GetID() -1) ) continue; |
| 483 | newAodTrack = gID >= 0 ? aodTrack1 : fAOD->GetTrack(trackMap->GetValue(-1-gID)); //Take those global track who corresponds to TPC only track |
| 484 | |
| 485 | Float_t dxy = 0., dz = 0.; |
| 486 | |
| 487 | dxy = aodTrack1->DCA(); |
| 488 | dz = aodTrack1->ZAtDCA(); |
| 489 | |
| 490 | Double_t pt = aodTrack1->Pt(); |
| 491 | Double_t eta = aodTrack1->Eta(); |
| 492 | Double_t nclus = aodTrack1->GetTPCClusterInfo(2,1); |
| 493 | Double_t chi2ndf = aodTrack1->Chi2perNDF(); |
| 494 | |
| 495 | /* |
| 496 | if( fabs(dxy) > fDCAxy ) continue; |
| 497 | if( fabs(dz) > fDCAz ) continue; |
| 498 | //Extra cut on DCA---( Similar to BF Task.. ) |
| 499 | if( fDCAxy !=0 && fDCAz !=0 ){ |
| 500 | if( TMath::Sqrt((dxy*dxy)/(fDCAxy*fDCAxy)+(dz*dz)/(fDCAz*fDCAz)) > 1. ) continue; |
| 501 | } |
| 502 | */ |
| 503 | if( pt < fPtLowerLimit || pt > fPtHigherLimit ) continue; |
| 504 | if( eta < fEtaLowerLimit || eta > fEtaHigherLimit ) continue; |
| 505 | if( nclus < fTPCNClus ) continue; |
| 506 | if( chi2ndf > fChi2perNDF ) continue; |
| 507 | |
| 508 | |
| 509 | fHistQA[6]->Fill(dxy); |
| 510 | fHistQA[7]->Fill(dz); |
| 511 | fHistQA[8]->Fill(pt); |
| 512 | fHistQA[9]->Fill(eta); |
| 513 | fHistQA[10]->Fill(aodTrack1->Phi()); |
| 514 | fHistQA[11]->Fill(nclus); |
| 515 | fHistQA[12]->Fill(chi2ndf); |
| 516 | fHistDCA->Fill(dxy,dz); |
| 517 | |
| 518 | Short_t gCharge = aodTrack1->Charge(); |
| 519 | |
| 520 | if(gCharge > 0) nPlusCharge += 1.; |
| 521 | if(gCharge < 0) nMinusCharge += 1.; |
| 522 | |
| 523 | if( fUsePid ) { |
| 524 | |
| 525 | gPid = (Int_t)fParticleSpecies; |
| 526 | |
| 527 | Double_t rap = newAodTrack->Y(AliPID::ParticleMass(fParticleSpecies)); |
| 528 | Double_t tpcSignal = newAodTrack->GetTPCsignal(); |
| 529 | //Double_t rap = aodTrack1->Y(AliPID::ParticleMass(fParticleSpecies)); |
| 530 | //Double_t tpcSignal = aodTrack1->GetTPCsignal(); |
| 531 | |
| 532 | if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut |
| 533 | |
| 534 | fTPCSig->Fill(newAodTrack->GetTPCmomentum(),tpcSignal); |
| 535 | |
| 536 | Float_t nsigmaTPCPID = -999.; |
| 537 | Float_t nsigmaTOFPID = -999.; |
| 538 | //Float_t nsigmaTPCTOFPID = -999.; |
| 539 | |
| 540 | nsigmaTPCPID = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,fParticleSpecies)); |
| 541 | nsigmaTOFPID = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(newAodTrack,fParticleSpecies)); |
| 542 | |
| 543 | if ( nsigmaTPCPID < fNSigmaCut ){ |
| 544 | |
| 545 | if (gCharge > 0) nPartile +=1.; |
| 546 | if( gCharge < 0 ) nAntiParticle +=1.; |
| 547 | |
| 548 | } |
| 549 | }//fUsepid----- |
| 550 | |
| 551 | }//--------- Track Loop to select with filterbit |
| 552 | |
| 553 | |
| 554 | |
| 555 | Double_t fContainerCh[3] = { fCentrality, nPlusCharge, nMinusCharge}; |
| 556 | Double_t fContainerPid[3] = { fCentrality, nPartile, nAntiParticle}; |
| 557 | |
| 558 | if( fUsePid ){ |
| 559 | |
| 560 | fTHnCentNplusNminusPid[gPid]->Fill(fContainerPid); |
| 561 | |
| 562 | // cout << "nCentrality "<< fCentrality <<", nParticle="<< nPartile << ", nMinusParticle=" << nAntiParticle << endl; |
| 563 | } |
| 564 | else{ |
| 565 | |
| 566 | fTHnCentNplusNminusCh->Fill(fContainerCh); |
| 567 | } |
| 568 | |
| 569 | //cout << "nCentrality "<< fCentrality <<", nPlusCharge="<< nPlusCharge << ", nMinusCharge=" << nMinusCharge << endl; |
| 570 | |
| 571 | fEventCounter->Fill(7); |
| 572 | |
| 573 | return; |
| 574 | |
| 575 | } |
| 576 | //-------------------------------------------------------------------------------------- |
| 577 | |
| 578 | |
| 579 | //-------------------------------------------------------------------------------------- |
| 580 | void AliEbyEHigherMomentsEffContTask::doMCAODEvent(){ |
| 581 | |
| 582 | |
| 583 | //--------- |
| 584 | Double_t nPlusCharge = 0.; |
| 585 | Double_t nMinusCharge = 0.; |
| 586 | Double_t nPartile = 0.; |
| 587 | Double_t nAntiParticle = 0.; |
| 588 | Int_t gPid = 0; |
| 589 | |
| 590 | //Connect to Anlaysis manager------ |
| 591 | AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager(); |
| 592 | if (!manager) { |
| 593 | cout<<"ERROR: Analysis manager not found."<<endl; |
| 594 | return; |
| 595 | } |
| 596 | |
| 597 | AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler()); |
| 598 | if (!inputHandler) { |
| 599 | cout<<"ERROR: Input handler not found."<<endl; |
| 600 | return; |
| 601 | } |
| 602 | |
| 603 | //AOD event------ |
| 604 | fAOD = dynamic_cast<AliAODEvent*>(InputEvent()); |
| 605 | if (!fAOD) |
| 606 | { |
| 607 | cout<< "ERROR 01: AOD not found " <<endl; |
| 608 | return; |
| 609 | } |
| 610 | |
| 611 | fPIDResponse =inputHandler->GetPIDResponse(); |
| 612 | |
| 613 | |
| 614 | if (!fPIDResponse){ |
| 615 | AliFatal("This Task needs the PID response attached to the inputHandler"); |
| 616 | return; |
| 617 | } |
| 618 | // -- Get MC header |
| 619 | // ------------------------------------------------------------------ |
| 620 | |
| 621 | fArrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName()); |
| 622 | if (!fArrayMC) { |
| 623 | AliFatal("Error: MC particles branch not found!\n"); |
| 624 | } |
| 625 | |
| 626 | AliAODMCHeader *mcHdr=NULL; |
| 627 | mcHdr=(AliAODMCHeader*)fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()); |
| 628 | if(!mcHdr) { |
| 629 | printf("MC header branch not found!\n"); |
| 630 | return; |
| 631 | } |
| 632 | |
| 633 | AliAODHeader *aodHeader = fAOD->GetHeader(); |
| 634 | |
| 635 | Int_t cent = -1; |
| 636 | cent = aodHeader->GetCentralityP()->GetCentralityClass10(fCentralityEstimator.Data()); |
| 637 | |
| 638 | if (cent == 0) |
| 639 | fCentrality = aodHeader->GetCentralityP()->GetCentralityClass5(fCentralityEstimator.Data()); |
| 640 | else if (cent == 10 || cent == -1.) |
| 641 | fCentrality = -1; |
| 642 | else if (cent > 0 && cent < 9) |
| 643 | fCentrality = cent + 1; |
| 644 | |
| 645 | if( fCentrality < 0 || fCentrality >= 10) return; |
| 646 | |
| 647 | fEventCounter->Fill(2); |
| 648 | |
| 649 | |
| 650 | |
| 651 | if(!ProperVertex()) return; |
| 652 | |
| 653 | Int_t nTracks = fAOD->GetNumberOfTracks(); |
| 654 | |
| 655 | //Creat Lables --- |
| 656 | fLabel = new Int_t*[2]; |
| 657 | fLabel[0] = new Int_t[nTracks]; //All charged hadrons---- |
| 658 | fLabel[1] = new Int_t[nTracks]; //For Pid----------- |
| 659 | //Initialize labels---- |
| 660 | |
| 661 | if(!fLabel[0]){ |
| 662 | AliError("Can't Get fLabel[0] "); |
| 663 | return; |
| 664 | } |
| 665 | if(!fLabel[1]){ |
| 666 | AliError("Can't Get fLabel[1] "); |
| 667 | return; |
| 668 | } |
| 669 | |
| 670 | for(Int_t i=0; i < 2; i++){ |
| 671 | for(Int_t j=0; j < nTracks; j++){ |
| 672 | fLabel[i][j] = 0; |
| 673 | } |
| 674 | } |
| 675 | |
| 676 | |
| 677 | TExMap *trackMap = new TExMap();//Mapping matrix---- |
| 678 | |
| 679 | for(Int_t i = 0; i < nTracks; i++) { |
| 680 | |
| 681 | AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(i)); |
| 682 | |
| 683 | if(!aodTrack) { |
| 684 | AliError(Form("ERROR: Could not retrieve AODtrack %d",i)); |
| 685 | continue; |
| 686 | } |
| 687 | |
| 688 | Double_t tpcSignalAll = aodTrack->GetTPCsignal(); |
| 689 | fTPCSigA->Fill(aodTrack->GetTPCmomentum(),tpcSignalAll); |
| 690 | |
| 691 | Int_t gID = aodTrack->GetID(); |
| 692 | |
| 693 | if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, i);//Global tracks----- |
| 694 | |
| 695 | }//1st track loop---- |
| 696 | |
| 697 | AliAODTrack* newAodTrack; |
| 698 | |
| 699 | for( Int_t j = 0; j < nTracks; j++ ){ |
| 700 | |
| 701 | AliAODTrack* aodTrack1 = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(j)); |
| 702 | |
| 703 | if(!aodTrack1) { |
| 704 | AliError(Form("ERROR: Could not retrieve AODtrack %d",j)); |
| 705 | continue; |
| 706 | } |
| 707 | |
| 708 | |
| 709 | if(!aodTrack1->TestFilterBit(fAODtrackCutBit)) continue; |
| 710 | |
| 711 | |
| 712 | Int_t gID = aodTrack1->GetID(); |
| 713 | |
| 714 | //if( aodTrack->GetID() != (-aodTrack1->GetID() -1) ) continue; |
| 715 | |
| 716 | newAodTrack = gID >= 0 ? aodTrack1 : fAOD->GetTrack(trackMap->GetValue(-1-gID)); //Take those global track who corresponds to TPC only track |
| 717 | |
| 718 | Float_t dxy = 0., dz = 0.; |
| 719 | |
| 720 | dxy = aodTrack1->DCA(); |
| 721 | dz = aodTrack1->ZAtDCA(); |
| 722 | |
| 723 | //cout << dxy << endl; |
| 724 | Double_t pt = aodTrack1->Pt(); |
| 725 | Double_t eta = aodTrack1->Eta(); |
| 726 | Double_t nclus = aodTrack1->GetTPCClusterInfo(2,1); |
| 727 | Double_t chi2ndf = aodTrack1->Chi2perNDF(); |
| 728 | |
| 729 | /* |
| 730 | if( fabs(dxy) > fDCAxy ) continue; |
| 731 | if( fabs(dz) > fDCAz ) continue; |
| 732 | //Extra cut on DCA---( Similar to BF Task.. ) |
| 733 | if( fDCAxy !=0 && fDCAz !=0 ){ |
| 734 | if( TMath::Sqrt((dxy*dxy)/(fDCAxy*fDCAxy)+(dz*dz)/(fDCAz*fDCAz)) > 1. ) continue; |
| 735 | } |
| 736 | */ |
| 737 | |
| 738 | |
| 739 | if( pt < fPtLowerLimit || pt > fPtHigherLimit ) continue; |
| 740 | if( eta < fEtaLowerLimit || eta > fEtaHigherLimit ) continue; |
| 741 | // if( nclus < fTPCNClus ) continue; |
| 742 | //if( chi2ndf > fChi2perNDF ) continue; |
| 743 | |
| 744 | |
| 745 | fHistQA[6]->Fill(dxy); |
| 746 | fHistQA[7]->Fill(dz); |
| 747 | fHistQA[8]->Fill(pt); |
| 748 | fHistQA[9]->Fill(eta); |
| 749 | fHistQA[10]->Fill(aodTrack1->Phi()); |
| 750 | fHistQA[11]->Fill(nclus); |
| 751 | fHistQA[12]->Fill(chi2ndf); |
| 752 | fHistDCA->Fill(dxy,dz); |
| 753 | |
| 754 | Short_t gCharge = aodTrack1->Charge(); |
| 755 | |
| 756 | if( gCharge == 0 ) continue; |
| 757 | |
| 758 | Int_t label = TMath::Abs(aodTrack1->GetLabel()); |
| 759 | //fill the labels-------- |
| 760 | fLabel[0][j] = label; |
| 761 | |
| 762 | if(gCharge > 0) nPlusCharge += 1.; |
| 763 | if(gCharge < 0) nMinusCharge += 1.; |
| 764 | |
| 765 | if( fUsePid ) { |
| 766 | |
| 767 | gPid = (Int_t)fParticleSpecies;; |
| 768 | Double_t rap = newAodTrack->Y(AliPID::ParticleMass(fParticleSpecies)); |
| 769 | Double_t tpcSignal = newAodTrack->GetTPCsignal(); |
| 770 | |
| 771 | if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut |
| 772 | |
| 773 | fTPCSig->Fill(newAodTrack->GetTPCmomentum(),tpcSignal); |
| 774 | |
| 775 | Float_t nsigmaTPCPID = -999.; |
| 776 | Float_t nsigmaTOFPID = -999.; |
| 777 | //Float_t nsigmaTPCTOFPID = -999.; |
| 778 | |
| 779 | nsigmaTPCPID = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,fParticleSpecies)); |
| 780 | nsigmaTOFPID = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(newAodTrack,fParticleSpecies)); |
| 781 | |
| 782 | if ( nsigmaTPCPID < fNSigmaCut ){ |
| 783 | |
| 784 | //Fill the labels---- |
| 785 | fLabel[1][j] = label; |
| 786 | |
| 787 | if (gCharge > 0) nPartile +=1.; |
| 788 | if( gCharge < 0 ) nAntiParticle +=1.; |
| 789 | |
| 790 | } |
| 791 | }//fUsepid----- |
| 792 | //Check Contamination------------------- |
| 793 | |
| 794 | if( fCheckCont ){ |
| 795 | CheckContTrackAOD(label, gCharge, j); |
| 796 | } |
| 797 | |
| 798 | }//--------- Track Loop to select with filterbit |
| 799 | |
| 800 | if(fEff) FillEffSparse(); |
| 801 | |
| 802 | //cout << "Cent=" << fCentrality << " MC-PlusChrg=" << nPlusCharge << " MC-MinusChrg=" << nMinusCharge << endl; |
| 803 | |
| 804 | //cout << "nCentrality "<< fCentrality <<", nParticle="<< nPartile << ", nMinusParticle=" << nAntiParticle << endl; |
| 805 | |
| 806 | Double_t fContainerCh[3] = { fCentrality, nPlusCharge, nMinusCharge}; |
| 807 | Double_t fContainerPid[3] = { fCentrality, nPartile, nAntiParticle}; |
| 808 | |
| 809 | if( fUsePid ){ |
| 810 | |
| 811 | fTHnCentNplusNminusPid[gPid]->Fill(fContainerPid); |
| 812 | |
| 813 | } |
| 814 | else{ |
| 815 | |
| 816 | fTHnCentNplusNminusCh->Fill(fContainerCh); |
| 817 | } |
| 818 | |
| 819 | fEventCounter->Fill(7); |
| 820 | |
| 821 | return; |
| 822 | |
| 823 | } |
| 824 | //--------------------------------------------------------------------------------------- |
| 825 | |
| 826 | //--------------------------------------------------------------------------------------- |
| 827 | Bool_t AliEbyEHigherMomentsEffContTask::ProperVertex(){ |
| 828 | |
| 829 | Bool_t IsVtx = kFALSE; |
| 830 | |
| 831 | const AliAODVertex *vertex = fAOD->GetPrimaryVertex(); |
| 832 | |
| 833 | if(vertex) { |
| 834 | Double32_t fCov[6]; |
| 835 | vertex->GetCovarianceMatrix(fCov); |
| 836 | if(vertex->GetNContributors() > 0) { |
| 837 | if(fCov[5] != 0) { |
| 838 | |
| 839 | Double_t lvx = vertex->GetX(); |
| 840 | Double_t lvy = vertex->GetY(); |
| 841 | Double_t lvz = vertex->GetZ(); |
| 842 | |
| 843 | fEventCounter->Fill(5); |
| 844 | |
| 845 | fHistQA[3]->Fill(lvx);fHistQA[4]->Fill(lvy);fHistQA[5]->Fill(lvz); |
| 846 | |
| 847 | if(TMath::Abs(lvx) < fVxMax) { |
| 848 | if(TMath::Abs(lvy) < fVyMax) { |
| 849 | if(TMath::Abs(lvz) < fVzMax) { |
| 850 | |
| 851 | fEventCounter->Fill(6); |
| 852 | fHistQA[0]->Fill(lvx);fHistQA[1]->Fill(lvy);fHistQA[2]->Fill(lvz); |
| 853 | |
| 854 | IsVtx = kTRUE; |
| 855 | |
| 856 | }//Z-Vertex cut--- |
| 857 | }//Y-vertex cut-- |
| 858 | }//X-vertex cut--- |
| 859 | }//Covariance------ |
| 860 | }//Contributors check--- |
| 861 | }//If vertex----------- |
| 862 | |
| 863 | return IsVtx; |
| 864 | } |
| 865 | //--------------------------------------------------------------------------------------- |
| 866 | |
| 867 | //--------------------------------------------------------------------------------------- |
| 868 | void AliEbyEHigherMomentsEffContTask::FillEffSparse(){ |
| 869 | |
| 870 | //For Efficiency------------------------------------- |
| 871 | |
| 872 | Int_t nTracks = fAOD->GetNumberOfTracks(); |
| 873 | Int_t gPdgCode = AliPID::ParticleCode(fParticleSpecies); |
| 874 | |
| 875 | Int_t nMCTrack = fArrayMC->GetEntriesFast(); |
| 876 | |
| 877 | |
| 878 | for (Int_t iMC = 0; iMC < nMCTrack; iMC++) { |
| 879 | |
| 880 | AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC); |
| 881 | |
| 882 | if(!partMC){ |
| 883 | AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC)); |
| 884 | continue; |
| 885 | } |
| 886 | |
| 887 | /* TDatabasePDG *db = TDatabasePDG::Instance(); |
| 888 | TParticlePDG *part = 0x0; |
| 889 | |
| 890 | if (!(part = db->GetParticle(partMC->PdgCode()))) continue; |
| 891 | if ((part = db->GetParticle(partMC->PdgCode()))->Charge() == 0.0) continue; |
| 892 | */ |
| 893 | if(!(partMC->PdgCode())) continue; |
| 894 | if(partMC->Charge() == 0) continue; |
| 895 | //-(1) Check for primary---- |
| 896 | if(!partMC->IsPhysicalPrimary()) continue; |
| 897 | |
| 898 | //-(2) Basic track cuts-------- |
| 899 | |
| 900 | if (partMC->Eta() < fEtaLowerLimit || partMC->Eta() > fEtaHigherLimit) continue; |
| 901 | if (partMC->Pt() < fPtLowerLimit || partMC->Pt() > fPtHigherLimit) continue; |
| 902 | |
| 903 | Double_t rap = partMC->Y(); |
| 904 | |
| 905 | if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut |
| 906 | |
| 907 | if(TMath::Abs(partMC->GetPdgCode()) != gPdgCode) continue; |
| 908 | |
| 909 | Short_t gCharge = partMC->Charge(); |
| 910 | |
| 911 | Float_t sign = (gCharge < 0 ) ? -1. : 1.; |
| 912 | |
| 913 | Float_t findable = 1.;// Float_t(fHelper->IsParticleFindable(idxMC)); XXX |
| 914 | |
| 915 | // (3) Initialize rec. variables----- |
| 916 | Float_t recStatus = 0.; |
| 917 | Float_t recPid = 0.; |
| 918 | // (4) Initialize track parameters------ |
| 919 | Float_t etaRec = 0.; |
| 920 | Float_t phiRec = 0.; |
| 921 | Float_t ptRec = 0.; |
| 922 | |
| 923 | |
| 924 | |
| 925 | // (5) Loop over all Labels-------- |
| 926 | |
| 927 | |
| 928 | for( Int_t iRec = 0; iRec < nTracks; iRec++ ){ |
| 929 | |
| 930 | if( iMC != fLabel[0][iRec] ) continue; |
| 931 | else{ |
| 932 | recStatus = 1.; |
| 933 | if( iMC == fLabel[1][iRec] ) |
| 934 | recPid = 1.; |
| 935 | |
| 936 | AliAODTrack* aodTrack = NULL; |
| 937 | |
| 938 | if(fAOD) |
| 939 | |
| 940 | aodTrack = fAOD->GetTrack(iRec); |
| 941 | |
| 942 | if(aodTrack){ |
| 943 | |
| 944 | etaRec = aodTrack->Eta(); |
| 945 | phiRec = aodTrack->Phi(); |
| 946 | ptRec = aodTrack->Pt(); |
| 947 | |
| 948 | }//aodTrack |
| 949 | |
| 950 | break; |
| 951 | } |
| 952 | |
| 953 | |
| 954 | |
| 955 | }//all charged track rec. check--- |
| 956 | |
| 957 | |
| 958 | //Fill ThnSparse---- |
| 959 | |
| 960 | Double_t EffContainer[15] = {fCentrality, partMC->Eta(), partMC->Y(), partMC->Phi(), partMC->Pt(), sign,findable, recStatus, recPid, etaRec, phiRec, ptRec, partMC->Eta()-etaRec, partMC->Phi()-phiRec, partMC->Pt()-ptRec}; |
| 961 | |
| 962 | fTHnEff->Fill(EffContainer); |
| 963 | |
| 964 | }//iMC Track loop----- |
| 965 | |
| 966 | return; |
| 967 | } |
| 968 | //--------------------------------------------------------------------------------------- |
| 969 | |
| 970 | //--------------------------------------------------------------------------------------- |
| 971 | void AliEbyEHigherMomentsEffContTask::CheckContTrackAOD(Int_t label, Float_t sign, Int_t idxTrack){ |
| 972 | |
| 973 | AliAODMCParticle* particle = (AliAODMCParticle*)fArrayMC->At(label); |
| 974 | |
| 975 | if (!particle) |
| 976 | return; |
| 977 | |
| 978 | Int_t gPdgCode = AliPID::ParticleCode(fParticleSpecies); |
| 979 | Bool_t isSecondaryFromWeakDecay = kFALSE; |
| 980 | Bool_t isSecondaryFromMaterial = kFALSE; |
| 981 | |
| 982 | // -- Check if correctly identified |
| 983 | if (particle->GetPdgCode() == (sign*gPdgCode)) { |
| 984 | |
| 985 | // Check if is physical primary -> all ok |
| 986 | if(particle->IsPhysicalPrimary()) |
| 987 | return; |
| 988 | |
| 989 | // -- Check if secondaries from material or weak decay |
| 990 | isSecondaryFromWeakDecay = particle->IsSecondaryFromWeakDecay(); |
| 991 | isSecondaryFromMaterial = particle->IsSecondaryFromMaterial(); |
| 992 | |
| 993 | } |
| 994 | |
| 995 | // -- Get MC pdg |
| 996 | Float_t contSign = 0.; |
| 997 | |
| 998 | // TDatabasePDG *db = TDatabasePDG::Instance(); |
| 999 | //TParticlePDG *part = 0x0; |
| 1000 | /* |
| 1001 | if(( part = db->GetParticle(particle->PdgCode()))->Charge() == 0.) contSign = 0.; |
| 1002 | else if((part = db->GetParticle(particle->PdgCode()))->Charge() < 0.) contSign = -1.; |
| 1003 | else if((part = db->GetParticle(particle->PdgCode()))->Charge() > 0.) contSign = 1.; |
| 1004 | */ |
| 1005 | |
| 1006 | if(particle->Charge() == 0.) contSign = 0.; |
| 1007 | if(particle->Charge() < 0. ) contSign = -1; |
| 1008 | if(particle->Charge() > 0. ) contSign = 1.; |
| 1009 | |
| 1010 | // -- Get contaminating particle |
| 1011 | Float_t contPart = 0; |
| 1012 | if (isSecondaryFromWeakDecay) contPart = 7; // probeParticle from WeakDecay |
| 1013 | else if (isSecondaryFromMaterial) contPart = 8; // probeParticle from Material |
| 1014 | else { |
| 1015 | if (TMath::Abs(particle->GetPdgCode()) == 211) contPart = 1; // pion |
| 1016 | else if (TMath::Abs(particle->GetPdgCode()) == 321) contPart = 2; // kaon |
| 1017 | else if (TMath::Abs(particle->GetPdgCode()) == 2212) contPart = 3; // proton |
| 1018 | else if (TMath::Abs(particle->GetPdgCode()) == 11) contPart = 4; // electron |
| 1019 | else if (TMath::Abs(particle->GetPdgCode()) == 13) contPart = 5; // muon |
| 1020 | else contPart = 6; // other |
| 1021 | } |
| 1022 | |
| 1023 | // -- Get Reconstructed values |
| 1024 | Float_t etaRec = 0.; |
| 1025 | Float_t phiRec = 0.; |
| 1026 | Float_t ptRec = 0.; |
| 1027 | |
| 1028 | AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(idxTrack)); |
| 1029 | |
| 1030 | if (aodTrack) { |
| 1031 | // if no track present (which should not happen) |
| 1032 | // -> pt = 0. , which is not in the looked at range |
| 1033 | |
| 1034 | // -- Get Reconstructed eta/phi/pt |
| 1035 | etaRec = aodTrack->Eta(); |
| 1036 | phiRec = aodTrack->Phi(); |
| 1037 | ptRec = aodTrack->Pt(); |
| 1038 | } |
| 1039 | |
| 1040 | // -- Fill THnSparse |
| 1041 | Double_t ContContainer[14] = { fCentrality, particle->Eta(), particle->Y(), particle->Phi(), particle->Pt(), sign, contPart, contSign, etaRec, phiRec, ptRec, particle->Eta()-etaRec, particle->Phi()-phiRec, particle->Pt()-ptRec }; |
| 1042 | |
| 1043 | fTHnCont->Fill(ContContainer); |
| 1044 | |
| 1045 | return; |
| 1046 | |
| 1047 | } |
| 1048 | //------------------------------------------------------------------------ |
| 1049 | void AliEbyEHigherMomentsEffContTask::Terminate( Option_t * ){ |
| 1050 | |
| 1051 | Info("AliEbyEHigherMomentsEffContTask"," Task Successfully finished"); |
| 1052 | |
| 1053 | } |
| 1054 | |
| 1055 | //------------------------------------------------------------------------ |
git://git.uio.no / u / mrichter / AliRoot.git / blame