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