| c683985a |
1 | #include "TChain.h" |
| 2 | #include "TList.h" |
| 3 | #include "TCanvas.h" |
| 4 | #include "TLorentzVector.h" |
| 5 | #include "TGraphErrors.h" |
| 6 | #include "TH1F.h" |
| 7 | #include "TH2F.h" |
| 8 | #include "TH3F.h" |
| 9 | #include "TH2D.h" |
| 10 | #include "TH3D.h" |
| 11 | #include "TArrayF.h" |
| 12 | #include "TF1.h" |
| 13 | #include "TRandom.h" |
| 6152cc3c |
14 | #include "TROOT.h" |
| c683985a |
15 | |
| 16 | #include "AliAnalysisTaskSE.h" |
| 17 | #include "AliAnalysisManager.h" |
| 18 | |
| 19 | #include "AliESDVertex.h" |
| 20 | #include "AliESDEvent.h" |
| 21 | #include "AliESDInputHandler.h" |
| 22 | #include "AliAODEvent.h" |
| 23 | #include "AliAODTrack.h" |
| 24 | #include "AliAODInputHandler.h" |
| 25 | #include "AliAODMCParticle.h" |
| 26 | #include "AliCollisionGeometry.h" |
| 27 | #include "AliGenEventHeader.h" |
| 28 | #include "AliMCEventHandler.h" |
| 29 | #include "AliMCEvent.h" |
| 30 | #include "AliStack.h" |
| 31 | #include "AliESDtrackCuts.h" |
| 32 | #include "AliEventplane.h" |
| 33 | #include "AliTHn.h" |
| 34 | #include "AliLog.h" |
| 35 | #include "AliAnalysisUtils.h" |
| 36 | |
| 37 | #include "AliEventPoolManager.h" |
| 38 | |
| 39 | #include "AliPID.h" |
| 40 | #include "AliPIDResponse.h" |
| 41 | #include "AliPIDCombined.h" |
| 42 | |
| 43 | #include "AliAnalysisTaskBFPsi.h" |
| 44 | #include "AliBalancePsi.h" |
| 45 | #include "AliAnalysisTaskTriggeredBF.h" |
| 46 | |
| 47 | |
| 48 | // Analysis task for the BF vs Psi code |
| 49 | // Authors: Panos.Christakoglou@nikhef.nl |
| 50 | |
| 51 | using std::cout; |
| 52 | using std::endl; |
| 53 | |
| 54 | ClassImp(AliAnalysisTaskBFPsi) |
| 55 | |
| 56 | //________________________________________________________________________ |
| 57 | AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name) |
| 58 | : AliAnalysisTaskSE(name), |
| 59 | fDebugLevel(kFALSE), |
| 60 | fArrayMC(0), //+++++++++++++ |
| 61 | fBalance(0), |
| 62 | fRunShuffling(kFALSE), |
| 63 | fShuffledBalance(0), |
| 64 | fRunMixing(kFALSE), |
| 65 | fRunMixingEventPlane(kFALSE), |
| 66 | fMixingTracks(50000), |
| 67 | fMixedBalance(0), |
| 68 | fPoolMgr(0), |
| 69 | fList(0), |
| 70 | fListBF(0), |
| 71 | fListBFS(0), |
| 72 | fListBFM(0), |
| 73 | fHistListPIDQA(0), |
| 74 | fHistEventStats(0), |
| 75 | fHistCentStats(0), |
| 76 | fHistCentStatsUsed(0), |
| 77 | fHistTriggerStats(0), |
| 78 | fHistTrackStats(0), |
| 79 | fHistVx(0), |
| 80 | fHistVy(0), |
| 81 | fHistVz(0), |
| 82 | fHistTPCvsVZEROMultiplicity(0), |
| 83 | fHistVZEROSignal(0), |
| 84 | fHistEventPlane(0), |
| 85 | fHistClus(0), |
| 86 | fHistDCA(0), |
| 87 | fHistChi2(0), |
| 88 | fHistPt(0), |
| 89 | fHistEta(0), |
| 90 | fHistRapidity(0), |
| 91 | fHistPhi(0), |
| 92 | fHistEtaPhiPos(0), |
| 93 | fHistEtaPhiNeg(0), |
| 94 | fHistPhiBefore(0), |
| 95 | fHistPhiAfter(0), |
| 96 | fHistPhiPos(0), |
| 97 | fHistPhiNeg(0), |
| 98 | fHistV0M(0), |
| 99 | fHistRefTracks(0), |
| 100 | fHistdEdxVsPTPCbeforePID(NULL), |
| 101 | fHistBetavsPTOFbeforePID(NULL), |
| 102 | fHistProbTPCvsPtbeforePID(NULL), |
| 103 | fHistProbTOFvsPtbeforePID(NULL), |
| 104 | fHistProbTPCTOFvsPtbeforePID(NULL), |
| 105 | fHistNSigmaTPCvsPtbeforePID(NULL), |
| 106 | fHistNSigmaTOFvsPtbeforePID(NULL), |
| 107 | fHistBetaVsdEdXbeforePID(NULL), //+++++++ |
| 108 | fHistNSigmaTPCTOFvsPtbeforePID(NULL), //++++++ |
| 109 | fHistdEdxVsPTPCafterPID(NULL), |
| 110 | fHistBetavsPTOFafterPID(NULL), |
| 111 | fHistProbTPCvsPtafterPID(NULL), |
| 112 | fHistProbTOFvsPtafterPID(NULL), |
| 113 | fHistProbTPCTOFvsPtafterPID(NULL), |
| 114 | fHistNSigmaTPCvsPtafterPID(NULL), |
| 115 | fHistNSigmaTOFvsPtafterPID(NULL), |
| 116 | fHistBetaVsdEdXafterPID(NULL), //+++++++ |
| 117 | fHistNSigmaTPCTOFvsPtafterPID(NULL), //+++++++ |
| 118 | fHistdEdxVsPTPCbeforePIDelectron(NULL), //+++++++ |
| 119 | fHistNSigmaTPCvsPtbeforePIDelectron(NULL), //+++++++ |
| 120 | fHistdEdxVsPTPCafterPIDelectron(NULL), //+++++++ |
| 121 | fHistNSigmaTPCvsPtafterPIDelectron(NULL), //+++++++ |
| 122 | fCentralityArrayBinsForCorrections(kCENTRALITY), |
| 123 | fCentralityWeights(0x0), |
| 124 | fPIDResponse(0x0), |
| 125 | fPIDCombined(0x0), |
| 126 | fParticleOfInterest(kPion), |
| 127 | fPidDetectorConfig(kTPCTOF), |
| 128 | fUsePID(kFALSE), |
| 129 | fUsePIDnSigma(kTRUE), |
| 130 | fUsePIDPropabilities(kFALSE), |
| 131 | fPIDNSigma(3.), |
| 132 | fMinAcceptedPIDProbability(0.8), |
| 133 | fElectronRejection(kFALSE), |
| 134 | fElectronOnlyRejection(kFALSE), |
| 135 | fElectronRejectionNSigma(-1.), |
| 136 | fElectronRejectionMinPt(0.), |
| 137 | fElectronRejectionMaxPt(1000.), |
| 138 | fESDtrackCuts(0), |
| 139 | fCentralityEstimator("V0M"), |
| 140 | fUseCentrality(kFALSE), |
| 141 | fCentralityPercentileMin(0.), |
| 142 | fCentralityPercentileMax(5.), |
| 143 | fImpactParameterMin(0.), |
| 144 | fImpactParameterMax(20.), |
| 145 | fMultiplicityEstimator("V0A"), |
| 146 | fUseMultiplicity(kFALSE), |
| 147 | fNumberOfAcceptedTracksMin(0), |
| 148 | fNumberOfAcceptedTracksMax(10000), |
| 149 | fHistNumberOfAcceptedTracks(0), |
| 150 | fHistMultiplicity(0), |
| 151 | fUseOfflineTrigger(kFALSE), |
| 152 | fCheckFirstEventInChunk(kFALSE), |
| 153 | fCheckPileUp(kFALSE), |
| 154 | fCheckPrimaryFlagAOD(kFALSE), |
| 155 | fUseMCforKinematics(kFALSE), |
| 156 | fVxMax(0.3), |
| 157 | fVyMax(0.3), |
| 158 | fVzMax(10.), |
| 159 | fnAODtrackCutBit(128), |
| 160 | fPtMin(0.3), |
| 161 | fPtMax(1.5), |
| c683985a |
162 | fEtaMin(-0.8), |
| 69d0337c |
163 | fEtaMax(0.8), |
| c683985a |
164 | fPhiMin(0.), |
| 165 | fPhiMax(360.), |
| c683985a |
166 | fDCAxyCut(-1), |
| 167 | fDCAzCut(-1), |
| 168 | fTPCchi2Cut(-1), |
| 169 | fNClustersTPCCut(-1), |
| 170 | fAcceptanceParameterization(0), |
| 171 | fDifferentialV2(0), |
| 172 | fUseFlowAfterBurner(kFALSE), |
| 173 | fExcludeResonancesInMC(kFALSE), |
| 174 | fExcludeElectronsInMC(kFALSE), |
| 175 | fUseMCPdgCode(kFALSE), |
| 176 | fPDGCodeToBeAnalyzed(-1), |
| 177 | fEventClass("EventPlane"), |
| 178 | fCustomBinning(""), |
| 179 | fHistVZEROAGainEqualizationMap(0), |
| 180 | fHistVZEROCGainEqualizationMap(0), |
| 181 | fHistVZEROChannelGainEqualizationMap(0) { |
| 182 | // Constructor |
| 183 | // Define input and output slots here |
| 184 | // Input slot #0 works with a TChain |
| 185 | |
| 186 | //======================================================correction |
| 187 | for (Int_t i=0; i<kCENTRALITY; i++){ |
| 188 | fHistCorrectionPlus[i] = NULL; |
| 189 | fHistCorrectionMinus[i] = NULL; |
| 190 | fCentralityArrayForCorrections[i] = -1.; |
| 191 | } |
| 192 | //=====================================================correction |
| 193 | |
| 194 | DefineInput(0, TChain::Class()); |
| 195 | // Output slot #0 writes into a TH1 container |
| 196 | DefineOutput(1, TList::Class()); |
| 197 | DefineOutput(2, TList::Class()); |
| 198 | DefineOutput(3, TList::Class()); |
| 199 | DefineOutput(4, TList::Class()); |
| 200 | DefineOutput(5, TList::Class()); |
| 201 | } |
| 202 | |
| 203 | //________________________________________________________________________ |
| 204 | AliAnalysisTaskBFPsi::~AliAnalysisTaskBFPsi() { |
| 205 | |
| 206 | // delete fBalance; |
| 207 | // delete fShuffledBalance; |
| 208 | // delete fList; |
| 209 | // delete fListBF; |
| 210 | // delete fListBFS; |
| 211 | |
| 212 | // delete fHistEventStats; |
| 213 | // delete fHistTrackStats; |
| 214 | // delete fHistVx; |
| 215 | // delete fHistVy; |
| 216 | // delete fHistVz; |
| 217 | |
| 218 | // delete fHistClus; |
| 219 | // delete fHistDCA; |
| 220 | // delete fHistChi2; |
| 221 | // delete fHistPt; |
| 222 | // delete fHistEta; |
| 223 | // delete fHistPhi; |
| 224 | // delete fHistEtaPhiPos; |
| 225 | // delete fHistEtaPhiNeg; |
| 226 | // delete fHistV0M; |
| 227 | } |
| 228 | |
| 229 | //________________________________________________________________________ |
| 230 | void AliAnalysisTaskBFPsi::UserCreateOutputObjects() { |
| 231 | // Create histograms |
| 232 | // Called once |
| 233 | |
| 234 | // global switch disabling the reference |
| 235 | // (to avoid "Replacing existing TH1" if several wagons are created in train) |
| 236 | Bool_t oldStatus = TH1::AddDirectoryStatus(); |
| 237 | TH1::AddDirectory(kFALSE); |
| 238 | |
| 239 | if(!fBalance) { |
| 240 | fBalance = new AliBalancePsi(); |
| 241 | fBalance->SetAnalysisLevel("ESD"); |
| 242 | fBalance->SetEventClass(fEventClass); |
| 243 | //fBalance->SetNumberOfBins(-1,16); |
| 244 | //fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.); |
| 245 | } |
| 246 | if(fRunShuffling) { |
| 247 | if(!fShuffledBalance) { |
| 248 | fShuffledBalance = new AliBalancePsi(); |
| 249 | fShuffledBalance->SetAnalysisLevel("ESD"); |
| 250 | fShuffledBalance->SetEventClass(fEventClass); |
| 251 | //fShuffledBalance->SetNumberOfBins(-1,16); |
| 252 | //fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.); |
| 253 | } |
| 254 | } |
| 255 | if(fRunMixing) { |
| 256 | if(!fMixedBalance) { |
| 257 | fMixedBalance = new AliBalancePsi(); |
| 258 | fMixedBalance->SetAnalysisLevel("ESD"); |
| 259 | fMixedBalance->SetEventClass(fEventClass); |
| 260 | } |
| 261 | } |
| 262 | |
| 263 | //QA list |
| 264 | fList = new TList(); |
| 265 | fList->SetName("listQA"); |
| 266 | fList->SetOwner(); |
| 267 | |
| 268 | //Balance Function list |
| 269 | fListBF = new TList(); |
| 270 | fListBF->SetName("listBF"); |
| 271 | fListBF->SetOwner(); |
| 272 | |
| 273 | if(fRunShuffling) { |
| 274 | fListBFS = new TList(); |
| 275 | fListBFS->SetName("listBFShuffled"); |
| 276 | fListBFS->SetOwner(); |
| 277 | } |
| 278 | |
| 279 | if(fRunMixing) { |
| 280 | fListBFM = new TList(); |
| 281 | fListBFM->SetName("listTriggeredBFMixed"); |
| 282 | fListBFM->SetOwner(); |
| 283 | } |
| 284 | |
| 285 | //PID QA list |
| 286 | if(fUsePID || fElectronRejection) { |
| 287 | fHistListPIDQA = new TList(); |
| 288 | fHistListPIDQA->SetName("listQAPID"); |
| 289 | fHistListPIDQA->SetOwner(); |
| 290 | } |
| 291 | |
| 292 | //Event stats. |
| 293 | TString gCutName[7] = {"Total","Offline trigger", |
| 294 | "Vertex","Analyzed","sel. Centrality","Not1stEvInChunk","No Pile-Up"}; |
| 295 | fHistEventStats = new TH2F("fHistEventStats", |
| 296 | "Event statistics;;Centrality percentile;N_{events}", |
| 297 | 7,0.5,7.5,220,-5,105); |
| 298 | for(Int_t i = 1; i <= 7; i++) |
| 299 | fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data()); |
| 300 | fList->Add(fHistEventStats); |
| 301 | |
| 302 | TString gCentName[13] = {"V0M","V0A","V0C","FMD","TRK","TKL","CL0","CL1","ZNA","ZPA","V0MvsFMD","TKLvsV0M","ZEMvsZDC"}; |
| 303 | fHistCentStats = new TH2F("fHistCentStats", |
| 304 | "Centrality statistics;;Cent percentile", |
| 305 | 13,-0.5,12.5,220,-5,105); |
| 306 | for(Int_t i = 1; i <= 13; i++){ |
| 307 | fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data()); |
| 308 | //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data()); |
| 309 | } |
| 310 | fList->Add(fHistCentStats); |
| 311 | |
| 312 | fHistCentStatsUsed = new TH2F("fHistCentStatsUsed","Centrality statistics;;Cent percentile", 1,-0.5,0.5,220,-5,105); |
| 313 | fHistCentStatsUsed->GetXaxis()->SetBinLabel(1,fCentralityEstimator.Data()); |
| 314 | fList->Add(fHistCentStatsUsed); |
| 315 | |
| 316 | fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",1025,0,1025); |
| 317 | fList->Add(fHistTriggerStats); |
| 318 | |
| 319 | fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",16,0,16); |
| 320 | fList->Add(fHistTrackStats); |
| 321 | |
| 322 | fHistNumberOfAcceptedTracks = new TH2F("fHistNumberOfAcceptedTracks",";N_{acc.};Centrality percentile;Entries",4001,-0.5,4000.5,220,-5,105); |
| 323 | fList->Add(fHistNumberOfAcceptedTracks); |
| 324 | |
| 325 | fHistMultiplicity = new TH1F("fHistMultiplicity",";N_{ch.};Entries",30001,-0.5,30000.5); |
| 326 | fList->Add(fHistMultiplicity); |
| 327 | |
| 328 | // Vertex distributions |
| 329 | fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5); |
| 330 | fList->Add(fHistVx); |
| 331 | fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5); |
| 332 | fList->Add(fHistVy); |
| 333 | fHistVz = new TH2F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Centrality percentile;Entries",100,-20.,20.,220,-5,105); |
| 334 | fList->Add(fHistVz); |
| 335 | |
| 336 | //TPC vs VZERO multiplicity |
| 337 | fHistTPCvsVZEROMultiplicity = new TH2F("fHistTPCvsVZEROMultiplicity","VZERO vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5); |
| 338 | if(fMultiplicityEstimator == "V0A") |
| 339 | fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO-A multiplicity (a.u.)"); |
| 340 | else if(fMultiplicityEstimator == "V0C") |
| 341 | fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO-C multiplicity (a.u.)"); |
| 342 | else |
| 343 | fHistTPCvsVZEROMultiplicity->GetXaxis()->SetTitle("VZERO multiplicity (a.u.)"); |
| 344 | fList->Add(fHistTPCvsVZEROMultiplicity); |
| 345 | |
| 346 | fHistVZEROSignal = new TH2F("fHistVZEROSignal","VZERO signal vs VZERO channel;VZERO channel; Signal (a.u.)",64,0.5,64.5,3001,-0.5,30000.5); |
| 347 | fList->Add(fHistVZEROSignal); |
| 348 | |
| 349 | //Event plane |
| 350 | fHistEventPlane = new TH2F("fHistEventPlane",";#Psi_{2} [deg.];Centrality percentile;Counts",100,0,360.,220,-5,105); |
| 351 | fList->Add(fHistEventPlane); |
| 352 | |
| 353 | // QA histograms |
| 354 | fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200); |
| 355 | fList->Add(fHistClus); |
| 356 | fHistChi2 = new TH2F("fHistChi2","Chi2/NDF distribution;#chi^{2}/ndf;Centrality percentile",200,0,10,220,-5,105); |
| 357 | fList->Add(fHistChi2); |
| 358 | fHistDCA = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5); |
| 359 | fList->Add(fHistDCA); |
| 360 | fHistPt = new TH2F("fHistPt","p_{T} distribution;p_{T} (GeV/c);Centrality percentile",200,0,10,220,-5,105); |
| 361 | fList->Add(fHistPt); |
| 362 | fHistEta = new TH2F("fHistEta","#eta distribution;#eta;Centrality percentile",200,-2,2,220,-5,105); |
| 363 | fList->Add(fHistEta); |
| 364 | fHistRapidity = new TH2F("fHistRapidity","y distribution;y;Centrality percentile",200,-2,2,220,-5,105); |
| 365 | fList->Add(fHistRapidity); |
| 366 | fHistPhi = new TH2F("fHistPhi","#phi distribution;#phi (rad);Centrality percentile",200,0.0,2.*TMath::Pi(),220,-5,105); |
| 367 | fList->Add(fHistPhi); |
| 368 | fHistEtaPhiPos = new TH3F("fHistEtaPhiPos","#eta-#phi distribution (+);#eta;#phi (rad);Centrality percentile",80,-2,2,72,0.0,2.*TMath::Pi(),220,-5,105); |
| 369 | fList->Add(fHistEtaPhiPos); |
| 370 | fHistEtaPhiNeg = new TH3F("fHistEtaPhiNeg","#eta-#phi distribution (-);#eta;#phi (rad);Centrality percentile",80,-2,2,72,0.0,2.*TMath::Pi(),220,-5,105); |
| 371 | fList->Add(fHistEtaPhiNeg); |
| 372 | fHistPhiBefore = new TH2F("fHistPhiBefore","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105); |
| 373 | fList->Add(fHistPhiBefore); |
| 374 | fHistPhiAfter = new TH2F("fHistPhiAfter","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105); |
| 375 | fList->Add(fHistPhiAfter); |
| 376 | fHistPhiPos = new TH2F("fHistPhiPos","#phi distribution for positive particles;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105); |
| 377 | fList->Add(fHistPhiPos); |
| 378 | fHistPhiNeg = new TH2F("fHistPhiNeg","#phi distribution for negative particles;#phi;Centrality percentile",200,0.,2.*TMath::Pi(),220,-5,105); |
| 379 | fList->Add(fHistPhiNeg); |
| 380 | fHistV0M = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000); |
| 381 | fList->Add(fHistV0M); |
| 382 | TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"}; |
| 383 | fHistRefTracks = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000); |
| 384 | for(Int_t i = 1; i <= 6; i++) |
| 385 | fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data()); |
| 386 | fList->Add(fHistRefTracks); |
| 387 | |
| 388 | // Balance function histograms |
| 389 | // Initialize histograms if not done yet (including the custom binning) |
| 390 | if(!fBalance->GetHistNp()){ |
| 391 | AliInfo("Histograms not yet initialized! --> Will be done now"); |
| 392 | fBalance->SetCustomBinning(fCustomBinning); |
| 393 | fBalance->InitHistograms(); |
| 394 | } |
| 395 | |
| 396 | if(fRunShuffling) { |
| 397 | if(!fShuffledBalance->GetHistNp()) { |
| 398 | AliInfo("Histograms (shuffling) not yet initialized! --> Will be done now"); |
| 399 | fShuffledBalance->SetCustomBinning(fCustomBinning); |
| 400 | fShuffledBalance->InitHistograms(); |
| 401 | } |
| 402 | } |
| 403 | |
| 404 | if(fRunMixing) { |
| 405 | if(!fMixedBalance->GetHistNp()) { |
| 406 | AliInfo("Histograms (mixing) not yet initialized! --> Will be done now"); |
| 407 | fMixedBalance->SetCustomBinning(fCustomBinning); |
| 408 | fMixedBalance->InitHistograms(); |
| 409 | } |
| 410 | } |
| 411 | |
| 412 | // QA histograms for different cuts |
| 413 | fList->Add(fBalance->GetQAHistHBTbefore()); |
| 414 | fList->Add(fBalance->GetQAHistHBTafter()); |
| 415 | fList->Add(fBalance->GetQAHistConversionbefore()); |
| 416 | fList->Add(fBalance->GetQAHistConversionafter()); |
| 417 | fList->Add(fBalance->GetQAHistPsiMinusPhi()); |
| 418 | fList->Add(fBalance->GetQAHistResonancesBefore()); |
| 419 | fList->Add(fBalance->GetQAHistResonancesRho()); |
| 420 | fList->Add(fBalance->GetQAHistResonancesK0()); |
| 421 | fList->Add(fBalance->GetQAHistResonancesLambda()); |
| 422 | fList->Add(fBalance->GetQAHistQbefore()); |
| 423 | fList->Add(fBalance->GetQAHistQafter()); |
| 424 | |
| 425 | //for(Int_t a = 0; a < ANALYSIS_TYPES; a++){ |
| 426 | fListBF->Add(fBalance->GetHistNp()); |
| 427 | fListBF->Add(fBalance->GetHistNn()); |
| 428 | fListBF->Add(fBalance->GetHistNpn()); |
| 429 | fListBF->Add(fBalance->GetHistNnn()); |
| 430 | fListBF->Add(fBalance->GetHistNpp()); |
| 431 | fListBF->Add(fBalance->GetHistNnp()); |
| 432 | |
| 433 | if(fRunShuffling) { |
| 434 | fListBFS->Add(fShuffledBalance->GetHistNp()); |
| 435 | fListBFS->Add(fShuffledBalance->GetHistNn()); |
| 436 | fListBFS->Add(fShuffledBalance->GetHistNpn()); |
| 437 | fListBFS->Add(fShuffledBalance->GetHistNnn()); |
| 438 | fListBFS->Add(fShuffledBalance->GetHistNpp()); |
| 439 | fListBFS->Add(fShuffledBalance->GetHistNnp()); |
| 440 | } |
| 441 | |
| 442 | if(fRunMixing) { |
| 443 | fListBFM->Add(fMixedBalance->GetHistNp()); |
| 444 | fListBFM->Add(fMixedBalance->GetHistNn()); |
| 445 | fListBFM->Add(fMixedBalance->GetHistNpn()); |
| 446 | fListBFM->Add(fMixedBalance->GetHistNnn()); |
| 447 | fListBFM->Add(fMixedBalance->GetHistNpp()); |
| 448 | fListBFM->Add(fMixedBalance->GetHistNnp()); |
| 449 | } |
| 450 | //} |
| 451 | |
| 452 | |
| 453 | // Event Mixing |
| 454 | if(fRunMixing){ |
| 455 | Int_t trackDepth = fMixingTracks; |
| 456 | Int_t poolsize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager |
| 457 | |
| 458 | // centrality bins |
| 459 | Double_t* centbins = NULL; |
| 460 | Int_t nCentralityBins; |
| 461 | if(fBalance->IsUseVertexBinning()){ |
| 462 | centbins = fBalance->GetBinning(fBalance->GetBinningString(), "centralityVertex", nCentralityBins); |
| 463 | } |
| 464 | else{ |
| 465 | centbins = fBalance->GetBinning(fBalance->GetBinningString(), "centrality", nCentralityBins); |
| 466 | } |
| 467 | |
| 468 | // multiplicity bins |
| 469 | Double_t* multbins = NULL; |
| 470 | Int_t nMultiplicityBins; |
| 471 | multbins = fBalance->GetBinning(fBalance->GetBinningString(), "multiplicity", nMultiplicityBins); |
| 472 | |
| 473 | // Zvtx bins |
| 474 | Double_t* vtxbins = NULL; |
| 475 | Int_t nVertexBins; |
| 476 | if(fBalance->IsUseVertexBinning()){ |
| 477 | vtxbins = fBalance->GetBinning(fBalance->GetBinningString(), "vertexVertex", nVertexBins); |
| 478 | } |
| 479 | else{ |
| 480 | vtxbins = fBalance->GetBinning(fBalance->GetBinningString(), "vertex", nVertexBins); |
| 481 | } |
| 482 | |
| 483 | // Event plane angle (Psi) bins |
| 484 | Double_t* psibins = NULL; |
| 485 | Int_t nPsiBins; |
| 486 | psibins = fBalance->GetBinning(fBalance->GetBinningString(), "eventPlane", nPsiBins); |
| 487 | |
| b78ae8b4 |
488 | |
| c683985a |
489 | // run the event mixing also in bins of event plane (statistics!) |
| 490 | if(fRunMixingEventPlane){ |
| 491 | if(fEventClass=="Multiplicity"){ |
| 492 | if(multbins && vtxbins && psibins){ |
| 493 | fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins, nPsiBins, psibins); |
| 494 | } |
| 495 | } |
| 496 | else{ |
| 497 | if(centbins && vtxbins && psibins){ |
| 498 | fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins); |
| 499 | } |
| 500 | } |
| 501 | } |
| 502 | else{ |
| 503 | if(fEventClass=="Multiplicity"){ |
| 504 | if(multbins && vtxbins){ |
| 505 | fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins); |
| 506 | } |
| 507 | } |
| 508 | else{ |
| 509 | if(centbins && vtxbins){ |
| 510 | fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins); |
| 511 | } |
| 512 | } |
| 513 | } |
| b78ae8b4 |
514 | |
| 515 | if(centbins) delete [] centbins; |
| 516 | if(multbins) delete [] multbins; |
| 517 | if(vtxbins) delete [] vtxbins; |
| 518 | |
| c683985a |
519 | // check pool manager |
| 520 | if(!fPoolMgr){ |
| 521 | AliError("Event Mixing required, but Pool Manager not initialized..."); |
| 522 | return; |
| 523 | } |
| 524 | } |
| 525 | |
| 526 | if(fESDtrackCuts) fList->Add(fESDtrackCuts); |
| 527 | |
| 528 | //====================PID========================// |
| 529 | if(fUsePID) { |
| 530 | fPIDCombined = new AliPIDCombined(); |
| 531 | fPIDCombined->SetDefaultTPCPriors(); |
| 532 | |
| 533 | fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000); |
| 534 | fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID); |
| 535 | |
| 536 | fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2); |
| 537 | fHistListPIDQA->Add(fHistBetavsPTOFbeforePID); |
| 538 | |
| 539 | fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0); |
| 540 | fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID); |
| 541 | |
| 542 | fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); |
| 543 | fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID); |
| 544 | |
| 545 | fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); |
| 546 | fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID); |
| 547 | |
| 548 | fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500); |
| 549 | fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID); |
| 550 | |
| 551 | fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, 0, 500); |
| 552 | fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID); |
| 553 | |
| 554 | fHistBetaVsdEdXbeforePID = new TH2D ("BetaVsdEdXbefore","BetaVsdEdXbefore", 1000, 0., 1000, 1000, 0, 1.2); |
| 555 | fHistListPIDQA->Add(fHistBetaVsdEdXbeforePID); //++++++++ |
| 556 | |
| 557 | fHistNSigmaTPCTOFvsPtbeforePID = new TH2D ("NSigmaTPCTOFvsPtbefore","NSigmaTPCTOFvsPtbefore", 1000, -10., 10., 1000, 0, 500); |
| 558 | fHistListPIDQA->Add(fHistNSigmaTPCTOFvsPtbeforePID); //++++++++ |
| 559 | |
| 560 | fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000); |
| 561 | fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID); |
| 562 | |
| 563 | fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2); |
| 564 | fHistListPIDQA->Add(fHistBetavsPTOFafterPID); |
| 565 | |
| 566 | fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2); |
| 567 | fHistListPIDQA->Add(fHistProbTPCvsPtafterPID); |
| 568 | |
| 569 | fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000, -10, 10, 1000, 0, 2); |
| 570 | fHistListPIDQA->Add(fHistProbTOFvsPtafterPID); |
| 571 | |
| 572 | fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0); |
| 573 | fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID); |
| 574 | |
| 575 | fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500); |
| 576 | fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID); |
| 577 | |
| 578 | fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, 0, 500); |
| 579 | fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID); |
| 580 | |
| 581 | fHistBetaVsdEdXafterPID = new TH2D ("BetaVsdEdXafter","BetaVsdEdXafter", 1000, 0., 1000, 1000, 0, 1.2); |
| 582 | fHistListPIDQA->Add(fHistBetaVsdEdXafterPID); //++++++++ |
| 583 | |
| 584 | fHistNSigmaTPCTOFvsPtafterPID = new TH2D ("NSigmaTPCTOFvsPtafter","NSigmaTPCTOFvsPtafter", 1000, -10., 10., 1000, 0, 500); |
| 585 | fHistListPIDQA->Add(fHistNSigmaTPCTOFvsPtafterPID); //++++++++ |
| 586 | } |
| 587 | |
| 588 | // for electron rejection only TPC nsigma histograms |
| 589 | if(fElectronRejection) { |
| 590 | |
| 591 | fHistdEdxVsPTPCbeforePIDelectron = new TH2D ("dEdxVsPTPCbeforeelectron","dEdxVsPTPCbeforeelectron", 1000, -10.0, 10.0, 1000, 0, 1000); |
| 592 | fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePIDelectron); |
| 593 | |
| 594 | fHistNSigmaTPCvsPtbeforePIDelectron = new TH2D ("NSigmaTPCvsPtbeforeelectron","NSigmaTPCvsPtbeforeelectron", 1000, -10, 10, 1000, 0, 500); |
| 595 | fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePIDelectron); |
| 596 | |
| 597 | fHistdEdxVsPTPCafterPIDelectron = new TH2D ("dEdxVsPTPCafterelectron","dEdxVsPTPCafterelectron", 1000, -10, 10, 1000, 0, 1000); |
| 598 | fHistListPIDQA->Add(fHistdEdxVsPTPCafterPIDelectron); |
| 599 | |
| 600 | fHistNSigmaTPCvsPtafterPIDelectron = new TH2D ("NSigmaTPCvsPtafterelectron","NSigmaTPCvsPtafterelectron", 1000, -10, 10, 1000, 0, 500); |
| 601 | fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPIDelectron); |
| 602 | } |
| 603 | //====================PID========================// |
| 604 | |
| 605 | // Post output data. |
| 606 | PostData(1, fList); |
| 607 | PostData(2, fListBF); |
| 608 | if(fRunShuffling) PostData(3, fListBFS); |
| 609 | if(fRunMixing) PostData(4, fListBFM); |
| 610 | if(fUsePID || fElectronRejection) PostData(5, fHistListPIDQA); //PID |
| 611 | |
| 612 | AliInfo("Finished setting up the Output"); |
| 613 | |
| 614 | TH1::AddDirectory(oldStatus); |
| b78ae8b4 |
615 | |
| c683985a |
616 | } |
| 617 | |
| 618 | |
| 619 | //________________________________________________________________________ |
| 620 | void AliAnalysisTaskBFPsi::SetInputCorrection(TString filename, |
| 621 | Int_t nCentralityBins, |
| 622 | Double_t *centralityArrayForCorrections) { |
| 623 | //Open files that will be used for correction |
| 624 | fCentralityArrayBinsForCorrections = nCentralityBins; |
| 625 | for (Int_t i=0; i<nCentralityBins; i++) |
| 626 | fCentralityArrayForCorrections[i] = centralityArrayForCorrections[i]; |
| 627 | |
| 628 | // No file specified -> run without corrections |
| 629 | if(!filename.Contains(".root")) { |
| 630 | AliInfo(Form("No correction file specified (= %s) --> run without corrections",filename.Data())); |
| 631 | return; |
| 632 | } |
| 633 | |
| 634 | //Open the input file |
| 635 | TFile *f = TFile::Open(filename); |
| 636 | if(!f->IsOpen()) { |
| 637 | AliInfo(Form("File %s not found --> run without corrections",filename.Data())); |
| 638 | return; |
| 639 | } |
| 640 | |
| 641 | //TString listEffName = ""; |
| 642 | for (Int_t iCent = 0; iCent < fCentralityArrayBinsForCorrections-1; iCent++) { |
| 643 | //Printf("iCent %d:",iCent); |
| 644 | TString histoName = "fHistCorrectionPlus"; |
| 645 | histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1])); |
| 646 | fHistCorrectionPlus[iCent]= dynamic_cast<TH3F *>(f->Get(histoName.Data())); |
| 647 | if(!fHistCorrectionPlus[iCent]) { |
| 648 | AliError(Form("fHist %s not found!!!",histoName.Data())); |
| 649 | return; |
| 650 | } |
| 651 | |
| 652 | histoName = "fHistCorrectionMinus"; |
| 653 | histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1])); |
| 654 | fHistCorrectionMinus[iCent] = dynamic_cast<TH3F *>(f->Get(histoName.Data())); |
| 655 | if(!fHistCorrectionMinus[iCent]) { |
| 656 | AliError(Form("fHist %s not found!!!",histoName.Data())); |
| 657 | return; |
| 658 | } |
| 659 | }//loop over centralities: ONLY the PbPb case is covered |
| c683985a |
660 | } |
| 661 | |
| 662 | //________________________________________________________________________ |
| 663 | void AliAnalysisTaskBFPsi::UserExec(Option_t *) { |
| 664 | // Main loop |
| 665 | // Called for each event |
| 666 | |
| 667 | TString gAnalysisLevel = fBalance->GetAnalysisLevel(); |
| 668 | Int_t gNumberOfAcceptedTracks = 0; |
| 669 | Double_t lMultiplicityVar = -999.; //-1 |
| 670 | Double_t gReactionPlane = -1.; |
| 671 | Float_t bSign = 0.; |
| 672 | |
| 673 | // get the event (for generator level: MCEvent()) |
| 674 | AliVEvent* eventMain = NULL; |
| 675 | if(gAnalysisLevel == "MC") { |
| 676 | eventMain = dynamic_cast<AliVEvent*>(MCEvent()); |
| 677 | } |
| 678 | else{ |
| 679 | eventMain = dynamic_cast<AliVEvent*>(InputEvent()); |
| 680 | // for HBT like cuts need magnetic field sign |
| 681 | bSign = (eventMain->GetMagneticField() > 0) ? 1 : -1; |
| 682 | } |
| 683 | if(!eventMain) { |
| 684 | AliError("eventMain not available"); |
| 685 | return; |
| 686 | } |
| 687 | |
| 688 | // PID Response task active? |
| 689 | if(fUsePID || fElectronRejection) { |
| 690 | fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse(); |
| 691 | if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler"); |
| 692 | } |
| 693 | |
| 694 | // check event cuts and fill event histograms |
| 695 | if((lMultiplicityVar = IsEventAccepted(eventMain)) < 0){ |
| 696 | return; |
| 697 | } |
| 698 | // get the reaction plane |
| 6152cc3c |
699 | if(fEventClass != "Multiplicity" && gAnalysisLevel!="AODnano") { |
| c683985a |
700 | gReactionPlane = GetEventPlane(eventMain); |
| 701 | fHistEventPlane->Fill(gReactionPlane,lMultiplicityVar); |
| 702 | if(gReactionPlane < 0){ |
| 703 | return; |
| 704 | } |
| 705 | } |
| 706 | |
| 707 | // get the accepted tracks in main event |
| 708 | TObjArray *tracksMain = GetAcceptedTracks(eventMain,lMultiplicityVar,gReactionPlane); |
| 709 | gNumberOfAcceptedTracks = tracksMain->GetEntriesFast(); |
| 710 | |
| 711 | //multiplicity cut (used in pp) |
| 712 | fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks,lMultiplicityVar); |
| 713 | |
| 714 | // store charges of all accepted tracks,shuffle and reassign(two extra loops) |
| 715 | TObjArray* tracksShuffled = NULL; |
| 716 | if(fRunShuffling){ |
| 717 | tracksShuffled = GetShuffledTracks(tracksMain,lMultiplicityVar); |
| 718 | } |
| 719 | |
| 720 | // Event mixing |
| 721 | if (fRunMixing) |
| 722 | { |
| 723 | // 1. First get an event pool corresponding in mult (cent) and |
| 724 | // zvertex to the current event. Once initialized, the pool |
| 725 | // should contain nMix (reduced) events. This routine does not |
| 726 | // pre-scan the chain. The first several events of every chain |
| 727 | // will be skipped until the needed pools are filled to the |
| 728 | // specified depth. If the pool categories are not too rare, this |
| 729 | // should not be a problem. If they are rare, you could lose` |
| 730 | // statistics. |
| 731 | |
| 732 | // 2. Collect the whole pool's content of tracks into one TObjArray |
| 733 | // (bgTracks), which is effectively a single background super-event. |
| 734 | |
| 735 | // 3. The reduced and bgTracks arrays must both be passed into |
| 736 | // FillCorrelations(). Also nMix should be passed in, so a weight |
| 737 | // of 1./nMix can be applied. |
| 738 | |
| 739 | AliEventPool* pool = fPoolMgr->GetEventPool(lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane); |
| 740 | |
| 741 | if (!pool){ |
| 742 | AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane)); |
| 743 | } |
| 744 | else{ |
| 745 | |
| 746 | //pool->SetDebug(1); |
| 747 | |
| 748 | if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5){ |
| 749 | |
| 750 | |
| 751 | Int_t nMix = pool->GetCurrentNEvents(); |
| 752 | //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl; |
| 753 | |
| 754 | //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2); |
| 755 | //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool()); |
| 756 | //if (pool->IsReady()) |
| 757 | //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3); |
| 758 | |
| 759 | // Fill mixed-event histos here |
| 760 | for (Int_t jMix=0; jMix<nMix; jMix++) |
| 761 | { |
| 762 | TObjArray* tracksMixed = pool->GetEvent(jMix); |
| 763 | fMixedBalance->CalculateBalance(gReactionPlane,tracksMain,tracksMixed,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ()); |
| 764 | } |
| 765 | } |
| 766 | |
| 767 | // Update the Event pool |
| 768 | pool->UpdatePool(tracksMain); |
| 769 | //pool->PrintInfo(); |
| 770 | |
| 771 | }//pool NULL check |
| 772 | }//run mixing |
| 773 | |
| 774 | // calculate balance function |
| 775 | fBalance->CalculateBalance(gReactionPlane,tracksMain,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ()); |
| 776 | |
| 777 | // calculate shuffled balance function |
| 778 | if(fRunShuffling && tracksShuffled != NULL) { |
| 779 | fShuffledBalance->CalculateBalance(gReactionPlane,tracksShuffled,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ()); |
| 780 | } |
| 781 | } |
| 782 | |
| 783 | //________________________________________________________________________ |
| 784 | Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){ |
| 785 | // Checks the Event cuts |
| 786 | // Fills Event statistics histograms |
| 787 | |
| 788 | Bool_t isSelectedMain = kTRUE; |
| 789 | Float_t gRefMultiplicity = -1.; |
| 790 | TString gAnalysisLevel = fBalance->GetAnalysisLevel(); |
| 791 | |
| 792 | AliMCEvent *mcevent = dynamic_cast<AliMCEvent*>(event); |
| 793 | |
| 794 | fHistEventStats->Fill(1,gRefMultiplicity); //all events |
| 795 | |
| 796 | // check first event in chunk (is not needed for new reconstructions) |
| 797 | if(fCheckFirstEventInChunk){ |
| 798 | AliAnalysisUtils ut; |
| 799 | if(ut.IsFirstEventInChunk(event)) |
| 800 | return -1.; |
| 801 | fHistEventStats->Fill(6,gRefMultiplicity); |
| 802 | } |
| 803 | // check for pile-up event |
| 804 | if(fCheckPileUp){ |
| 805 | AliAnalysisUtils ut; |
| 806 | ut.SetUseMVPlpSelection(kTRUE); |
| 807 | ut.SetUseOutOfBunchPileUp(kTRUE); |
| 808 | if(ut.IsPileUpEvent(event)) |
| 809 | return -1.; |
| 810 | fHistEventStats->Fill(7,gRefMultiplicity); |
| 811 | } |
| 812 | |
| 813 | // Event trigger bits |
| 814 | fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()); |
| 815 | if(fUseOfflineTrigger) |
| 816 | isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected(); |
| 817 | |
| 818 | if(isSelectedMain) { |
| 819 | fHistEventStats->Fill(2,gRefMultiplicity); //triggered events |
| 820 | |
| 821 | // Event Vertex MC |
| 822 | if(gAnalysisLevel == "MC") { |
| 823 | if(!event) { |
| 824 | AliError("mcEvent not available"); |
| 825 | return 0x0; |
| 826 | } |
| 827 | |
| 828 | if(mcevent){ |
| 829 | AliGenEventHeader *header = dynamic_cast<AliGenEventHeader*>(mcevent->GenEventHeader()); |
| 830 | if(header){ |
| 831 | TArrayF gVertexArray; |
| 832 | header->PrimaryVertex(gVertexArray); |
| 833 | //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)", |
| 834 | //gVertexArray.At(0), |
| 835 | //gVertexArray.At(1), |
| 836 | //gVertexArray.At(2)); |
| 837 | fHistEventStats->Fill(3,gRefMultiplicity); //events with a proper vertex |
| 838 | if(TMath::Abs(gVertexArray.At(0)) < fVxMax) { |
| 839 | if(TMath::Abs(gVertexArray.At(1)) < fVyMax) { |
| 840 | if(TMath::Abs(gVertexArray.At(2)) < fVzMax) { |
| 841 | fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events |
| 842 | |
| 843 | // get the reference multiplicty or centrality |
| 844 | gRefMultiplicity = GetRefMultiOrCentrality(event); |
| 845 | |
| 846 | fHistVx->Fill(gVertexArray.At(0)); |
| 847 | fHistVy->Fill(gVertexArray.At(1)); |
| 848 | fHistVz->Fill(gVertexArray.At(2),gRefMultiplicity); |
| 849 | |
| 850 | // take only events inside centrality class |
| 851 | if(fUseCentrality) { |
| 852 | if((fImpactParameterMin < gRefMultiplicity) && (fImpactParameterMax > gRefMultiplicity)){ |
| 853 | fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality |
| 854 | return gRefMultiplicity; |
| 855 | }//centrality class |
| 856 | } |
| 857 | // take events only within the same multiplicity class |
| 858 | else if(fUseMultiplicity) { |
| 859 | if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) { |
| 860 | fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity |
| 861 | return gRefMultiplicity; |
| 862 | } |
| 863 | }//multiplicity range |
| 864 | }//Vz cut |
| 865 | }//Vy cut |
| 866 | }//Vx cut |
| 867 | }//header |
| 868 | }//MC event object |
| 869 | }//MC |
| 870 | |
| 871 | // Event Vertex AOD, ESD, ESDMC |
| 872 | else{ |
| 873 | const AliVVertex *vertex = event->GetPrimaryVertex(); |
| 874 | |
| 875 | if(vertex) { |
| 876 | Double32_t fCov[6]; |
| 877 | vertex->GetCovarianceMatrix(fCov); |
| 878 | if(vertex->GetNContributors() > 0) { |
| 879 | if(fCov[5] != 0) { |
| 880 | fHistEventStats->Fill(3,gRefMultiplicity); //proper vertex |
| 881 | if(TMath::Abs(vertex->GetX()) < fVxMax) { |
| 882 | if(TMath::Abs(vertex->GetY()) < fVyMax) { |
| 883 | if(TMath::Abs(vertex->GetZ()) < fVzMax) { |
| 884 | fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events |
| 885 | |
| 886 | // get the reference multiplicty or centrality |
| 887 | gRefMultiplicity = GetRefMultiOrCentrality(event); |
| 888 | |
| 889 | fHistVx->Fill(vertex->GetX()); |
| 890 | fHistVy->Fill(vertex->GetY()); |
| 891 | fHistVz->Fill(vertex->GetZ(),gRefMultiplicity); |
| 892 | |
| 893 | // take only events inside centrality class |
| 894 | if(fUseCentrality) { |
| 895 | if((gRefMultiplicity > fCentralityPercentileMin) && (gRefMultiplicity < fCentralityPercentileMax)){ |
| 896 | |
| 897 | // centrality weighting (optional for 2011 if central and semicentral triggers are used) |
| 898 | if (fCentralityWeights && !AcceptEventCentralityWeight(gRefMultiplicity)){ |
| 899 | AliInfo(Form("Rejecting event because of centrality weighting: %f", gRefMultiplicity)); |
| 900 | return -1; |
| 901 | } |
| 902 | |
| 903 | fHistEventStats->Fill(5,gRefMultiplicity); //events with correct centrality |
| 904 | return gRefMultiplicity; |
| 905 | }//centrality class |
| 906 | } |
| 907 | // take events only within the same multiplicity class |
| 908 | else if(fUseMultiplicity) { |
| 909 | //if(fDebugLevel) |
| 910 | //Printf("N(min): %.0f, N(max): %.0f - N(ref): %.0f",fNumberOfAcceptedTracksMin, |
| 911 | //fNumberOfAcceptedTracksMax,gRefMultiplicity); |
| 912 | |
| 913 | if((gRefMultiplicity > fNumberOfAcceptedTracksMin) && (gRefMultiplicity < fNumberOfAcceptedTracksMax)) { |
| 914 | fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity |
| 915 | return gRefMultiplicity; |
| 916 | } |
| 917 | }//multiplicity range |
| 918 | }//Vz cut |
| 919 | }//Vy cut |
| 920 | }//Vx cut |
| 921 | }//proper vertex resolution |
| 922 | }//proper number of contributors |
| 923 | }//vertex object valid |
| 924 | }//triggered event |
| 925 | }//AOD,ESD,ESDMC |
| 926 | |
| 927 | // in all other cases return -1 (event not accepted) |
| 928 | return -1; |
| 929 | } |
| 930 | |
| 931 | |
| 932 | //________________________________________________________________________ |
| 933 | Double_t AliAnalysisTaskBFPsi::GetRefMultiOrCentrality(AliVEvent *event){ |
| 934 | // Checks the Event cuts |
| 935 | // Fills Event statistics histograms |
| 6152cc3c |
936 | |
| c683985a |
937 | Float_t gCentrality = -1.; |
| 938 | Double_t gMultiplicity = -1.; |
| 939 | TString gAnalysisLevel = fBalance->GetAnalysisLevel(); |
| 940 | |
| 941 | |
| 942 | // calculate centrality always (not only in centrality mode) |
| 943 | if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ) { //centrality in AOD header //++++++++++++++ |
| 944 | AliAODHeader *header = (AliAODHeader*) event->GetHeader(); |
| 945 | if(header){ |
| 946 | gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data()); |
| 947 | |
| 948 | // QA for centrality estimators |
| 949 | fHistCentStats->Fill(0.,header->GetCentralityP()->GetCentralityPercentile("V0M")); |
| 950 | fHistCentStats->Fill(1.,header->GetCentralityP()->GetCentralityPercentile("V0A")); |
| 951 | fHistCentStats->Fill(2.,header->GetCentralityP()->GetCentralityPercentile("V0C")); |
| 952 | fHistCentStats->Fill(3.,header->GetCentralityP()->GetCentralityPercentile("FMD")); |
| 953 | fHistCentStats->Fill(4.,header->GetCentralityP()->GetCentralityPercentile("TRK")); |
| 954 | fHistCentStats->Fill(5.,header->GetCentralityP()->GetCentralityPercentile("TKL")); |
| 955 | fHistCentStats->Fill(6.,header->GetCentralityP()->GetCentralityPercentile("CL0")); |
| 956 | fHistCentStats->Fill(7.,header->GetCentralityP()->GetCentralityPercentile("CL1")); |
| 957 | fHistCentStats->Fill(8.,header->GetCentralityP()->GetCentralityPercentile("ZNA")); |
| 958 | fHistCentStats->Fill(9.,header->GetCentralityP()->GetCentralityPercentile("ZPA")); |
| 959 | fHistCentStats->Fill(10.,header->GetCentralityP()->GetCentralityPercentile("V0MvsFMD")); |
| 960 | fHistCentStats->Fill(11.,header->GetCentralityP()->GetCentralityPercentile("TKLvsV0M")); |
| 961 | fHistCentStats->Fill(12.,header->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC")); |
| 962 | |
| 963 | // Centrality estimator USED ++++++++++++++++++++++++++++++ |
| 964 | fHistCentStatsUsed->Fill(0.,header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data())); |
| 965 | |
| 966 | // centrality QA (V0M) |
| 967 | fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C()); |
| 968 | |
| 969 | // centrality QA (reference tracks) |
| 970 | fHistRefTracks->Fill(0.,header->GetRefMultiplicity()); |
| 971 | fHistRefTracks->Fill(1.,header->GetRefMultiplicityPos()); |
| 972 | fHistRefTracks->Fill(2.,header->GetRefMultiplicityNeg()); |
| 973 | fHistRefTracks->Fill(3.,header->GetTPConlyRefMultiplicity()); |
| 974 | fHistRefTracks->Fill(4.,header->GetNumberOfITSClusters(0)); |
| 975 | fHistRefTracks->Fill(5.,header->GetNumberOfITSClusters(1)); |
| 976 | fHistRefTracks->Fill(6.,header->GetNumberOfITSClusters(2)); |
| 977 | fHistRefTracks->Fill(7.,header->GetNumberOfITSClusters(3)); |
| 978 | fHistRefTracks->Fill(8.,header->GetNumberOfITSClusters(4)); |
| 979 | |
| 980 | }//AOD header |
| 981 | }//AOD |
| 6152cc3c |
982 | |
| 983 | // calculate centrality always (not only in centrality mode) |
| 984 | else if(gAnalysisLevel == "AODnano" ) { //centrality via JF workaround |
| 985 | |
| 986 | AliAODHeader *header = (AliAODHeader*) event->GetHeader(); |
| 987 | if(header){ |
| 988 | gCentrality = (Float_t) gROOT->ProcessLine(Form("100.0 + 100.0 * ((AliNanoAODHeader*) %p)->GetCentrality(\"%s\")", header,fCentralityEstimator.Data())) / 100 - 1.0; |
| 989 | |
| 990 | // QA histogram |
| 991 | fHistCentStatsUsed->Fill(0.,gCentrality); |
| 992 | |
| 993 | }//AOD header |
| 994 | }//AODnano |
| c683985a |
995 | |
| 996 | else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs |
| 997 | AliCentrality *centrality = event->GetCentrality(); |
| 998 | gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data()); |
| 999 | |
| 1000 | // QA for centrality estimators |
| 1001 | fHistCentStats->Fill(0.,centrality->GetCentralityPercentile("V0M")); |
| 1002 | fHistCentStats->Fill(1.,centrality->GetCentralityPercentile("V0A")); |
| 1003 | fHistCentStats->Fill(2.,centrality->GetCentralityPercentile("V0C")); |
| 1004 | fHistCentStats->Fill(3.,centrality->GetCentralityPercentile("FMD")); |
| 1005 | fHistCentStats->Fill(4.,centrality->GetCentralityPercentile("TRK")); |
| 1006 | fHistCentStats->Fill(5.,centrality->GetCentralityPercentile("TKL")); |
| 1007 | fHistCentStats->Fill(6.,centrality->GetCentralityPercentile("CL0")); |
| 1008 | fHistCentStats->Fill(7.,centrality->GetCentralityPercentile("CL1")); |
| 1009 | fHistCentStats->Fill(8.,centrality->GetCentralityPercentile("ZNA")); |
| 1010 | fHistCentStats->Fill(9.,centrality->GetCentralityPercentile("ZPA")); |
| 1011 | fHistCentStats->Fill(10.,centrality->GetCentralityPercentile("V0MvsFMD")); |
| 1012 | fHistCentStats->Fill(11.,centrality->GetCentralityPercentile("TKLvsV0M")); |
| 1013 | fHistCentStats->Fill(12.,centrality->GetCentralityPercentile("ZEMvsZDC")); |
| 1014 | |
| 1015 | // Centrality estimator USED ++++++++++++++++++++++++++++++ |
| 1016 | fHistCentStatsUsed->Fill(0.,centrality->GetCentralityPercentile(fCentralityEstimator.Data())); |
| 1017 | |
| 1018 | // centrality QA (V0M) |
| 1019 | fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C()); |
| 1020 | }//ESD |
| 1021 | |
| 1022 | else if(gAnalysisLevel == "MC"){ |
| 1023 | Double_t gImpactParameter = 0.; |
| 1024 | AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event); |
| 1025 | if(gMCEvent){ |
| 1026 | AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader()); |
| 1027 | if(headerH){ |
| 1028 | gImpactParameter = headerH->ImpactParameter(); |
| 1029 | gCentrality = gImpactParameter; |
| 1030 | }//MC header |
| 1031 | }//MC event cast |
| 1032 | }//MC |
| 1033 | |
| 1034 | else{ |
| 1035 | gCentrality = -1.; |
| 1036 | } |
| 1037 | |
| 1038 | // calculate reference multiplicity always (not only in multiplicity mode) |
| 1039 | if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ |
| 1040 | AliESDEvent* gESDEvent = dynamic_cast<AliESDEvent*>(event); |
| 1041 | if(gESDEvent){ |
| 1042 | gMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(gESDEvent, AliESDtrackCuts::kTrackletsITSTPC,0.5); |
| 1043 | fHistMultiplicity->Fill(gMultiplicity); |
| 1044 | }//AliESDevent cast |
| 1045 | }//ESD mode |
| 1046 | |
| 1047 | else if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ){ |
| 1048 | AliAODHeader *header = (AliAODHeader*) event->GetHeader(); |
| 1049 | if ((fMultiplicityEstimator == "V0M")|| |
| 1050 | (fMultiplicityEstimator == "V0A")|| |
| 1051 | (fMultiplicityEstimator == "V0C") || |
| 1052 | (fMultiplicityEstimator == "TPC")) { |
| 1053 | gMultiplicity = GetReferenceMultiplicityFromAOD(event); |
| 1054 | if(fDebugLevel) Printf("Reference multiplicity (calculated): %.0f",gMultiplicity); |
| 1055 | } |
| 1056 | else { |
| 1057 | if(header) |
| 1058 | gMultiplicity = header->GetRefMultiplicity(); |
| 1059 | if(fDebugLevel) Printf("Reference multiplicity (AOD header): %.0f",gMultiplicity); |
| 1060 | } |
| 1061 | fHistMultiplicity->Fill(gMultiplicity); |
| 1062 | }//AOD mode |
| 1063 | else if(gAnalysisLevel == "MC") { |
| 1064 | AliMCEvent* gMCEvent = dynamic_cast<AliMCEvent*>(event); |
| 1065 | //Calculating the multiplicity as the number of charged primaries |
| 1066 | //within \pm 0.8 in eta and pT > 0.1 GeV/c |
| 1067 | for(Int_t iParticle = 0; iParticle < gMCEvent->GetNumberOfPrimaries(); iParticle++) { |
| 1068 | AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iParticle)); |
| 1069 | if (!track) { |
| 1070 | AliError(Form("Could not receive particle %d", iParticle)); |
| 1071 | continue; |
| 1072 | } |
| 1073 | |
| 1074 | //exclude non stable particles |
| 1075 | if(!(gMCEvent->IsPhysicalPrimary(iParticle))) continue; |
| 1076 | |
| 1077 | //++++++++++++++++ |
| 1078 | if (fMultiplicityEstimator == "V0M") { |
| 1079 | if((track->Eta() > 5.1 || track->Eta() < 2.8)&&(track->Eta() < -3.7 || track->Eta() > -1.7)) |
| 1080 | continue;} |
| 1081 | else if (fMultiplicityEstimator == "V0A") { |
| 1082 | if(track->Eta() > 5.1 || track->Eta() < 2.8) continue;} |
| 1083 | else if (fMultiplicityEstimator == "V0C") { |
| 1084 | if(track->Eta() > -1.7 || track->Eta() < -3.7) continue;} |
| 1085 | else if (fMultiplicityEstimator == "TPC") { |
| 1086 | if(track->Eta() < fEtaMin || track->Eta() > fEtaMax) continue; |
| 1087 | if(track->Pt() < fPtMin || track->Pt() > fPtMax) continue; |
| 1088 | } |
| 1089 | else{ |
| 1090 | if(track->Pt() < fPtMin || track->Pt() > fPtMax) continue; |
| 1091 | if(track->Eta() < fEtaMin || track->Eta() > fEtaMax) continue; |
| 1092 | } |
| 1093 | //++++++++++++++++ |
| 1094 | |
| 1095 | if(track->Charge() == 0) continue; |
| 1096 | |
| 1097 | gMultiplicity += 1; |
| 1098 | }//loop over primaries |
| 1099 | fHistMultiplicity->Fill(gMultiplicity); |
| 1100 | }//MC mode |
| 1101 | else{ |
| 1102 | gMultiplicity = -1; |
| 1103 | } |
| 1104 | |
| 1105 | |
| 1106 | // decide what should be returned only here |
| 1107 | Double_t lReturnVal = -100; |
| 1108 | if(fEventClass=="Multiplicity"){ |
| 1109 | lReturnVal = gMultiplicity; |
| 1110 | }else if(fEventClass=="Centrality"){ |
| 1111 | lReturnVal = gCentrality; |
| 1112 | } |
| 1113 | return lReturnVal; |
| 1114 | } |
| 1115 | |
| 1116 | //________________________________________________________________________ |
| 1117 | Double_t AliAnalysisTaskBFPsi::GetReferenceMultiplicityFromAOD(AliVEvent *event){ |
| 1118 | //Function that returns the reference multiplicity from AODs (data or reco MC) |
| 1119 | //Different ref. mult. implemented: V0M, V0A, V0C, TPC |
| 1120 | Double_t gRefMultiplicity = 0., gRefMultiplicityTPC = 0.; |
| 1121 | Double_t gRefMultiplicityVZERO = 0., gRefMultiplicityVZEROA = 0., gRefMultiplicityVZEROC = 0.; |
| 1122 | |
| 1123 | AliAODHeader *header = dynamic_cast<AliAODHeader *>(event->GetHeader()); |
| 1124 | if(!header) { |
| 1125 | Printf("ERROR: AOD header not available"); |
| 1126 | return -999; |
| 1127 | } |
| 1128 | Int_t gRunNumber = header->GetRunNumber(); |
| 1129 | |
| 1130 | // Loop over tracks in event |
| 1131 | for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) { |
| 1132 | AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks)); |
| 1133 | if (!aodTrack) { |
| 1134 | AliError(Form("Could not receive track %d", iTracks)); |
| 1135 | continue; |
| 1136 | } |
| 1137 | |
| 1138 | // AOD track cuts |
| 1139 | if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue; |
| 1140 | |
| 1141 | if(aodTrack->Charge() == 0) continue; |
| 1142 | // Kinematics cuts from ESD track cuts |
| 1143 | if( aodTrack->Pt() < fPtMin || aodTrack->Pt() > fPtMax) continue; |
| 1144 | if( aodTrack->Eta() < fEtaMin || aodTrack->Eta() > fEtaMax) continue; |
| 1145 | |
| 1146 | //=================PID (so far only for electron rejection)==========================// |
| 1147 | if(fElectronRejection) { |
| 1148 | // get the electron nsigma |
| 1149 | Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron)); |
| 1150 | |
| 1151 | // check only for given momentum range |
| 1152 | if( aodTrack->Pt() > fElectronRejectionMinPt && aodTrack->Pt() < fElectronRejectionMaxPt ){ |
| 1153 | //look only at electron nsigma |
| 1154 | if(!fElectronOnlyRejection) { |
| 1155 | //Make the decision based on the n-sigma of electrons |
| 1156 | if(nSigma < fElectronRejectionNSigma) continue; |
| 1157 | } |
| 1158 | else { |
| 1159 | Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion)); |
| 1160 | Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon)); |
| 1161 | Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton)); |
| 1162 | |
| 1163 | //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species) |
| 1164 | if(nSigma < fElectronRejectionNSigma |
| 1165 | && nSigmaPions > fElectronRejectionNSigma |
| 1166 | && nSigmaKaons > fElectronRejectionNSigma |
| 1167 | && nSigmaProtons > fElectronRejectionNSigma ) continue; |
| 1168 | } |
| 1169 | } |
| 1170 | }//electron rejection |
| 1171 | |
| 1172 | gRefMultiplicityTPC += 1.0; |
| 1173 | }// track loop |
| 1174 | |
| 1175 | //VZERO segmentation in two detectors (0-31: VZERO-C, 32-63: VZERO-A) |
| 1176 | for(Int_t iChannel = 0; iChannel < 64; iChannel++) { |
| 1177 | fHistVZEROSignal->Fill(iChannel,event->GetVZEROEqMultiplicity(iChannel)); |
| 1178 | |
| 1179 | if(iChannel < 32) |
| 1180 | gRefMultiplicityVZEROC += event->GetVZEROEqMultiplicity(iChannel); |
| 1181 | else if(iChannel >= 32) |
| 1182 | gRefMultiplicityVZEROA += event->GetVZEROEqMultiplicity(iChannel); |
| 1183 | }//loop over PMTs |
| 1184 | |
| 1185 | //Equalization of gain |
| 1186 | Double_t gFactorA = GetEqualizationFactor(gRunNumber,"A"); |
| 1187 | if(gFactorA != 0) |
| 1188 | gRefMultiplicityVZEROA /= gFactorA; |
| 1189 | Double_t gFactorC = GetEqualizationFactor(gRunNumber,"C"); |
| 1190 | if(gFactorC != 0) |
| 1191 | gRefMultiplicityVZEROC /= gFactorC; |
| 1192 | if((gFactorA != 0)&&(gFactorC != 0)) |
| 1193 | gRefMultiplicityVZERO = (gRefMultiplicityVZEROA/gFactorA)+(gRefMultiplicityVZEROC/gFactorC); |
| 1194 | |
| 1195 | if(fDebugLevel) |
| 1196 | Printf("VZERO multiplicity: %.0f - TPC multiplicity: %.0f",gRefMultiplicityVZERO,gRefMultiplicityTPC); |
| 1197 | |
| 1198 | fHistTPCvsVZEROMultiplicity->Fill(gRefMultiplicityVZERO,gRefMultiplicityTPC); |
| 1199 | |
| 1200 | if(fMultiplicityEstimator == "TPC") |
| 1201 | gRefMultiplicity = gRefMultiplicityTPC; |
| 1202 | else if(fMultiplicityEstimator == "V0M") |
| 1203 | gRefMultiplicity = gRefMultiplicityVZERO; |
| 1204 | else if(fMultiplicityEstimator == "V0A") |
| 1205 | gRefMultiplicity = gRefMultiplicityVZEROA; |
| 1206 | else if(fMultiplicityEstimator == "V0C") |
| 1207 | gRefMultiplicity = gRefMultiplicityVZEROC; |
| 1208 | |
| 1209 | return gRefMultiplicity; |
| 1210 | } |
| 1211 | |
| 1212 | //________________________________________________________________________ |
| 1213 | Double_t AliAnalysisTaskBFPsi::GetEventPlane(AliVEvent *event){ |
| 1214 | // Get the event plane |
| 1215 | |
| 1216 | TString gAnalysisLevel = fBalance->GetAnalysisLevel(); |
| 1217 | |
| 1218 | Float_t gVZEROEventPlane = -10.; |
| 1219 | Float_t gReactionPlane = -10.; |
| 1220 | Double_t qxTot = 0.0, qyTot = 0.0; |
| 1221 | |
| 1222 | //MC: from reaction plane |
| 1223 | if(gAnalysisLevel == "MC"){ |
| 1224 | if(!event) { |
| 1225 | AliError("mcEvent not available"); |
| 1226 | return 0x0; |
| 1227 | } |
| 1228 | |
| 1229 | AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event); |
| 1230 | if(gMCEvent){ |
| 1231 | AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader()); |
| 1232 | if (headerH) { |
| 1233 | gReactionPlane = headerH->ReactionPlaneAngle(); |
| 1234 | //gReactionPlane *= TMath::RadToDeg(); |
| 1235 | }//MC header |
| 1236 | }//MC event cast |
| 1237 | }//MC |
| 1238 | |
| 1239 | // AOD,ESD,ESDMC: from VZERO Event Plane |
| 1240 | else{ |
| 1241 | |
| 1242 | AliEventplane *ep = event->GetEventplane(); |
| 1243 | if(ep){ |
| 1244 | gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot); |
| 1245 | if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi(); |
| 1246 | //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg(); |
| 1247 | gReactionPlane = gVZEROEventPlane; |
| 1248 | } |
| 1249 | }//AOD,ESD,ESDMC |
| 1250 | |
| 1251 | return gReactionPlane; |
| 1252 | } |
| 1253 | |
| 1254 | //________________________________________________________________________ |
| 1255 | Double_t AliAnalysisTaskBFPsi::GetTrackbyTrackCorrectionMatrix( Double_t vEta, |
| 1256 | Double_t vPhi, |
| 1257 | Double_t vPt, |
| 1258 | Short_t vCharge, |
| 1259 | Double_t gCentrality) { |
| 1260 | // -- Get efficiency correction of particle dependent on (eta, phi, pt, charge, centrality) |
| 1261 | |
| 1262 | Double_t correction = 1.; |
| c683985a |
1263 | Int_t gCentralityInt = -1; |
| 967e4db5 |
1264 | |
| c683985a |
1265 | for (Int_t i=0; i<fCentralityArrayBinsForCorrections-1; i++){ |
| 1266 | if((fCentralityArrayForCorrections[i] <= gCentrality)&&(gCentrality <= fCentralityArrayForCorrections[i+1])){ |
| 1267 | gCentralityInt = i; |
| 1268 | break; |
| 1269 | } |
| 1270 | } |
| 1271 | |
| 1272 | // centrality not in array --> no correction |
| 1273 | if(gCentralityInt < 0){ |
| 1274 | correction = 1.; |
| 1275 | } |
| 1276 | else{ |
| 1277 | |
| 1278 | //Printf("//=============CENTRALITY=============// %d:",gCentralityInt); |
| 1279 | |
| 1280 | if(fHistCorrectionPlus[gCentralityInt]){ |
| 1281 | if (vCharge > 0) { |
| 967e4db5 |
1282 | correction = fHistCorrectionPlus[gCentralityInt]->GetBinContent(fHistCorrectionPlus[gCentralityInt]->FindBin(vEta,vPt,vPhi)); |
| c683985a |
1283 | //Printf("CORRECTIONplus: %.2f | Centrality %d",correction,gCentralityInt); |
| 1284 | } |
| 1285 | if (vCharge < 0) { |
| 967e4db5 |
1286 | correction = fHistCorrectionPlus[gCentralityInt]->GetBinContent(fHistCorrectionMinus[gCentralityInt]->FindBin(vEta,vPt,vPhi)); |
| 1287 | //Printf("CORRECTIONminus: %.2f | Centrality %d",correction,gCentralityInt); |
| c683985a |
1288 | } |
| 1289 | } |
| 1290 | else { |
| 1291 | correction = 1.; |
| 1292 | } |
| 1293 | }//centrality in array |
| 1294 | |
| 1295 | if (correction == 0.) { |
| 1296 | AliError(Form("Should not happen : bin content = 0. >> eta: %.2f | phi : %.2f | pt : %.2f | cent %d",vEta, vPhi, vPt, gCentralityInt)); |
| 1297 | correction = 1.; |
| 1298 | } |
| 1299 | |
| 1300 | return correction; |
| 1301 | } |
| 1302 | |
| 1303 | //________________________________________________________________________ |
| 1304 | TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t gCentrality, Double_t gReactionPlane){ |
| 1305 | // Returns TObjArray with tracks after all track cuts (only for AOD!) |
| 1306 | // Fills QA histograms |
| 1307 | |
| 1308 | TString gAnalysisLevel = fBalance->GetAnalysisLevel(); |
| 1309 | |
| 1310 | //output TObjArray holding all good tracks |
| 1311 | TObjArray* tracksAccepted = new TObjArray; |
| 1312 | tracksAccepted->SetOwner(kTRUE); |
| 1313 | |
| 1314 | Short_t vCharge; |
| 1315 | Double_t vEta; |
| 1316 | Double_t vY; |
| 1317 | Double_t vPhi; |
| 1318 | Double_t vPt; |
| 1319 | |
| 1320 | |
| 1321 | if(gAnalysisLevel == "AOD") { // handling of TPC only tracks different in AOD and ESD |
| 1322 | // Loop over tracks in event |
| 1323 | |
| 1324 | for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) { |
| 1325 | AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks)); |
| 1326 | if (!aodTrack) { |
| 1327 | AliError(Form("Could not receive track %d", iTracks)); |
| 1328 | continue; |
| 1329 | } |
| 1330 | |
| 1331 | // AOD track cuts |
| 1332 | |
| 1333 | // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C |
| 1334 | // take only TPC only tracks |
| 1335 | //fHistTrackStats->Fill(aodTrack->GetFilterMap()); |
| 1336 | for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){ |
| 1337 | fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit)); |
| 1338 | } |
| 1339 | |
| 1340 | if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue; |
| 1341 | |
| 1342 | |
| 1343 | // additional check on kPrimary flag |
| 1344 | if(fCheckPrimaryFlagAOD){ |
| 1345 | if(aodTrack->GetType() != AliAODTrack::kPrimary) |
| 1346 | continue; |
| 1347 | } |
| 1348 | |
| 1349 | |
| 1350 | vCharge = aodTrack->Charge(); |
| 1351 | vEta = aodTrack->Eta(); |
| 1352 | vY = aodTrack->Y(); |
| 1353 | vPhi = aodTrack->Phi();// * TMath::RadToDeg(); |
| 1354 | vPt = aodTrack->Pt(); |
| 1355 | |
| 1356 | //===========================PID (so far only for electron rejection)===============================// |
| 1357 | if(fElectronRejection) { |
| 1358 | |
| 1359 | // get the electron nsigma |
| 1360 | Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron)); |
| 1361 | |
| 1362 | //Fill QA before the PID |
| 1363 | fHistdEdxVsPTPCbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal()); |
| 1364 | fHistNSigmaTPCvsPtbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); |
| 1365 | //end of QA-before pid |
| 1366 | |
| 1367 | // check only for given momentum range |
| 1368 | if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){ |
| 1369 | |
| 1370 | //look only at electron nsigma |
| 1371 | if(!fElectronOnlyRejection){ |
| 1372 | |
| 1373 | //Make the decision based on the n-sigma of electrons |
| 1374 | if(nSigma < fElectronRejectionNSigma) continue; |
| 1375 | } |
| 1376 | else{ |
| 1377 | |
| 1378 | Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion)); |
| 1379 | Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon)); |
| 1380 | Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton)); |
| 1381 | |
| 1382 | //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species) |
| 1383 | if(nSigma < fElectronRejectionNSigma |
| 1384 | && nSigmaPions > fElectronRejectionNSigma |
| 1385 | && nSigmaKaons > fElectronRejectionNSigma |
| 1386 | && nSigmaProtons > fElectronRejectionNSigma ) continue; |
| 1387 | } |
| 1388 | } |
| 1389 | |
| 1390 | //Fill QA after the PID |
| 1391 | fHistdEdxVsPTPCafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal()); |
| 1392 | fHistNSigmaTPCvsPtafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); |
| 1393 | |
| 1394 | } |
| 1395 | //===========================end of PID (so far only for electron rejection)===============================// |
| 1396 | |
| 1397 | //+++++++++++++++++++++++++++++// |
| 1398 | //===========================PID===============================// |
| 1399 | if(fUsePID) { |
| 1400 | Double_t prob[AliPID::kSPECIES]={0.}; |
| 1401 | Double_t probTPC[AliPID::kSPECIES]={0.}; |
| 1402 | Double_t probTOF[AliPID::kSPECIES]={0.}; |
| 1403 | Double_t probTPCTOF[AliPID::kSPECIES]={0.}; |
| 1404 | |
| 1405 | Double_t nSigma = 0.; |
| 1406 | Double_t nSigmaTPC = 0.; //++++ |
| 1407 | Double_t nSigmaTOF = 0.; //++++ |
| 1408 | Double_t nSigmaTPCTOF = 0.; //++++ |
| 1409 | UInt_t detUsedTPC = 0; |
| 1410 | UInt_t detUsedTOF = 0; |
| 1411 | UInt_t detUsedTPCTOF = 0; |
| 1412 | |
| 1413 | nSigmaTPC = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)fParticleOfInterest)); |
| 1414 | nSigmaTOF = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(aodTrack,(AliPID::EParticleType)fParticleOfInterest)); |
| 1415 | nSigmaTPCTOF = TMath::Sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF); |
| 1416 | |
| 1417 | //Decide what detector configuration we want to use |
| 1418 | switch(fPidDetectorConfig) { |
| 1419 | case kTPCpid: |
| 1420 | fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC); |
| 1421 | nSigma = nSigmaTPC; |
| 1422 | detUsedTPC = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPC); |
| 1423 | for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) |
| 1424 | prob[iSpecies] = probTPC[iSpecies]; |
| 1425 | break; |
| 1426 | case kTOFpid: |
| 1427 | fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF); |
| 1428 | nSigma = nSigmaTOF; |
| 1429 | detUsedTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTOF); |
| 1430 | for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) |
| 1431 | prob[iSpecies] = probTOF[iSpecies]; |
| 1432 | break; |
| 1433 | case kTPCTOF: |
| 1434 | fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC); |
| 1435 | nSigma = nSigmaTPCTOF; |
| 1436 | detUsedTPCTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPCTOF); |
| 1437 | for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) |
| 1438 | prob[iSpecies] = probTPCTOF[iSpecies]; |
| 1439 | break; |
| 1440 | default: |
| 1441 | break; |
| 1442 | }//end switch: define detector mask |
| 1443 | |
| 1444 | //Filling the PID QA |
| 1445 | Double_t tofTime = -999., length = 999., tof = -999.; |
| 1446 | Double_t c = TMath::C()*1.E-9;// m/ns |
| 1447 | Double_t beta = -999.; |
| 1448 | if ( (aodTrack->IsOn(AliAODTrack::kTOFin)) && |
| 1449 | (aodTrack->IsOn(AliAODTrack::kTIME)) ) { |
| 1450 | tofTime = aodTrack->GetTOFsignal();//in ps |
| 1451 | length = aodTrack->GetIntegratedLength(); |
| 1452 | tof = tofTime*1E-3; // ns |
| 1453 | |
| 1454 | if (tof <= 0) { |
| 1455 | //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n"); |
| 1456 | continue; |
| 1457 | } |
| 1458 | if (length <= 0){ |
| 1459 | //printf("WARNING: track with negative length found!Skipping this track for PID checks\n"); |
| 1460 | continue; |
| 1461 | } |
| 1462 | |
| 1463 | length = length*0.01; // in meters |
| 1464 | tof = tof*c; |
| 1465 | beta = length/tof; |
| 1466 | |
| 1467 | fHistBetavsPTOFbeforePID ->Fill(aodTrack->P()*aodTrack->Charge(),beta); |
| 1468 | fHistProbTOFvsPtbeforePID ->Fill(aodTrack->Pt(),probTOF[fParticleOfInterest]); |
| 1469 | fHistNSigmaTOFvsPtbeforePID ->Fill(aodTrack->Pt(),nSigmaTOF); |
| 1470 | }//TOF signal |
| 1471 | |
| 1472 | fHistdEdxVsPTPCbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal()); |
| 1473 | fHistProbTPCvsPtbeforePID -> Fill(aodTrack->Pt(),probTPC[fParticleOfInterest]); |
| 1474 | fHistNSigmaTPCvsPtbeforePID -> Fill(aodTrack->Pt(),nSigmaTPC); |
| 1475 | |
| 1476 | fHistProbTPCTOFvsPtbeforePID -> Fill(aodTrack->Pt(),probTPCTOF[fParticleOfInterest]); |
| 1477 | |
| 1478 | //combined TPC&TOF |
| 1479 | fHistBetaVsdEdXbeforePID->Fill(aodTrack->GetTPCsignal(),beta); //+++++++++ |
| 1480 | fHistNSigmaTPCTOFvsPtbeforePID -> Fill(aodTrack->Pt(),nSigmaTPCTOF); |
| 1481 | Printf("NSIGMA: %lf - nSigmaTOF: %lf - nSigmaTPC: %lf - nSigmaTPCTOF: %lf",nSigma,nSigmaTOF,nSigmaTPC,nSigmaTPCTOF); |
| 1482 | |
| 1483 | //end of QA-before pid |
| 1484 | |
| 1485 | if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) { |
| 1486 | //Make the decision based on the n-sigma |
| 1487 | if(fUsePIDnSigma) { |
| 1488 | if(nSigma > fPIDNSigma) continue; |
| 1489 | |
| 1490 | fHistNSigmaTOFvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTOF); |
| 1491 | fHistNSigmaTPCvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTPC); |
| 1492 | fHistNSigmaTPCTOFvsPtafterPID ->Fill(aodTrack->Pt(),nSigmaTPCTOF); |
| 1493 | } |
| 1494 | //Make the decision based on the bayesian |
| 1495 | else if(fUsePIDPropabilities) { |
| 1496 | if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue; |
| 1497 | if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue; |
| 1498 | |
| 1499 | fHistProbTOFvsPtafterPID ->Fill(aodTrack->Pt(),probTOF[fParticleOfInterest]); |
| 1500 | fHistProbTPCvsPtafterPID ->Fill(aodTrack->Pt(),probTPC[fParticleOfInterest]); |
| 1501 | fHistProbTPCTOFvsPtafterPID ->Fill(aodTrack->Pt(),probTPCTOF[fParticleOfInterest]); |
| 1502 | |
| 1503 | } |
| 1504 | |
| 1505 | //Fill QA after the PID |
| 1506 | fHistBetavsPTOFafterPID ->Fill(aodTrack->P()*aodTrack->Charge(),beta); |
| 1507 | fHistdEdxVsPTPCafterPID ->Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal()); |
| 1508 | fHistBetaVsdEdXafterPID->Fill(aodTrack->GetTPCsignal(),beta); //+++++++++ |
| 1509 | } |
| 1510 | } |
| 1511 | //===========================PID===============================// |
| 1512 | //+++++++++++++++++++++++++++++// |
| 1513 | |
| 1514 | |
| 1515 | Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on) |
| 1516 | Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on) |
| 1517 | |
| 1518 | |
| 1519 | // Kinematics cuts from ESD track cuts |
| 1520 | if( vPt < fPtMin || vPt > fPtMax) continue; |
| 1521 | if( vEta < fEtaMin || vEta > fEtaMax) continue; |
| 1522 | |
| 1523 | // Extra DCA cuts (for systematic studies [!= -1]) |
| 1524 | if( fDCAxyCut != -1 && fDCAzCut != -1){ |
| 1525 | if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){ |
| 1526 | continue; // 2D cut |
| 1527 | } |
| 1528 | } |
| 1529 | |
| 1530 | // Extra TPC cuts (for systematic studies [!= -1]) |
| 1531 | if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){ |
| 1532 | continue; |
| 1533 | } |
| 1534 | if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){ |
| 1535 | continue; |
| 1536 | } |
| 1537 | |
| 1538 | // fill QA histograms |
| 1539 | fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls()); |
| 1540 | fHistDCA->Fill(dcaZ,dcaXY); |
| 1541 | fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality); |
| 1542 | fHistPt->Fill(vPt,gCentrality); |
| 1543 | fHistEta->Fill(vEta,gCentrality); |
| 1544 | fHistRapidity->Fill(vY,gCentrality); |
| 1545 | if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality); |
| 1546 | else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality); |
| 1547 | fHistPhi->Fill(vPhi,gCentrality); |
| 1548 | if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality); |
| 1549 | else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality); |
| 1550 | |
| 1551 | //=======================================correction |
| 1552 | Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality); |
| 1553 | //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality); |
| 1554 | |
| 1555 | // add the track to the TObjArray |
| 1556 | tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); |
| 1557 | }//track loop |
| 1558 | }// AOD analysis |
| 1559 | |
| 69d0337c |
1560 | |
| 6152cc3c |
1561 | // nano AODs |
| 1562 | else if(gAnalysisLevel == "AODnano") { // not fully supported yet (PID missing) |
| 1563 | // Loop over tracks in event |
| 1564 | |
| 1565 | for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) { |
| 69d0337c |
1566 | AliVTrack* aodTrack = dynamic_cast<AliVTrack *>(event->GetTrack(iTracks)); |
| 6152cc3c |
1567 | if (!aodTrack) { |
| 1568 | AliError(Form("Could not receive track %d", iTracks)); |
| 1569 | continue; |
| 1570 | } |
| 1571 | |
| 69d0337c |
1572 | // AOD track cuts (not needed) |
| 1573 | //if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue; |
| 6152cc3c |
1574 | |
| 1575 | vCharge = aodTrack->Charge(); |
| 1576 | vEta = aodTrack->Eta(); |
| 69d0337c |
1577 | vY = -999.; |
| 6152cc3c |
1578 | vPhi = aodTrack->Phi();// * TMath::RadToDeg(); |
| 1579 | vPt = aodTrack->Pt(); |
| 69d0337c |
1580 | |
| 6152cc3c |
1581 | |
| 1582 | // Kinematics cuts from ESD track cuts |
| 1583 | if( vPt < fPtMin || vPt > fPtMax) continue; |
| 1584 | if( vEta < fEtaMin || vEta > fEtaMax) continue; |
| 1585 | |
| 69d0337c |
1586 | |
| 6152cc3c |
1587 | // fill QA histograms |
| 6152cc3c |
1588 | fHistPt->Fill(vPt,gCentrality); |
| 1589 | fHistEta->Fill(vEta,gCentrality); |
| 1590 | fHistRapidity->Fill(vY,gCentrality); |
| 1591 | if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality); |
| 1592 | else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality); |
| 1593 | fHistPhi->Fill(vPhi,gCentrality); |
| 1594 | if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality); |
| 1595 | else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality); |
| 1596 | |
| 1597 | //=======================================correction |
| 1598 | Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality); |
| 1599 | //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality); |
| 1600 | |
| 1601 | // add the track to the TObjArray |
| 1602 | tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); |
| 1603 | }//track loop |
| 1604 | }// AOD nano analysis |
| 1605 | |
| 1606 | |
| c683985a |
1607 | //============================================================================================================== |
| 1608 | else if(gAnalysisLevel == "MCAOD") { |
| 1609 | |
| 1610 | AliMCEvent* mcEvent = MCEvent(); |
| 1611 | if (!mcEvent) { |
| 1612 | AliError("ERROR: Could not retrieve MC event"); |
| 1613 | } |
| 1614 | else{ |
| 1615 | |
| 1616 | for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) { |
| 1617 | AliAODMCParticle *aodTrack = (AliAODMCParticle*) mcEvent->GetTrack(iTracks); |
| 1618 | if (!aodTrack) { |
| 1619 | AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks)); |
| 1620 | continue; |
| 1621 | } |
| 1622 | |
| 1623 | if(!aodTrack->IsPhysicalPrimary()) continue; |
| 1624 | |
| 1625 | vCharge = aodTrack->Charge(); |
| 1626 | vEta = aodTrack->Eta(); |
| 1627 | vY = aodTrack->Y(); |
| 1628 | vPhi = aodTrack->Phi();// * TMath::RadToDeg(); |
| 1629 | vPt = aodTrack->Pt(); |
| 1630 | |
| 1631 | // Kinematics cuts from ESD track cuts |
| 1632 | if( vPt < fPtMin || vPt > fPtMax) continue; |
| 1633 | if( vEta < fEtaMin || vEta > fEtaMax) continue; |
| 1634 | |
| 1635 | // Remove neutral tracks |
| 1636 | if( vCharge == 0 ) continue; |
| 1637 | |
| 1638 | //Exclude resonances |
| 1639 | if(fExcludeResonancesInMC) { |
| 1640 | |
| 1641 | Bool_t kExcludeParticle = kFALSE; |
| 1642 | Int_t gMotherIndex = aodTrack->GetMother(); |
| 1643 | if(gMotherIndex != -1) { |
| 1644 | AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex)); |
| 1645 | if(motherTrack) { |
| 1646 | Int_t pdgCodeOfMother = motherTrack->GetPdgCode(); |
| 1647 | //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) { |
| 1648 | //if(pdgCodeOfMother == 113) { |
| 1649 | if(pdgCodeOfMother == 113 // rho0 |
| 1650 | || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+ |
| 1651 | // || pdgCodeOfMother == 221 // eta |
| 1652 | // || pdgCodeOfMother == 331 // eta' |
| 1653 | // || pdgCodeOfMother == 223 // omega |
| 1654 | // || pdgCodeOfMother == 333 // phi |
| 1655 | || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0 |
| 1656 | // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0* |
| 1657 | // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+* |
| 1658 | || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda |
| 1659 | || pdgCodeOfMother == 111 // pi0 Dalitz |
| 1660 | || pdgCodeOfMother == 22 // photon |
| 1661 | ) { |
| 1662 | kExcludeParticle = kTRUE; |
| 1663 | } |
| 1664 | } |
| 1665 | } |
| 1666 | |
| 1667 | //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi |
| 1668 | if(kExcludeParticle) continue; |
| 1669 | } |
| 1670 | |
| 1671 | //Exclude electrons with PDG |
| 1672 | if(fExcludeElectronsInMC) { |
| 1673 | |
| 1674 | if(TMath::Abs(aodTrack->GetPdgCode()) == 11) continue; |
| 1675 | |
| 1676 | } |
| 1677 | |
| 1678 | // fill QA histograms |
| 1679 | fHistPt->Fill(vPt,gCentrality); |
| 1680 | fHistEta->Fill(vEta,gCentrality); |
| 1681 | fHistRapidity->Fill(vY,gCentrality); |
| 1682 | if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality); |
| 1683 | else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality); |
| 1684 | fHistPhi->Fill(vPhi,gCentrality); |
| 1685 | if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality); |
| 1686 | else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality); |
| 1687 | |
| 1688 | //=======================================correction |
| 1689 | Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality); |
| 1690 | //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality); |
| 1691 | |
| 1692 | // add the track to the TObjArray |
| 1693 | tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); |
| 1694 | }//aodTracks |
| 1695 | }//MC event |
| 1696 | }//MCAOD |
| 1697 | //============================================================================================================== |
| 1698 | |
| 1699 | //============================================================================================================== |
| 1700 | else if(gAnalysisLevel == "MCAODrec") { |
| 1701 | |
| 1702 | /* fAOD = dynamic_cast<AliAODEvent*>(InputEvent()); |
| 1703 | if (!fAOD) { |
| 1704 | printf("ERROR: fAOD not available\n"); |
| 1705 | return; |
| 1706 | }*/ |
| 1707 | |
| 1708 | fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName())); |
| 1709 | if (!fArrayMC) { |
| 1710 | AliError("No array of MC particles found !!!"); |
| 1711 | } |
| 1712 | |
| 1713 | AliMCEvent* mcEvent = MCEvent(); |
| 1714 | if (!mcEvent) { |
| 1715 | AliError("ERROR: Could not retrieve MC event"); |
| 1716 | return tracksAccepted; |
| 1717 | } |
| 1718 | |
| 1719 | for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) { |
| 1720 | AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks)); |
| 1721 | if (!aodTrack) { |
| 1722 | AliError(Form("Could not receive track %d", iTracks)); |
| 1723 | continue; |
| 1724 | } |
| 1725 | |
| 1726 | for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){ |
| 1727 | fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit)); |
| 1728 | } |
| 1729 | if(!aodTrack->TestFilterBit(fnAODtrackCutBit)) continue; |
| 1730 | |
| 1731 | vCharge = aodTrack->Charge(); |
| 1732 | vEta = aodTrack->Eta(); |
| 1733 | vY = aodTrack->Y(); |
| 1734 | vPhi = aodTrack->Phi();// * TMath::RadToDeg(); |
| 1735 | vPt = aodTrack->Pt(); |
| 1736 | |
| 1737 | //===========================use MC information for Kinematics===============================// |
| 1738 | if(fUseMCforKinematics){ |
| 1739 | |
| 1740 | Int_t label = TMath::Abs(aodTrack->GetLabel()); |
| 1741 | AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label); |
| 1742 | |
| 1743 | if(AODmcTrack){ |
| 1744 | vCharge = AODmcTrack->Charge(); |
| 1745 | vEta = AODmcTrack->Eta(); |
| 1746 | vY = AODmcTrack->Y(); |
| 1747 | vPhi = AODmcTrack->Phi();// * TMath::RadToDeg(); |
| 1748 | vPt = AODmcTrack->Pt(); |
| 1749 | } |
| 1750 | else{ |
| 1751 | AliDebug(1, "no MC particle for this track"); |
| 1752 | continue; |
| 1753 | } |
| 1754 | } |
| 1755 | //===========================end of use MC information for Kinematics========================// |
| 1756 | |
| 1757 | |
| 1758 | //===========================PID (so far only for electron rejection)===============================// |
| 1759 | if(fElectronRejection) { |
| 1760 | |
| 1761 | // get the electron nsigma |
| 1762 | Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron)); |
| 1763 | |
| 1764 | //Fill QA before the PID |
| 1765 | fHistdEdxVsPTPCbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal()); |
| 1766 | fHistNSigmaTPCvsPtbeforePIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); |
| 1767 | //end of QA-before pid |
| 1768 | |
| 1769 | // check only for given momentum range |
| 1770 | if( vPt > fElectronRejectionMinPt && vPt < fElectronRejectionMaxPt ){ |
| 1771 | |
| 1772 | //look only at electron nsigma |
| 1773 | if(!fElectronOnlyRejection){ |
| 1774 | |
| 1775 | //Make the decision based on the n-sigma of electrons |
| 1776 | if(nSigma < fElectronRejectionNSigma) continue; |
| 1777 | } |
| 1778 | else{ |
| 1779 | |
| 1780 | Double_t nSigmaPions = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kPion)); |
| 1781 | Double_t nSigmaKaons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kKaon)); |
| 1782 | Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kProton)); |
| 1783 | |
| 1784 | //Make the decision based on the n-sigma of electrons exclusively ( = track not in nsigma region for other species) |
| 1785 | if(nSigma < fElectronRejectionNSigma |
| 1786 | && nSigmaPions > fElectronRejectionNSigma |
| 1787 | && nSigmaKaons > fElectronRejectionNSigma |
| 1788 | && nSigmaProtons > fElectronRejectionNSigma ) continue; |
| 1789 | } |
| 1790 | } |
| 1791 | |
| 1792 | //Fill QA after the PID |
| 1793 | fHistdEdxVsPTPCafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal()); |
| 1794 | fHistNSigmaTPCvsPtafterPIDelectron -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); |
| 1795 | |
| 1796 | } |
| 1797 | //===========================end of PID (so far only for electron rejection)===============================// |
| 1798 | |
| 1799 | Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on) |
| 1800 | Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on) |
| 1801 | |
| 1802 | // Kinematics cuts from ESD track cuts |
| 1803 | if( vPt < fPtMin || vPt > fPtMax) continue; |
| 1804 | if( vEta < fEtaMin || vEta > fEtaMax) continue; |
| 1805 | |
| 1806 | // Extra DCA cuts (for systematic studies [!= -1]) |
| 1807 | if( fDCAxyCut != -1 && fDCAzCut != -1){ |
| 1808 | if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){ |
| 1809 | continue; // 2D cut |
| 1810 | } |
| 1811 | } |
| 1812 | |
| 1813 | // Extra TPC cuts (for systematic studies [!= -1]) |
| 1814 | if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){ |
| 1815 | continue; |
| 1816 | } |
| 1817 | if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){ |
| 1818 | continue; |
| 1819 | } |
| 1820 | |
| 1821 | //Exclude resonances |
| 1822 | if(fExcludeResonancesInMC) { |
| 1823 | |
| 1824 | Bool_t kExcludeParticle = kFALSE; |
| 1825 | |
| 1826 | Int_t label = TMath::Abs(aodTrack->GetLabel()); |
| 1827 | AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label); |
| 1828 | |
| 1829 | if (AODmcTrack){ |
| 1830 | //if (AODmcTrack->IsPhysicalPrimary()){ |
| 1831 | |
| 1832 | Int_t gMotherIndex = AODmcTrack->GetMother(); |
| 1833 | if(gMotherIndex != -1) { |
| 1834 | AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex)); |
| 1835 | if(motherTrack) { |
| 1836 | Int_t pdgCodeOfMother = motherTrack->GetPdgCode(); |
| 1837 | if(pdgCodeOfMother == 113 // rho0 |
| 1838 | || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+ |
| 1839 | // || pdgCodeOfMother == 221 // eta |
| 1840 | // || pdgCodeOfMother == 331 // eta' |
| 1841 | // || pdgCodeOfMother == 223 // omega |
| 1842 | // || pdgCodeOfMother == 333 // phi |
| 1843 | || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0 |
| 1844 | // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0* |
| 1845 | // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+* |
| 1846 | || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda |
| 1847 | || pdgCodeOfMother == 111 // pi0 Dalitz |
| 1848 | || pdgCodeOfMother == 22 // photon |
| 1849 | ) { |
| 1850 | kExcludeParticle = kTRUE; |
| 1851 | } |
| 1852 | } |
| 1853 | } |
| 1854 | } |
| 1855 | //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi |
| 1856 | if(kExcludeParticle) continue; |
| 1857 | } |
| 1858 | |
| 1859 | //Exclude electrons with PDG |
| 1860 | if(fExcludeElectronsInMC) { |
| 1861 | |
| 1862 | Int_t label = TMath::Abs(aodTrack->GetLabel()); |
| 1863 | AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label); |
| 1864 | |
| 1865 | if (AODmcTrack){ |
| 1866 | if(TMath::Abs(AODmcTrack->GetPdgCode()) == 11) continue; |
| 1867 | } |
| 1868 | } |
| 1869 | |
| 1870 | // fill QA histograms |
| 1871 | fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls()); |
| 1872 | fHistDCA->Fill(dcaZ,dcaXY); |
| 1873 | fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality); |
| 1874 | fHistPt->Fill(vPt,gCentrality); |
| 1875 | fHistEta->Fill(vEta,gCentrality); |
| 1876 | fHistRapidity->Fill(vY,gCentrality); |
| 1877 | if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality); |
| 1878 | else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality); |
| 1879 | fHistPhi->Fill(vPhi,gCentrality); |
| 1880 | if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality); |
| 1881 | else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality); |
| 1882 | |
| 1883 | //=======================================correction |
| 1884 | Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality); |
| 1885 | //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality); |
| 1886 | |
| 1887 | // add the track to the TObjArray |
| 1888 | tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); |
| 1889 | }//track loop |
| 1890 | }//MCAODrec |
| 1891 | //============================================================================================================== |
| 1892 | |
| 1893 | else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD") { // handling of TPC only tracks different in AOD and ESD |
| 1894 | |
| 1895 | AliESDtrack *trackTPC = NULL; |
| 1896 | AliMCParticle *mcTrack = NULL; |
| 1897 | |
| 1898 | AliMCEvent* mcEvent = NULL; |
| 1899 | |
| 1900 | // for MC ESDs use also MC information |
| 1901 | if(gAnalysisLevel == "MCESD"){ |
| 1902 | mcEvent = MCEvent(); |
| 1903 | if (!mcEvent) { |
| 1904 | AliError("mcEvent not available"); |
| 1905 | return tracksAccepted; |
| 1906 | } |
| 1907 | } |
| 1908 | |
| 1909 | // Loop over tracks in event |
| 1910 | for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) { |
| 1911 | AliESDtrack* track = dynamic_cast<AliESDtrack *>(event->GetTrack(iTracks)); |
| 1912 | if (!track) { |
| 1913 | AliError(Form("Could not receive track %d", iTracks)); |
| 1914 | continue; |
| 1915 | } |
| 1916 | |
| 1917 | // for MC ESDs use also MC information --> MC track not used further??? |
| 1918 | if(gAnalysisLevel == "MCESD"){ |
| 1919 | Int_t label = TMath::Abs(track->GetLabel()); |
| 1920 | if(label > mcEvent->GetNumberOfTracks()) continue; |
| 1921 | if(label > mcEvent->GetNumberOfPrimaries()) continue; |
| 1922 | |
| 1923 | mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label)); |
| 1924 | if(!mcTrack) continue; |
| 1925 | } |
| 1926 | |
| 1927 | // take only TPC only tracks |
| 1928 | trackTPC = new AliESDtrack(); |
| 1929 | if(!track->FillTPCOnlyTrack(*trackTPC)) continue; |
| 1930 | |
| 1931 | //ESD track cuts |
| 1932 | if(fESDtrackCuts) |
| 1933 | if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue; |
| 1934 | |
| 1935 | // fill QA histograms |
| 1936 | Float_t b[2]; |
| 1937 | Float_t bCov[3]; |
| 1938 | trackTPC->GetImpactParameters(b,bCov); |
| 1939 | if (bCov[0]<=0 || bCov[2]<=0) { |
| 1940 | AliDebug(1, "Estimated b resolution lower or equal zero!"); |
| 1941 | bCov[0]=0; bCov[2]=0; |
| 1942 | } |
| 1943 | |
| 1944 | Int_t nClustersTPC = -1; |
| 1945 | nClustersTPC = trackTPC->GetTPCNclsIter1(); // TPC standalone |
| 1946 | //nClustersTPC = track->GetTPCclusters(0); // global track |
| 1947 | Float_t chi2PerClusterTPC = -1; |
| 1948 | if (nClustersTPC!=0) { |
| 1949 | chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone |
| 1950 | //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track |
| 1951 | } |
| 1952 | |
| 1953 | //===========================PID===============================// |
| 1954 | if(fUsePID) { |
| 1955 | Double_t prob[AliPID::kSPECIES]={0.}; |
| 1956 | Double_t probTPC[AliPID::kSPECIES]={0.}; |
| 1957 | Double_t probTOF[AliPID::kSPECIES]={0.}; |
| 1958 | Double_t probTPCTOF[AliPID::kSPECIES]={0.}; |
| 1959 | |
| 1960 | Double_t nSigma = 0.; |
| 1961 | UInt_t detUsedTPC = 0; |
| 1962 | UInt_t detUsedTOF = 0; |
| 1963 | UInt_t detUsedTPCTOF = 0; |
| 1964 | |
| 1965 | //Decide what detector configuration we want to use |
| 1966 | switch(fPidDetectorConfig) { |
| 1967 | case kTPCpid: |
| 1968 | fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC); |
| 1969 | nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest)); |
| 1970 | detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC); |
| 1971 | for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) |
| 1972 | prob[iSpecies] = probTPC[iSpecies]; |
| 1973 | break; |
| 1974 | case kTOFpid: |
| 1975 | fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF); |
| 1976 | nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest)); |
| 1977 | detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF); |
| 1978 | for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) |
| 1979 | prob[iSpecies] = probTOF[iSpecies]; |
| 1980 | break; |
| 1981 | case kTPCTOF: |
| 1982 | fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC); |
| 1983 | detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF); |
| 1984 | for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) |
| 1985 | prob[iSpecies] = probTPCTOF[iSpecies]; |
| 1986 | break; |
| 1987 | default: |
| 1988 | break; |
| 1989 | }//end switch: define detector mask |
| 1990 | |
| 1991 | //Filling the PID QA |
| 1992 | Double_t tofTime = -999., length = 999., tof = -999.; |
| 1993 | Double_t c = TMath::C()*1.E-9;// m/ns |
| 1994 | Double_t beta = -999.; |
| 1995 | Double_t nSigmaTOFForParticleOfInterest = -999.; |
| 1996 | if ( (track->IsOn(AliESDtrack::kTOFin)) && |
| 1997 | (track->IsOn(AliESDtrack::kTIME)) ) { |
| 1998 | tofTime = track->GetTOFsignal();//in ps |
| 1999 | length = track->GetIntegratedLength(); |
| 2000 | tof = tofTime*1E-3; // ns |
| 2001 | |
| 2002 | if (tof <= 0) { |
| 2003 | //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n"); |
| 2004 | continue; |
| 2005 | } |
| 2006 | if (length <= 0){ |
| 2007 | //printf("WARNING: track with negative length found!Skipping this track for PID checks\n"); |
| 2008 | continue; |
| 2009 | } |
| 2010 | |
| 2011 | length = length*0.01; // in meters |
| 2012 | tof = tof*c; |
| 2013 | beta = length/tof; |
| 2014 | |
| 2015 | nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest); |
| 2016 | fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta); |
| 2017 | fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]); |
| 2018 | fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest); |
| 2019 | }//TOF signal |
| 2020 | |
| 2021 | |
| 2022 | Double_t nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest); |
| 2023 | fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal()); |
| 2024 | fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); |
| 2025 | fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); |
| 2026 | fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]); |
| 2027 | //end of QA-before pid |
| 2028 | |
| 2029 | if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) { |
| 2030 | //Make the decision based on the n-sigma |
| 2031 | if(fUsePIDnSigma) { |
| 2032 | if(nSigma > fPIDNSigma) continue;} |
| 2033 | |
| 2034 | //Make the decision based on the bayesian |
| 2035 | else if(fUsePIDPropabilities) { |
| 2036 | if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue; |
| 2037 | if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue; |
| 2038 | } |
| 2039 | |
| 2040 | //Fill QA after the PID |
| 2041 | fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta); |
| 2042 | fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]); |
| 2043 | fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest); |
| 2044 | |
| 2045 | fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal()); |
| 2046 | fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); |
| 2047 | fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]); |
| 2048 | fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); |
| 2049 | } |
| 2050 | } |
| 2051 | //===========================PID===============================// |
| 2052 | vCharge = trackTPC->Charge(); |
| 2053 | vY = trackTPC->Y(); |
| 2054 | vEta = trackTPC->Eta(); |
| 2055 | vPhi = trackTPC->Phi();// * TMath::RadToDeg(); |
| 2056 | vPt = trackTPC->Pt(); |
| 2057 | fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC); |
| 2058 | fHistDCA->Fill(b[1],b[0]); |
| 2059 | fHistChi2->Fill(chi2PerClusterTPC,gCentrality); |
| 2060 | fHistPt->Fill(vPt,gCentrality); |
| 2061 | fHistEta->Fill(vEta,gCentrality); |
| 2062 | fHistPhi->Fill(vPhi,gCentrality); |
| 2063 | if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality); |
| 2064 | else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality); |
| 2065 | fHistRapidity->Fill(vY,gCentrality); |
| 2066 | if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality); |
| 2067 | else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality); |
| 2068 | |
| 2069 | //=======================================correction |
| 2070 | Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality); |
| 2071 | //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality); |
| 2072 | |
| 2073 | // add the track to the TObjArray |
| 2074 | tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); |
| 2075 | |
| 2076 | delete trackTPC; |
| 2077 | }//track loop |
| 2078 | }// ESD analysis |
| 2079 | |
| 2080 | else if(gAnalysisLevel == "MC"){ |
| 2081 | if(!event) { |
| 2082 | AliError("mcEvent not available"); |
| 2083 | return 0x0; |
| 2084 | } |
| 2085 | |
| 2086 | AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event); |
| 2087 | if(gMCEvent) { |
| 2088 | // Loop over tracks in event |
| 2089 | for (Int_t iTracks = 0; iTracks < gMCEvent->GetNumberOfPrimaries(); iTracks++) { |
| 2090 | AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iTracks)); |
| 2091 | if (!track) { |
| 2092 | AliError(Form("Could not receive particle %d", iTracks)); |
| 2093 | continue; |
| 2094 | } |
| 2095 | |
| 2096 | //exclude non stable particles |
| 2097 | if(!(gMCEvent->IsPhysicalPrimary(iTracks))) continue; |
| 2098 | |
| 2099 | vCharge = track->Charge(); |
| 2100 | vEta = track->Eta(); |
| 2101 | vPt = track->Pt(); |
| 2102 | vY = track->Y(); |
| 2103 | |
| 2104 | if( vPt < fPtMin || vPt > fPtMax) |
| 2105 | continue; |
| 2106 | if (!fUsePID) { |
| 2107 | if( vEta < fEtaMin || vEta > fEtaMax) continue; |
| 2108 | } |
| 2109 | else if (fUsePID){ |
| 2110 | if( vY < fEtaMin || vY > fEtaMax) continue; |
| 2111 | } |
| 2112 | |
| 2113 | // Remove neutral tracks |
| 2114 | if( vCharge == 0 ) continue; |
| 2115 | |
| 2116 | //analyze one set of particles |
| 2117 | if(fUseMCPdgCode) { |
| 2118 | TParticle *particle = track->Particle(); |
| 2119 | if(!particle) continue; |
| 2120 | |
| 2121 | Int_t gPdgCode = particle->GetPdgCode(); |
| 2122 | if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode)) |
| 2123 | continue; |
| 2124 | } |
| 2125 | |
| 2126 | //Use the acceptance parameterization |
| 2127 | if(fAcceptanceParameterization) { |
| 2128 | Double_t gRandomNumber = gRandom->Rndm(); |
| 2129 | if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt())) |
| 2130 | continue; |
| 2131 | } |
| 2132 | |
| 2133 | //Exclude resonances |
| 2134 | if(fExcludeResonancesInMC) { |
| 2135 | TParticle *particle = track->Particle(); |
| 2136 | if(!particle) continue; |
| 2137 | |
| 2138 | Bool_t kExcludeParticle = kFALSE; |
| 2139 | Int_t gMotherIndex = particle->GetFirstMother(); |
| 2140 | if(gMotherIndex != -1) { |
| 2141 | AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(event->GetTrack(gMotherIndex)); |
| 2142 | if(motherTrack) { |
| 2143 | TParticle *motherParticle = motherTrack->Particle(); |
| 2144 | if(motherParticle) { |
| 2145 | Int_t pdgCodeOfMother = motherParticle->GetPdgCode(); |
| 2146 | //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) { |
| 2147 | if(pdgCodeOfMother == 113 // rho0 |
| 2148 | || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+ |
| 2149 | // || pdgCodeOfMother == 221 // eta |
| 2150 | // || pdgCodeOfMother == 331 // eta' |
| 2151 | // || pdgCodeOfMother == 223 // omega |
| 2152 | // || pdgCodeOfMother == 333 // phi |
| 2153 | || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0 |
| 2154 | // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0* |
| 2155 | // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+* |
| 2156 | || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda |
| 2157 | || pdgCodeOfMother == 111 // pi0 Dalitz |
| 2158 | ) { |
| 2159 | kExcludeParticle = kTRUE; |
| 2160 | } |
| 2161 | } |
| 2162 | } |
| 2163 | } |
| 2164 | |
| 2165 | //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi |
| 2166 | if(kExcludeParticle) continue; |
| 2167 | } |
| 2168 | |
| 2169 | //Exclude electrons with PDG |
| 2170 | if(fExcludeElectronsInMC) { |
| 2171 | |
| 2172 | TParticle *particle = track->Particle(); |
| 2173 | |
| 2174 | if (particle){ |
| 2175 | if(TMath::Abs(particle->GetPdgCode()) == 11) continue; |
| 2176 | } |
| 2177 | } |
| 2178 | |
| 2179 | vPhi = track->Phi(); |
| 2180 | //Printf("phi (before): %lf",vPhi); |
| 2181 | |
| 2182 | fHistPt->Fill(vPt,gCentrality); |
| 2183 | fHistEta->Fill(vEta,gCentrality); |
| 2184 | fHistPhi->Fill(vPhi,gCentrality); |
| 2185 | if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality); |
| 2186 | else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality); |
| 2187 | //fHistPhi->Fill(vPhi*TMath::RadToDeg(),gCentrality); |
| 2188 | fHistRapidity->Fill(vY,gCentrality); |
| 2189 | //if(vCharge > 0) fHistPhiPos->Fill(vPhi*TMath::RadToDeg(),gCentrality); |
| 2190 | //else if(vCharge < 0) fHistPhiNeg->Fill(vPhi*TMath::RadToDeg(),gCentrality); |
| 2191 | if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality); |
| 2192 | else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality); |
| 2193 | |
| 2194 | //Flow after burner |
| 2195 | if(fUseFlowAfterBurner) { |
| 2196 | Double_t precisionPhi = 0.001; |
| 2197 | Int_t maxNumberOfIterations = 100; |
| 2198 | |
| 2199 | Double_t phi0 = vPhi; |
| 2200 | Double_t gV2 = fDifferentialV2->Eval(vPt); |
| 2201 | |
| 2202 | for (Int_t j = 0; j < maxNumberOfIterations; j++) { |
| 2203 | Double_t phiprev = vPhi; |
| 2204 | Double_t fl = vPhi - phi0 + gV2*TMath::Sin(2.*(vPhi - gReactionPlane*TMath::DegToRad())); |
| 2205 | Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(vPhi - gReactionPlane*TMath::DegToRad())); |
| 2206 | vPhi -= fl/fp; |
| 2207 | if (TMath::AreEqualAbs(phiprev,vPhi,precisionPhi)) break; |
| 2208 | } |
| 2209 | //Printf("phi (after): %lf\n",vPhi); |
| 2210 | Double_t vDeltaphiBefore = phi0 - gReactionPlane*TMath::DegToRad(); |
| 2211 | if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi(); |
| 2212 | fHistPhiBefore->Fill(vDeltaphiBefore,gCentrality); |
| 2213 | |
| 2214 | Double_t vDeltaphiAfter = vPhi - gReactionPlane*TMath::DegToRad(); |
| 2215 | if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi(); |
| 2216 | fHistPhiAfter->Fill(vDeltaphiAfter,gCentrality); |
| 2217 | |
| 2218 | } |
| 2219 | |
| 2220 | //vPhi *= TMath::RadToDeg(); |
| 2221 | |
| 2222 | //=======================================correction |
| 2223 | Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality); |
| 2224 | //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality); |
| 2225 | |
| 2226 | tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); |
| 2227 | } //track loop |
| 2228 | }//MC event object |
| 2229 | }//MC |
| 2230 | |
| 2231 | return tracksAccepted; |
| 2232 | } |
| 2233 | |
| 2234 | //________________________________________________________________________ |
| 2235 | TObjArray* AliAnalysisTaskBFPsi::GetShuffledTracks(TObjArray *tracks, Double_t gCentrality){ |
| 2236 | // Clones TObjArray and returns it with tracks after shuffling the charges |
| 2237 | |
| 2238 | TObjArray* tracksShuffled = new TObjArray; |
| 2239 | tracksShuffled->SetOwner(kTRUE); |
| 2240 | |
| 2241 | vector<Short_t> *chargeVector = new vector<Short_t>; //original charge of accepted tracks |
| 2242 | |
| 2243 | for (Int_t i=0; i<tracks->GetEntriesFast(); i++) |
| 2244 | { |
| 2245 | AliVParticle* track = (AliVParticle*) tracks->At(i); |
| 2246 | chargeVector->push_back(track->Charge()); |
| 2247 | } |
| 2248 | |
| 2249 | random_shuffle(chargeVector->begin(), chargeVector->end()); |
| 2250 | |
| 2251 | for(Int_t i = 0; i < tracks->GetEntriesFast(); i++){ |
| 2252 | AliVParticle* track = (AliVParticle*) tracks->At(i); |
| 2253 | //==============================correction |
| 2254 | Double_t correction = GetTrackbyTrackCorrectionMatrix(track->Eta(), track->Phi(),track->Pt(), chargeVector->at(i), gCentrality); |
| 2255 | //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality); |
| 2256 | tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i), correction)); |
| 2257 | } |
| 2258 | |
| 2259 | delete chargeVector; |
| 2260 | |
| 2261 | return tracksShuffled; |
| 2262 | } |
| 2263 | |
| 2264 | //________________________________________________________________________ |
| 2265 | void AliAnalysisTaskBFPsi::SetVZEROCalibrationFile(const char* filename, |
| 2266 | const char* lhcPeriod) { |
| 2267 | //Function to setup the VZERO gain equalization |
| 2268 | //============Get the equilization map============// |
| 2269 | TFile *calibrationFile = TFile::Open(filename); |
| 2270 | if((!calibrationFile)||(!calibrationFile->IsOpen())) { |
| 2271 | Printf("No calibration file found!!!"); |
| 2272 | return; |
| 2273 | } |
| 2274 | |
| 2275 | TList *list = dynamic_cast<TList *>(calibrationFile->Get(lhcPeriod)); |
| 2276 | if(!list) { |
| 2277 | Printf("Calibration TList not found!!!"); |
| 2278 | return; |
| 2279 | } |
| 2280 | |
| 2281 | fHistVZEROAGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROAGainEqualizationMap")); |
| 2282 | if(!fHistVZEROAGainEqualizationMap) { |
| 2283 | Printf("VZERO-A calibration object not found!!!"); |
| 2284 | return; |
| 2285 | } |
| 2286 | fHistVZEROCGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROCGainEqualizationMap")); |
| 2287 | if(!fHistVZEROCGainEqualizationMap) { |
| 2288 | Printf("VZERO-C calibration object not found!!!"); |
| 2289 | return; |
| 2290 | } |
| 2291 | |
| 2292 | fHistVZEROChannelGainEqualizationMap = dynamic_cast<TH2F *>(list->FindObject("gHistVZEROChannelGainEqualizationMap")); |
| 2293 | if(!fHistVZEROChannelGainEqualizationMap) { |
| 2294 | Printf("VZERO channel calibration object not found!!!"); |
| 2295 | return; |
| 2296 | } |
| 2297 | } |
| 2298 | |
| 2299 | //________________________________________________________________________ |
| 2300 | Double_t AliAnalysisTaskBFPsi::GetChannelEqualizationFactor(Int_t run, |
| 2301 | Int_t channel) { |
| 2302 | // |
| 2303 | if(!fHistVZEROAGainEqualizationMap) return 1.0; |
| 2304 | |
| 2305 | for(Int_t iBinX = 1; iBinX <= fHistVZEROChannelGainEqualizationMap->GetNbinsX(); iBinX++) { |
| 2306 | Int_t gRunNumber = atoi(fHistVZEROChannelGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX)); |
| 2307 | if(gRunNumber == run) |
| 2308 | return fHistVZEROChannelGainEqualizationMap->GetBinContent(iBinX,channel+1); |
| 2309 | } |
| 2310 | |
| 2311 | return 1.0; |
| 2312 | } |
| 2313 | |
| 2314 | //________________________________________________________________________ |
| 2315 | Double_t AliAnalysisTaskBFPsi::GetEqualizationFactor(Int_t run, |
| 2316 | const char* side) { |
| 2317 | // |
| 2318 | if(!fHistVZEROAGainEqualizationMap) return 1.0; |
| 2319 | |
| 2320 | TString gVZEROSide = side; |
| 2321 | for(Int_t iBinX = 1; iBinX < fHistVZEROAGainEqualizationMap->GetNbinsX(); iBinX++) { |
| 2322 | Int_t gRunNumber = atoi(fHistVZEROAGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX)); |
| 2323 | //cout<<"Looking for run "<<run<<" - current run: "<<gRunNumber<<endl; |
| 2324 | if(gRunNumber == run) { |
| 2325 | if(gVZEROSide == "A") |
| 2326 | return fHistVZEROAGainEqualizationMap->GetBinContent(iBinX); |
| 2327 | else if(gVZEROSide == "C") |
| 2328 | return fHistVZEROCGainEqualizationMap->GetBinContent(iBinX); |
| 2329 | } |
| 2330 | } |
| 2331 | |
| 2332 | return 1.0; |
| 2333 | } |
| 2334 | |
| 2335 | //____________________________________________________________________ |
| 2336 | Bool_t AliAnalysisTaskBFPsi::AcceptEventCentralityWeight(Double_t centrality) |
| 2337 | { |
| 2338 | // copied from AliAnalysisTaskPhiCorrelations |
| 2339 | // |
| 2340 | // rejects "randomly" events such that the centrality gets flat |
| 2341 | // uses fCentralityWeights histogram |
| 2342 | |
| 2343 | // TODO code taken and adapted from AliRDHFCuts; waiting for general class AliCentralityFlattening |
| 2344 | |
| 2345 | Double_t weight = fCentralityWeights->GetBinContent(fCentralityWeights->FindBin(centrality)); |
| 2346 | Double_t centralityDigits = centrality*100. - (Int_t)(centrality*100.); |
| 2347 | |
| 2348 | Bool_t result = kFALSE; |
| 2349 | if (centralityDigits < weight) |
| 2350 | result = kTRUE; |
| 2351 | |
| 2352 | AliInfo(Form("Centrality: %f; Digits: %f; Weight: %f; Result: %d", centrality, centralityDigits, weight, result)); |
| 2353 | |
| 2354 | return result; |
| 2355 | } |
| 2356 | |
| 2357 | //________________________________________________________________________ |
| 2358 | void AliAnalysisTaskBFPsi::FinishTaskOutput(){ |
| 2359 | //Printf("END BF"); |
| 2360 | |
| 2361 | if (!fBalance) { |
| 2362 | AliError("fBalance not available"); |
| 2363 | return; |
| 2364 | } |
| 2365 | if(fRunShuffling) { |
| 2366 | if (!fShuffledBalance) { |
| 2367 | AliError("fShuffledBalance not available"); |
| 2368 | return; |
| 2369 | } |
| 2370 | } |
| 2371 | |
| 2372 | } |
| 2373 | |
| 2374 | //________________________________________________________________________ |
| 2375 | void AliAnalysisTaskBFPsi::Terminate(Option_t *) { |
| 2376 | // Draw result to the screen |
| 2377 | // Called once at the end of the query |
| 2378 | |
| 2379 | // not implemented ... |
| 2380 | |
| 2381 | } |
| 2382 | |
git://git.uio.no / u / mrichter / AliRoot.git / blame