| b8057da6 |
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 "TArrayF.h" |
| 9 | #include "TF1.h" |
| 10 | #include "TRandom.h" |
| 11 | |
| 12 | #include "AliAnalysisTaskSE.h" |
| 13 | #include "AliAnalysisManager.h" |
| 14 | |
| 15 | #include "AliESDVertex.h" |
| 16 | #include "AliESDEvent.h" |
| 17 | #include "AliESDInputHandler.h" |
| 18 | #include "AliAODEvent.h" |
| 19 | #include "AliAODTrack.h" |
| 20 | #include "AliAODInputHandler.h" |
| 21 | #include "AliGenEventHeader.h" |
| 22 | #include "AliGenHijingEventHeader.h" |
| 23 | #include "AliMCEventHandler.h" |
| 24 | #include "AliMCEvent.h" |
| 25 | #include "AliStack.h" |
| 26 | #include "AliESDtrackCuts.h" |
| 27 | |
| 28 | #include "TH2D.h" |
| 29 | #include "AliPID.h" |
| 30 | #include "AliPIDResponse.h" |
| 31 | #include "AliPIDCombined.h" |
| 32 | |
| 33 | #include "AliAnalysisTaskBF.h" |
| 34 | #include "AliBalance.h" |
| 35 | |
| 36 | |
| 37 | // Analysis task for the BF code |
| 38 | // Authors: Panos.Christakoglou@nikhef.nl |
| 39 | |
| 40 | ClassImp(AliAnalysisTaskBF) |
| 41 | |
| 42 | //________________________________________________________________________ |
| 43 | AliAnalysisTaskBF::AliAnalysisTaskBF(const char *name) |
| 44 | : AliAnalysisTaskSE(name), |
| 45 | fBalance(0), |
| 46 | fRunShuffling(kFALSE), |
| 47 | fShuffledBalance(0), |
| 48 | fList(0), |
| 49 | fListBF(0), |
| 50 | fListBFS(0), |
| 51 | fHistListPIDQA(0), |
| 52 | fHistEventStats(0), |
| 53 | fHistCentStats(0), |
| 54 | fHistTriggerStats(0), |
| 55 | fHistTrackStats(0), |
| 56 | fHistVx(0), |
| 57 | fHistVy(0), |
| 58 | fHistVz(0), |
| 59 | fHistClus(0), |
| 60 | fHistDCA(0), |
| 61 | fHistChi2(0), |
| 62 | fHistPt(0), |
| 63 | fHistEta(0), |
| 64 | fHistRapidity(0), |
| 65 | fHistPhi(0), |
| 66 | fHistPhiBefore(0), |
| 67 | fHistPhiAfter(0), |
| 68 | fHistPhiPos(0), |
| 69 | fHistPhiNeg(0), |
| 70 | fHistV0M(0), |
| 71 | fHistRefTracks(0), |
| 72 | fHistdEdxVsPTPCbeforePID(NULL), |
| 73 | fHistBetavsPTOFbeforePID(NULL), |
| 74 | fHistProbTPCvsPtbeforePID(NULL), |
| 75 | fHistProbTOFvsPtbeforePID(NULL), |
| 76 | fHistProbTPCTOFvsPtbeforePID(NULL), |
| 77 | fHistNSigmaTPCvsPtbeforePID(NULL), |
| 78 | fHistNSigmaTOFvsPtbeforePID(NULL), |
| 79 | fHistdEdxVsPTPCafterPID(NULL), |
| 80 | fHistBetavsPTOFafterPID(NULL), |
| 81 | fHistProbTPCvsPtafterPID(NULL), |
| 82 | fHistProbTOFvsPtafterPID(NULL), |
| 83 | fHistProbTPCTOFvsPtafterPID(NULL), |
| 84 | fHistNSigmaTPCvsPtafterPID(NULL), |
| 85 | fHistNSigmaTOFvsPtafterPID(NULL), |
| 86 | fPIDResponse(0x0), |
| 87 | fPIDCombined(0x0), |
| 88 | fParticleOfInterest(kPion), |
| 89 | fPidDetectorConfig(kTPCTOF), |
| 90 | fUsePID(kFALSE), |
| 91 | fUsePIDnSigma(kTRUE), |
| 92 | fUsePIDPropabilities(kFALSE), |
| 93 | fPIDNSigma(3.), |
| 94 | fMinAcceptedPIDProbability(0.8), |
| 95 | fESDtrackCuts(0), |
| 96 | fCentralityEstimator("V0M"), |
| 97 | fUseCentrality(kFALSE), |
| 98 | fCentralityPercentileMin(0.), |
| 99 | fCentralityPercentileMax(5.), |
| 100 | fImpactParameterMin(0.), |
| 101 | fImpactParameterMax(20.), |
| 102 | fUseMultiplicity(kFALSE), |
| 103 | fNumberOfAcceptedTracksMin(0), |
| 104 | fNumberOfAcceptedTracksMax(10000), |
| 105 | fHistNumberOfAcceptedTracks(0), |
| 106 | fUseOfflineTrigger(kFALSE), |
| 107 | fVxMax(0.3), |
| 108 | fVyMax(0.3), |
| 109 | fVzMax(10.), |
| 9fd4b54e |
110 | fAODtrackCutBit(128), |
| b8057da6 |
111 | fPtMin(0.3), |
| 112 | fPtMax(1.5), |
| 113 | fEtaMin(-0.8), |
| 114 | fEtaMax(-0.8), |
| 115 | fDCAxyCut(-1), |
| 116 | fDCAzCut(-1), |
| 117 | fTPCchi2Cut(-1), |
| 118 | fNClustersTPCCut(-1), |
| 119 | fAcceptanceParameterization(0), |
| 120 | fDifferentialV2(0), |
| 121 | fUseFlowAfterBurner(kFALSE), |
| 122 | fExcludeResonancesInMC(kFALSE), |
| 123 | fUseMCPdgCode(kFALSE), |
| 124 | fPDGCodeToBeAnalyzed(-1) { |
| 125 | // Constructor |
| 126 | // Define input and output slots here |
| 127 | // Input slot #0 works with a TChain |
| 128 | DefineInput(0, TChain::Class()); |
| 129 | // Output slot #0 writes into a TH1 container |
| 130 | DefineOutput(1, TList::Class()); |
| 131 | DefineOutput(2, TList::Class()); |
| 132 | DefineOutput(3, TList::Class()); |
| 133 | DefineOutput(4, TList::Class()); |
| 134 | } |
| 135 | |
| 136 | //________________________________________________________________________ |
| 137 | AliAnalysisTaskBF::~AliAnalysisTaskBF() { |
| 138 | |
| 139 | // delete fBalance; |
| 140 | // delete fShuffledBalance; |
| 141 | // delete fList; |
| 142 | // delete fListBF; |
| 143 | // delete fListBFS; |
| 144 | |
| 145 | // delete fHistEventStats; |
| 146 | // delete fHistTrackStats; |
| 147 | // delete fHistVx; |
| 148 | // delete fHistVy; |
| 149 | // delete fHistVz; |
| 150 | |
| 151 | // delete fHistClus; |
| 152 | // delete fHistDCA; |
| 153 | // delete fHistChi2; |
| 154 | // delete fHistPt; |
| 155 | // delete fHistEta; |
| 156 | // delete fHistPhi; |
| 157 | // delete fHistV0M; |
| 158 | } |
| 159 | |
| 160 | //________________________________________________________________________ |
| 161 | void AliAnalysisTaskBF::UserCreateOutputObjects() { |
| 162 | // Create histograms |
| 163 | // Called once |
| 164 | if(!fBalance) { |
| 165 | fBalance = new AliBalance(); |
| 166 | fBalance->SetAnalysisLevel("ESD"); |
| 167 | //fBalance->SetNumberOfBins(-1,16); |
| 168 | fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6); |
| 169 | } |
| 170 | if(fRunShuffling) { |
| 171 | if(!fShuffledBalance) { |
| 172 | fShuffledBalance = new AliBalance(); |
| 173 | fShuffledBalance->SetAnalysisLevel("ESD"); |
| 174 | //fShuffledBalance->SetNumberOfBins(-1,16); |
| 175 | fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6); |
| 176 | } |
| 177 | } |
| 178 | |
| 179 | //QA list |
| 180 | fList = new TList(); |
| 181 | fList->SetName("listQA"); |
| 182 | fList->SetOwner(); |
| 183 | |
| 184 | //Balance Function list |
| 185 | fListBF = new TList(); |
| 186 | fListBF->SetName("listBF"); |
| 187 | fListBF->SetOwner(); |
| 188 | |
| 189 | if(fRunShuffling) { |
| 190 | fListBFS = new TList(); |
| 191 | fListBFS->SetName("listBFShuffled"); |
| 192 | fListBFS->SetOwner(); |
| 193 | } |
| 194 | |
| 195 | //PID QA list |
| 196 | if(fUsePID) { |
| 197 | fHistListPIDQA = new TList(); |
| 198 | fHistListPIDQA->SetName("listQAPID"); |
| 199 | fHistListPIDQA->SetOwner(); |
| 200 | } |
| 201 | |
| 202 | //Event stats. |
| 203 | TString gCutName[4] = {"Total","Offline trigger", |
| 204 | "Vertex","Analyzed"}; |
| 205 | fHistEventStats = new TH1F("fHistEventStats", |
| 206 | "Event statistics;;N_{events}", |
| 207 | 4,0.5,4.5); |
| 208 | for(Int_t i = 1; i <= 4; i++) |
| 209 | fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data()); |
| 210 | fList->Add(fHistEventStats); |
| 211 | |
| 212 | TString gCentName[9] = {"V0M","FMD","TRK","TKL","CL0","CL1","V0MvsFMD","TKLvsV0M","ZEMvsZDC"}; |
| 213 | fHistCentStats = new TH2F("fHistCentStats", |
| 214 | "Centrality statistics;;Cent percentile", |
| 215 | 9,-0.5,8.5,220,-5,105); |
| 216 | for(Int_t i = 1; i <= 9; i++) |
| 217 | fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data()); |
| 218 | fList->Add(fHistCentStats); |
| 219 | |
| 220 | fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",130,0,130); |
| 221 | fList->Add(fHistTriggerStats); |
| 222 | |
| 223 | fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",130,0,130); |
| 224 | fList->Add(fHistTrackStats); |
| 225 | |
| 226 | fHistNumberOfAcceptedTracks = new TH1F("fHistNumberOfAcceptedTracks",";N_{acc.};Entries",4001,-0.5,4000.5); |
| 227 | fList->Add(fHistNumberOfAcceptedTracks); |
| 228 | |
| 229 | // Vertex distributions |
| 230 | fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5); |
| 231 | fList->Add(fHistVx); |
| 232 | fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5); |
| 233 | fList->Add(fHistVy); |
| 234 | fHistVz = new TH1F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Entries",100,-20.,20.); |
| 235 | fList->Add(fHistVz); |
| 236 | |
| 237 | // QA histograms |
| 238 | fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200); |
| 239 | fList->Add(fHistClus); |
| 240 | fHistChi2 = new TH1F("fHistChi2","Chi2/NDF distribution",200,0,10); |
| 241 | fList->Add(fHistChi2); |
| 242 | fHistDCA = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5); |
| 243 | fList->Add(fHistDCA); |
| 244 | fHistPt = new TH1F("fHistPt","p_{T} distribution",200,0,10); |
| 245 | fList->Add(fHistPt); |
| 246 | fHistEta = new TH1F("fHistEta","#eta distribution",200,-2,2); |
| 247 | fList->Add(fHistEta); |
| 248 | fHistRapidity = new TH1F("fHistRapidity","y distribution",200,-2,2); |
| 249 | fList->Add(fHistRapidity); |
| 250 | fHistPhi = new TH1F("fHistPhi","#phi distribution",200,-20,380); |
| 251 | fList->Add(fHistPhi); |
| 252 | fHistPhiBefore = new TH1F("fHistPhiBefore","#phi distribution",200,0.,2*TMath::Pi()); |
| 253 | fList->Add(fHistPhiBefore); |
| 254 | fHistPhiAfter = new TH1F("fHistPhiAfter","#phi distribution",200,0.,2*TMath::Pi()); |
| 255 | fList->Add(fHistPhiAfter); |
| 256 | fHistPhiPos = new TH1F("fHistPhiPos","#phi distribution for positive particles",200,-20,380); |
| 257 | fList->Add(fHistPhiPos); |
| 258 | fHistPhiNeg = new TH1F("fHistPhiNeg","#phi distribution for negative particles",200,-20,380); |
| 259 | fList->Add(fHistPhiNeg); |
| 260 | fHistV0M = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000); |
| 261 | fList->Add(fHistV0M); |
| 262 | TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"}; |
| 263 | fHistRefTracks = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000); |
| 264 | for(Int_t i = 1; i <= 6; i++) |
| 265 | fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data()); |
| 266 | fList->Add(fHistRefTracks); |
| 267 | |
| 6de53600 |
268 | // QA histograms for HBTinspired and Conversion cuts |
| 269 | fList->Add(fBalance->GetQAHistHBTbefore()); |
| 270 | fList->Add(fBalance->GetQAHistHBTafter()); |
| 271 | fList->Add(fBalance->GetQAHistConversionbefore()); |
| 272 | fList->Add(fBalance->GetQAHistConversionafter()); |
| 273 | |
| b8057da6 |
274 | // Balance function histograms |
| 275 | // Initialize histograms if not done yet |
| 276 | if(!fBalance->GetHistNp(0)){ |
| 277 | AliWarning("Histograms not yet initialized! --> Will be done now"); |
| 278 | AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction"); |
| 279 | fBalance->InitHistograms(); |
| 280 | } |
| 281 | |
| 282 | if(fRunShuffling) { |
| 283 | if(!fShuffledBalance->GetHistNp(0)) { |
| 284 | AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now"); |
| 285 | AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction"); |
| 286 | fShuffledBalance->InitHistograms(); |
| 287 | } |
| 288 | } |
| 289 | |
| 290 | for(Int_t a = 0; a < ANALYSIS_TYPES; a++){ |
| 291 | fListBF->Add(fBalance->GetHistNp(a)); |
| 292 | fListBF->Add(fBalance->GetHistNn(a)); |
| 293 | fListBF->Add(fBalance->GetHistNpn(a)); |
| 294 | fListBF->Add(fBalance->GetHistNnn(a)); |
| 295 | fListBF->Add(fBalance->GetHistNpp(a)); |
| 296 | fListBF->Add(fBalance->GetHistNnp(a)); |
| 297 | |
| 298 | if(fRunShuffling) { |
| 299 | fListBFS->Add(fShuffledBalance->GetHistNp(a)); |
| 300 | fListBFS->Add(fShuffledBalance->GetHistNn(a)); |
| 301 | fListBFS->Add(fShuffledBalance->GetHistNpn(a)); |
| 302 | fListBFS->Add(fShuffledBalance->GetHistNnn(a)); |
| 303 | fListBFS->Add(fShuffledBalance->GetHistNpp(a)); |
| 304 | fListBFS->Add(fShuffledBalance->GetHistNnp(a)); |
| 305 | } |
| 6de53600 |
306 | } |
| b8057da6 |
307 | |
| 308 | if(fESDtrackCuts) fList->Add(fESDtrackCuts); |
| 309 | |
| 310 | //====================PID========================// |
| 311 | if(fUsePID) { |
| 312 | fPIDCombined = new AliPIDCombined(); |
| 313 | fPIDCombined->SetDefaultTPCPriors(); |
| 314 | |
| 315 | fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000); |
| 316 | fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID); //addition |
| 317 | |
| 318 | fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2); |
| 319 | fHistListPIDQA->Add(fHistBetavsPTOFbeforePID); //addition |
| 320 | |
| 321 | fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0); |
| 322 | fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID); //addition |
| 323 | |
| 324 | fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); |
| 325 | fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID); //addition |
| 326 | |
| 327 | fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); |
| 328 | fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID); //addition |
| 329 | |
| 330 | fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500); |
| 331 | fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID); //addition |
| 332 | |
| 333 | fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, 0, 500); |
| 334 | fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID); //addition |
| 335 | |
| 336 | fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000); |
| 337 | fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID); //addition |
| 338 | |
| 339 | fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2); |
| 340 | fHistListPIDQA->Add(fHistBetavsPTOFafterPID); //addition |
| 341 | |
| 342 | fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2); |
| 343 | fHistListPIDQA->Add(fHistProbTPCvsPtafterPID); //addition |
| 344 | |
| 345 | fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000, -10, 10, 1000, 0, 2); |
| 346 | fHistListPIDQA->Add(fHistProbTOFvsPtafterPID); //addition |
| 347 | |
| 348 | fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0); |
| 349 | fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID); //addition |
| 350 | |
| 351 | fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500); |
| 352 | fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID); //addition |
| 353 | |
| 354 | fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, 0, 500); |
| 355 | fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID); //addition |
| 356 | } |
| 357 | //====================PID========================// |
| 358 | |
| 359 | // Post output data. |
| 360 | PostData(1, fList); |
| 361 | PostData(2, fListBF); |
| 362 | if(fRunShuffling) PostData(3, fListBFS); |
| 363 | if(fUsePID) PostData(4, fHistListPIDQA); //PID |
| 364 | } |
| 365 | |
| 366 | //________________________________________________________________________ |
| 367 | void AliAnalysisTaskBF::UserExec(Option_t *) { |
| 368 | // Main loop |
| 369 | // Called for each event |
| 370 | TString gAnalysisLevel = fBalance->GetAnalysisLevel(); |
| 371 | |
| 9fd4b54e |
372 | AliESDtrack *trackTPC = NULL; |
| b8057da6 |
373 | |
| 374 | Int_t gNumberOfAcceptedTracks = 0; |
| 375 | Float_t fCentrality = 0.; |
| 376 | |
| 377 | // for HBT like cuts need magnetic field sign |
| 378 | Float_t bSign = 0; // only used in AOD so far |
| 379 | |
| 380 | // vector holding the charges/kinematics of all tracks (charge,y,eta,phi,p0,p1,p2,pt,E) |
| 381 | vector<Double_t> *chargeVectorShuffle[9]; // this will be shuffled |
| 382 | vector<Double_t> *chargeVector[9]; // original charge |
| 383 | for(Int_t i = 0; i < 9; i++){ |
| 384 | chargeVectorShuffle[i] = new vector<Double_t>; |
| 385 | chargeVector[i] = new vector<Double_t>; |
| 386 | } |
| 387 | |
| 9fd4b54e |
388 | Double_t vCharge; |
| 389 | Double_t vY; |
| 390 | Double_t vEta; |
| 391 | Double_t vPhi; |
| 392 | Double_t vP[3]; |
| 393 | Double_t vPt; |
| 394 | Double_t vE; |
| b8057da6 |
395 | |
| 396 | if(fUsePID) { |
| 397 | fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse(); |
| 398 | if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler"); |
| 399 | } |
| 400 | |
| 401 | //ESD analysis |
| 402 | if(gAnalysisLevel == "ESD") { |
| 403 | AliESDEvent* gESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE |
| 404 | if (!gESD) { |
| 405 | Printf("ERROR: gESD not available"); |
| 406 | return; |
| 407 | } |
| 408 | |
| 409 | // store offline trigger bits |
| 410 | fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()); |
| 411 | |
| 412 | // event selection done in AliAnalysisTaskSE::Exec() --> this is not used |
| 413 | fHistEventStats->Fill(1); //all events |
| 414 | Bool_t isSelected = kTRUE; |
| 415 | if(fUseOfflineTrigger) |
| 416 | isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected(); |
| 417 | if(isSelected) { |
| 418 | fHistEventStats->Fill(2); //triggered events |
| 419 | |
| 420 | if(fUseCentrality) { |
| 421 | //Centrality stuff |
| 422 | AliCentrality *centrality = gESD->GetCentrality(); |
| 423 | |
| 424 | fCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data()); |
| 425 | |
| 426 | // take only events inside centrality class |
| 427 | if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin, |
| 428 | fCentralityPercentileMax, |
| 429 | fCentralityEstimator.Data())) |
| 430 | return; |
| 431 | |
| 432 | // centrality QA (V0M) |
| 433 | fHistV0M->Fill(gESD->GetVZEROData()->GetMTotV0A(), gESD->GetVZEROData()->GetMTotV0C()); |
| 434 | } |
| 435 | |
| 436 | const AliESDVertex *vertex = gESD->GetPrimaryVertex(); |
| 437 | if(vertex) { |
| 438 | if(vertex->GetNContributors() > 0) { |
| 439 | if(vertex->GetZRes() != 0) { |
| 440 | fHistEventStats->Fill(3); //events with a proper vertex |
| 441 | if(TMath::Abs(vertex->GetXv()) < fVxMax) { |
| 442 | if(TMath::Abs(vertex->GetYv()) < fVyMax) { |
| 443 | if(TMath::Abs(vertex->GetZv()) < fVzMax) { |
| 444 | fHistEventStats->Fill(4); //analayzed events |
| 445 | fHistVx->Fill(vertex->GetXv()); |
| 446 | fHistVy->Fill(vertex->GetYv()); |
| 447 | fHistVz->Fill(vertex->GetZv()); |
| 448 | |
| 449 | //Printf("There are %d tracks in this event", gESD->GetNumberOfTracks()); |
| 450 | for (Int_t iTracks = 0; iTracks < gESD->GetNumberOfTracks(); iTracks++) { |
| 451 | AliESDtrack* track = dynamic_cast<AliESDtrack *>(gESD->GetTrack(iTracks)); |
| 452 | if (!track) { |
| 453 | Printf("ERROR: Could not receive track %d", iTracks); |
| 454 | continue; |
| 455 | } |
| 456 | |
| 457 | // take only TPC only tracks |
| 9fd4b54e |
458 | trackTPC = new AliESDtrack(); |
| 459 | if(!track->FillTPCOnlyTrack(*trackTPC)) continue; |
| b8057da6 |
460 | |
| 461 | //ESD track cuts |
| 462 | if(fESDtrackCuts) |
| 9fd4b54e |
463 | if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue; |
| b8057da6 |
464 | |
| 465 | // fill QA histograms |
| 466 | Float_t b[2]; |
| 467 | Float_t bCov[3]; |
| 9fd4b54e |
468 | trackTPC->GetImpactParameters(b,bCov); |
| b8057da6 |
469 | if (bCov[0]<=0 || bCov[2]<=0) { |
| 470 | AliDebug(1, "Estimated b resolution lower or equal zero!"); |
| 471 | bCov[0]=0; bCov[2]=0; |
| 472 | } |
| 473 | |
| 474 | Int_t nClustersTPC = -1; |
| 9fd4b54e |
475 | nClustersTPC = trackTPC->GetTPCNclsIter1(); // TPC standalone |
| b8057da6 |
476 | //nClustersTPC = track->GetTPCclusters(0); // global track |
| 477 | Float_t chi2PerClusterTPC = -1; |
| 478 | if (nClustersTPC!=0) { |
| 9fd4b54e |
479 | chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone |
| b8057da6 |
480 | //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track |
| 481 | } |
| 482 | |
| 483 | //===========================PID===============================// |
| 484 | if(fUsePID) { |
| 485 | Double_t prob[AliPID::kSPECIES]={0.}; |
| 486 | Double_t probTPC[AliPID::kSPECIES]={0.}; |
| 487 | Double_t probTOF[AliPID::kSPECIES]={0.}; |
| 488 | Double_t probTPCTOF[AliPID::kSPECIES]={0.}; |
| 489 | |
| 490 | Double_t nSigma = 0.; |
| 491 | UInt_t detUsedTPC = 0; |
| 492 | UInt_t detUsedTOF = 0; |
| 493 | UInt_t detUsedTPCTOF = 0; |
| 494 | |
| 495 | //Decide what detector configuration we want to use |
| 496 | switch(fPidDetectorConfig) { |
| 497 | case kTPCpid: |
| 498 | fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC); |
| 499 | nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest)); |
| 500 | detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC); |
| 501 | for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) |
| 502 | prob[iSpecies] = probTPC[iSpecies]; |
| 503 | break; |
| 504 | case kTOFpid: |
| 505 | fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF); |
| 506 | nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest)); |
| 507 | detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF); |
| 508 | for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) |
| 509 | prob[iSpecies] = probTOF[iSpecies]; |
| 510 | break; |
| 511 | case kTPCTOF: |
| 512 | fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC); |
| 513 | detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF); |
| 514 | for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) |
| 515 | prob[iSpecies] = probTPCTOF[iSpecies]; |
| 516 | break; |
| 517 | default: |
| 518 | break; |
| 519 | }//end switch: define detector mask |
| 520 | |
| 521 | //Filling the PID QA |
| 522 | Double_t tofTime = -999., length = 999., tof = -999.; |
| 523 | Double_t c = TMath::C()*1.E-9;// m/ns |
| 524 | Double_t beta = -999.; |
| 525 | Double_t nSigmaTOFForParticleOfInterest = -999.; |
| 526 | if ( (track->IsOn(AliESDtrack::kTOFin)) && |
| 527 | (track->IsOn(AliESDtrack::kTIME)) ) { |
| 528 | tofTime = track->GetTOFsignal();//in ps |
| 529 | length = track->GetIntegratedLength(); |
| 530 | tof = tofTime*1E-3; // ns |
| 531 | |
| 532 | if (tof <= 0) { |
| 533 | //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n"); |
| 534 | continue; |
| 535 | } |
| 536 | if (length <= 0){ |
| 537 | //printf("WARNING: track with negative length found!Skipping this track for PID checks\n"); |
| 538 | continue; |
| 539 | } |
| 540 | |
| 541 | length = length*0.01; // in meters |
| 542 | tof = tof*c; |
| 543 | beta = length/tof; |
| 544 | |
| 545 | nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest); |
| 546 | fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta); |
| 547 | fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]); |
| 548 | fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest); |
| 549 | }//TOF signal |
| 550 | |
| 551 | |
| 552 | Double_t nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest); |
| 553 | fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal()); |
| 554 | fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); |
| 555 | fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); |
| 556 | fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]); |
| 557 | //end of QA-before pid |
| 558 | |
| 559 | if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) { |
| 560 | //Make the decision based on the n-sigma |
| 561 | if(fUsePIDnSigma) { |
| 562 | if(nSigma > fPIDNSigma) continue;} |
| 563 | |
| 564 | //Make the decision based on the bayesian |
| 565 | else if(fUsePIDPropabilities) { |
| 566 | if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue; |
| 567 | if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue; |
| 568 | } |
| 569 | |
| 570 | //Fill QA after the PID |
| 571 | fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta); |
| 572 | fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]); |
| 573 | fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest); |
| 574 | |
| 575 | fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal()); |
| 576 | fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); |
| 577 | fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]); |
| 578 | fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); |
| 579 | } |
| 580 | |
| 581 | PostData(4, fHistListPIDQA); |
| 582 | } |
| 583 | //===========================PID===============================// |
| 9fd4b54e |
584 | vCharge = trackTPC->Charge(); |
| 585 | vY = trackTPC->Y(); |
| 586 | vEta = trackTPC->Eta(); |
| 587 | vPhi = trackTPC->Phi() * TMath::RadToDeg(); |
| 588 | vE = trackTPC->E(); |
| 589 | vPt = trackTPC->Pt(); |
| 590 | trackTPC->PxPyPz(vP); |
| 591 | fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC); |
| b8057da6 |
592 | fHistDCA->Fill(b[1],b[0]); |
| 593 | fHistChi2->Fill(chi2PerClusterTPC); |
| 9fd4b54e |
594 | fHistPt->Fill(vPt); |
| 595 | fHistEta->Fill(vEta); |
| 596 | fHistPhi->Fill(vPhi); |
| 597 | fHistRapidity->Fill(vY); |
| 598 | if(vCharge > 0) fHistPhiPos->Fill(vPhi); |
| 599 | else if(vCharge < 0) fHistPhiNeg->Fill(vPhi); |
| b8057da6 |
600 | |
| 601 | // fill charge vector |
| 9fd4b54e |
602 | chargeVector[0]->push_back(vCharge); |
| 603 | chargeVector[1]->push_back(vY); |
| 604 | chargeVector[2]->push_back(vEta); |
| 605 | chargeVector[3]->push_back(vPhi); |
| 606 | chargeVector[4]->push_back(vP[0]); |
| 607 | chargeVector[5]->push_back(vP[1]); |
| 608 | chargeVector[6]->push_back(vP[2]); |
| 609 | chargeVector[7]->push_back(vPt); |
| 610 | chargeVector[8]->push_back(vE); |
| b8057da6 |
611 | |
| 612 | if(fRunShuffling) { |
| 9fd4b54e |
613 | chargeVectorShuffle[0]->push_back(vCharge); |
| 614 | chargeVectorShuffle[1]->push_back(vY); |
| 615 | chargeVectorShuffle[2]->push_back(vEta); |
| 616 | chargeVectorShuffle[3]->push_back(vPhi); |
| 617 | chargeVectorShuffle[4]->push_back(vP[0]); |
| 618 | chargeVectorShuffle[5]->push_back(vP[1]); |
| 619 | chargeVectorShuffle[6]->push_back(vP[2]); |
| 620 | chargeVectorShuffle[7]->push_back(vPt); |
| 621 | chargeVectorShuffle[8]->push_back(vE); |
| b8057da6 |
622 | } |
| 623 | |
| 9fd4b54e |
624 | delete trackTPC; |
| b8057da6 |
625 | |
| 626 | } //track loop |
| 627 | }//Vz cut |
| 628 | }//Vy cut |
| 629 | }//Vx cut |
| 630 | }//proper vertex resolution |
| 631 | }//proper number of contributors |
| 632 | }//vertex object valid |
| 633 | }//triggered event |
| 634 | }//ESD analysis |
| 635 | |
| 636 | //AOD analysis (vertex and track cuts also here!!!!) |
| 637 | else if(gAnalysisLevel == "AOD") { |
| 638 | AliAODEvent* gAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // from TaskSE |
| 639 | if(!gAOD) { |
| 640 | Printf("ERROR: gAOD not available"); |
| 641 | return; |
| 642 | } |
| 643 | |
| 644 | // for HBT like cuts need magnetic field sign |
| 645 | bSign = (gAOD->GetMagneticField() > 0) ? 1 : -1; |
| 646 | |
| 647 | AliAODHeader *aodHeader = gAOD->GetHeader(); |
| 648 | |
| 649 | // store offline trigger bits |
| 650 | fHistTriggerStats->Fill(aodHeader->GetOfflineTrigger()); |
| 651 | |
| 652 | // event selection done in AliAnalysisTaskSE::Exec() --> this is not used |
| 653 | fHistEventStats->Fill(1); //all events |
| 654 | Bool_t isSelected = kTRUE; |
| 655 | if(fUseOfflineTrigger) |
| 656 | isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected(); |
| 657 | if(isSelected) { |
| 658 | fHistEventStats->Fill(2); //triggered events |
| 659 | |
| 660 | //Centrality stuff (centrality in AOD header) |
| 661 | if(fUseCentrality) { |
| 662 | fCentrality = aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data()); |
| 15de2c8a |
663 | |
| 664 | // in OLD AODs (i.e. AOD049) fCentrality can be == 0 |
| 665 | if(fCentrality == 0) |
| 666 | return; |
| b8057da6 |
667 | |
| 668 | // QA for centrality estimators |
| 669 | fHistCentStats->Fill(0.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0M")); |
| 670 | fHistCentStats->Fill(1.,aodHeader->GetCentralityP()->GetCentralityPercentile("FMD")); |
| 671 | fHistCentStats->Fill(2.,aodHeader->GetCentralityP()->GetCentralityPercentile("TRK")); |
| 672 | fHistCentStats->Fill(3.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKL")); |
| 673 | fHistCentStats->Fill(4.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL0")); |
| 674 | fHistCentStats->Fill(5.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL1")); |
| 675 | fHistCentStats->Fill(6.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0MvsFMD")); |
| 676 | fHistCentStats->Fill(7.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKLvsV0M")); |
| 677 | fHistCentStats->Fill(8.,aodHeader->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC")); |
| 678 | |
| 679 | // take only events inside centrality class |
| 680 | if((fCentrality < fCentralityPercentileMin) || (fCentrality > fCentralityPercentileMax)) |
| 681 | return; |
| 682 | |
| 683 | // centrality QA (V0M) |
| 684 | fHistV0M->Fill(gAOD->GetVZEROData()->GetMTotV0A(), gAOD->GetVZEROData()->GetMTotV0C()); |
| 685 | |
| 686 | // centrality QA (reference tracks) |
| 687 | fHistRefTracks->Fill(0.,aodHeader->GetRefMultiplicity()); |
| 688 | fHistRefTracks->Fill(1.,aodHeader->GetRefMultiplicityPos()); |
| 689 | fHistRefTracks->Fill(2.,aodHeader->GetRefMultiplicityNeg()); |
| 690 | fHistRefTracks->Fill(3.,aodHeader->GetTPConlyRefMultiplicity()); |
| 691 | fHistRefTracks->Fill(4.,aodHeader->GetNumberOfITSClusters(0)); |
| 692 | fHistRefTracks->Fill(5.,aodHeader->GetNumberOfITSClusters(1)); |
| 693 | fHistRefTracks->Fill(6.,aodHeader->GetNumberOfITSClusters(2)); |
| 694 | fHistRefTracks->Fill(7.,aodHeader->GetNumberOfITSClusters(3)); |
| 695 | fHistRefTracks->Fill(8.,aodHeader->GetNumberOfITSClusters(4)); |
| 696 | } |
| 697 | |
| 698 | const AliAODVertex *vertex = gAOD->GetPrimaryVertex(); |
| 699 | |
| 700 | if(vertex) { |
| 701 | Double32_t fCov[6]; |
| 702 | vertex->GetCovarianceMatrix(fCov); |
| 703 | |
| 704 | if(vertex->GetNContributors() > 0) { |
| 705 | if(fCov[5] != 0) { |
| 706 | fHistEventStats->Fill(3); //events with a proper vertex |
| 707 | if(TMath::Abs(vertex->GetX()) < fVxMax) { |
| 708 | if(TMath::Abs(vertex->GetY()) < fVyMax) { |
| 709 | if(TMath::Abs(vertex->GetZ()) < fVzMax) { |
| 710 | fHistEventStats->Fill(4); //analyzed events |
| 711 | fHistVx->Fill(vertex->GetX()); |
| 712 | fHistVy->Fill(vertex->GetY()); |
| 713 | fHistVz->Fill(vertex->GetZ()); |
| 714 | |
| 715 | //===========================================// |
| 716 | TExMap *trackMap = new TExMap(); |
| 717 | for (Int_t iTracks = 0; iTracks < gAOD->GetNumberOfTracks(); iTracks++) { |
| 718 | AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(gAOD->GetTrack(iTracks)); |
| 719 | if (!aodTrack) { |
| 720 | Printf("ERROR: Could not receive track %d", iTracks); |
| 721 | continue; |
| 722 | } |
| 723 | Int_t gID = aodTrack->GetID(); |
| 9fd4b54e |
724 | //if (!aodTrack->TestFilterBit(fAODtrackCutBit)) trackMap->Add(gID, iTracks); |
| e1d6ee8e |
725 | if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, iTracks); |
| b8057da6 |
726 | } |
| 727 | AliAODTrack* newAodTrack; |
| 728 | //===========================================// |
| 729 | |
| 730 | //Printf("There are %d tracks in this event", gAOD->GetNumberOfTracks()); |
| 731 | for (Int_t iTracks = 0; iTracks < gAOD->GetNumberOfTracks(); iTracks++) { |
| 732 | AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(gAOD->GetTrack(iTracks)); |
| 733 | if (!aodTrack) { |
| 734 | Printf("ERROR: Could not receive track %d", iTracks); |
| 735 | continue; |
| 736 | } |
| 737 | |
| e1d6ee8e |
738 | // AOD track cuts |
| b8057da6 |
739 | // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C |
| b8057da6 |
740 | //===========================================// |
| 741 | // take only TPC only tracks |
| 742 | fHistTrackStats->Fill(aodTrack->GetFilterMap()); |
| 9fd4b54e |
743 | if(!aodTrack->TestFilterBit(fAODtrackCutBit)) continue; |
| e1d6ee8e |
744 | |
| b8057da6 |
745 | Int_t gID = aodTrack->GetID(); |
| 746 | newAodTrack = gID >= 0 ? aodTrack : gAOD->GetTrack(trackMap->GetValue(-1-gID)); |
| e1d6ee8e |
747 | //Printf("Label: %d - Pt: %lf (old) - %d - Pt: %lf(new)",gID,aodTrack->Pt(), newAodTrack->GetID(), newAodTrack->Pt()); |
| b8057da6 |
748 | //===========================================// |
| 749 | |
| 750 | //fHistTrackStats->Fill(aodTrack->GetFilterMap()); |
| 9fd4b54e |
751 | //if(!aodTrack->TestFilterBit(fAODtrackCutBit)) continue; |
| b8057da6 |
752 | |
| 9fd4b54e |
753 | vCharge = aodTrack->Charge(); |
| 754 | vY = aodTrack->Y(); |
| 755 | vEta = aodTrack->Eta(); |
| 756 | vPhi = aodTrack->Phi() * TMath::RadToDeg(); |
| 757 | vE = aodTrack->E(); |
| 758 | vPt = aodTrack->Pt(); |
| 759 | aodTrack->PxPyPz(vP); |
| b8057da6 |
760 | |
| 9fd4b54e |
761 | Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on) |
| 762 | Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on) |
| b8057da6 |
763 | |
| 764 | |
| 765 | // Kinematics cuts from ESD track cuts |
| 9fd4b54e |
766 | if( vPt < fPtMin || vPt > fPtMax) continue; |
| b8057da6 |
767 | |
| 768 | if (!fUsePID) { |
| 9fd4b54e |
769 | if( vEta < fEtaMin || vEta > fEtaMax) continue; |
| b8057da6 |
770 | } |
| 771 | |
| 772 | else if (fUsePID){ |
| 9fd4b54e |
773 | if( vY < fEtaMin || vY > fEtaMax) continue; |
| b8057da6 |
774 | } |
| 775 | |
| 776 | // Extra DCA cuts (for systematic studies [!= -1]) |
| 777 | if( fDCAxyCut != -1 && fDCAzCut != -1){ |
| 9fd4b54e |
778 | if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){ |
| b8057da6 |
779 | continue; // 2D cut |
| 780 | } |
| 781 | } |
| 782 | |
| 783 | // Extra TPC cuts (for systematic studies [!= -1]) |
| 784 | if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){ |
| 785 | continue; |
| 786 | } |
| 787 | if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){ |
| 788 | continue; |
| 789 | } |
| 790 | |
| 791 | //===============================================PID==================================// |
| 792 | if(fUsePID) { |
| 793 | Double_t prob[AliPID::kSPECIES]={0.}; |
| 794 | Double_t probTPC[AliPID::kSPECIES]={0.}; |
| 795 | Double_t probTOF[AliPID::kSPECIES]={0.}; |
| 796 | Double_t probTPCTOF[AliPID::kSPECIES]={0.}; |
| 797 | |
| 798 | Double_t nSigma = 0.; |
| 799 | UInt_t detUsedTPC = 0; |
| 800 | UInt_t detUsedTOF = 0; |
| 801 | UInt_t detUsedTPCTOF = 0; |
| 802 | |
| 803 | //Decide what detector configuration we want to use |
| 804 | switch(fPidDetectorConfig) { |
| 805 | case kTPCpid: |
| 806 | fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC); |
| 807 | nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)fParticleOfInterest)); |
| 808 | //detUsedTPC = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPC); |
| 809 | detUsedTPC = (AliPIDResponse::kDetTPC); |
| 810 | for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) |
| 811 | prob[iSpecies] = probTPC[iSpecies]; |
| 812 | break; |
| 813 | case kTOFpid: |
| 814 | fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF); |
| 815 | nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(newAodTrack,(AliPID::EParticleType)fParticleOfInterest)); |
| 816 | //detUsedTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTOF); |
| 817 | detUsedTPC = (AliPIDResponse::kDetTPC); |
| 818 | for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) |
| 819 | prob[iSpecies] = probTOF[iSpecies]; |
| 820 | break; |
| 821 | case kTPCTOF: |
| 822 | fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC); |
| 823 | //detUsedTPCTOF = fPIDCombined->ComputeProbabilities(newAodTrack, fPIDResponse, probTPCTOF); |
| 824 | detUsedTPC = (AliPIDResponse::kDetTPC); |
| 825 | for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) |
| 826 | prob[iSpecies] = probTPCTOF[iSpecies]; |
| 827 | break; |
| 828 | default: |
| 829 | break; |
| 830 | }//end switch: define detector mask |
| 831 | |
| 832 | //Filling the PID QA |
| 833 | Double_t tofTime = -999., tof = -999.; //length = 999., tof = -999.; |
| 834 | //Double_t c = TMath::C()*1.E-9;// m/ns |
| 835 | //Double_t beta = -999.; |
| 836 | Double_t nSigmaTOFForParticleOfInterest = -999.; |
| 837 | if ( (newAodTrack->IsOn(AliAODTrack::kTOFin)) && |
| 838 | (newAodTrack->IsOn(AliAODTrack::kTIME)) ) { |
| 839 | tofTime = newAodTrack->GetTOFsignal();//in ps |
| 840 | //length = newAodTrack->GetIntegratedLength(); |
| 841 | tof = tofTime*1E-3; // ns |
| 842 | |
| 843 | if (tof <= 0) { |
| 844 | //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n"); |
| 845 | continue; |
| 846 | } |
| 847 | //if (length <= 0){ |
| 848 | //printf("WARNING: track with negative length found!Skipping this track for PID checks\n"); |
| 849 | //continue; |
| 850 | //} |
| 851 | |
| 852 | //length = length*0.01; // in meters |
| 853 | //tof = tof*c; |
| 854 | //beta = length/tof; |
| 855 | |
| 856 | nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(newAodTrack,(AliPID::EParticleType)fParticleOfInterest); |
| 857 | //fHistBetavsPTOFbeforePID ->Fill(aodTrack->P()*aodTrack->Charge(),beta); |
| 858 | fHistProbTOFvsPtbeforePID ->Fill(newAodTrack->Pt(),probTOF[fParticleOfInterest]); |
| 859 | fHistNSigmaTOFvsPtbeforePID ->Fill(newAodTrack->Pt(),nSigmaTOFForParticleOfInterest); |
| 860 | }//TOF signal |
| 861 | |
| 862 | |
| 863 | Double_t nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)fParticleOfInterest); |
| 864 | fHistdEdxVsPTPCbeforePID -> Fill(newAodTrack->P()*newAodTrack->Charge(),newAodTrack->GetTPCsignal()); |
| 865 | fHistProbTPCvsPtbeforePID -> Fill(newAodTrack->Pt(),probTPC[fParticleOfInterest]); |
| 866 | fHistNSigmaTPCvsPtbeforePID -> Fill(newAodTrack->Pt(),nSigmaTPCForParticleOfInterest); |
| 867 | fHistProbTPCTOFvsPtbeforePID -> Fill(newAodTrack->Pt(),probTPCTOF[fParticleOfInterest]); |
| 868 | //end of QA-before pid |
| 869 | |
| 870 | if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) { |
| 871 | //Make the decision based on the n-sigma |
| 872 | if(fUsePIDnSigma) { |
| 873 | if(nSigma > fPIDNSigma) continue;} |
| 874 | |
| 875 | //Make the decision based on the bayesian |
| 876 | else if(fUsePIDPropabilities) { |
| 877 | if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue; |
| 878 | if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue; |
| 879 | } |
| 880 | |
| 881 | //Fill QA after the PID |
| 882 | //fHistBetavsPTOFafterPID ->Fill(newAodTrack->P()*newAodTrack->Charge(),beta); |
| 883 | fHistProbTOFvsPtafterPID ->Fill(newAodTrack->Pt(),probTOF[fParticleOfInterest]); |
| 884 | fHistNSigmaTOFvsPtafterPID ->Fill(newAodTrack->Pt(),nSigmaTOFForParticleOfInterest); |
| 885 | |
| 886 | fHistdEdxVsPTPCafterPID -> Fill(newAodTrack->P()*newAodTrack->Charge(),newAodTrack->GetTPCsignal()); |
| 887 | fHistProbTPCvsPtafterPID -> Fill(newAodTrack->Pt(),probTPC[fParticleOfInterest]); |
| 888 | fHistProbTPCTOFvsPtafterPID -> Fill(newAodTrack->Pt(),probTPCTOF[fParticleOfInterest]); |
| 889 | fHistNSigmaTPCvsPtafterPID -> Fill(newAodTrack->Pt(),nSigmaTPCForParticleOfInterest); |
| 890 | } |
| 891 | |
| 892 | PostData(4, fHistListPIDQA); |
| 893 | } |
| 894 | |
| 895 | //=========================================================PID=================================================================// |
| 896 | |
| 897 | // fill QA histograms |
| 898 | fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls()); |
| 9fd4b54e |
899 | fHistDCA->Fill(dcaZ,dcaXY); |
| b8057da6 |
900 | fHistChi2->Fill(aodTrack->Chi2perNDF()); |
| 9fd4b54e |
901 | fHistPt->Fill(vPt); |
| 902 | fHistEta->Fill(vEta); |
| 903 | fHistPhi->Fill(vPhi); |
| 904 | fHistRapidity->Fill(vY); |
| 905 | if(vCharge > 0) fHistPhiPos->Fill(vPhi); |
| 906 | else if(vCharge < 0) fHistPhiNeg->Fill(vPhi); |
| b8057da6 |
907 | |
| 908 | // fill charge vector |
| 9fd4b54e |
909 | chargeVector[0]->push_back(vCharge); |
| 910 | chargeVector[1]->push_back(vY); |
| 911 | chargeVector[2]->push_back(vEta); |
| 912 | chargeVector[3]->push_back(vPhi); |
| 913 | chargeVector[4]->push_back(vP[0]); |
| 914 | chargeVector[5]->push_back(vP[1]); |
| 915 | chargeVector[6]->push_back(vP[2]); |
| 916 | chargeVector[7]->push_back(vPt); |
| 917 | chargeVector[8]->push_back(vE); |
| b8057da6 |
918 | |
| 919 | if(fRunShuffling) { |
| 9fd4b54e |
920 | chargeVectorShuffle[0]->push_back(vCharge); |
| 921 | chargeVectorShuffle[1]->push_back(vY); |
| 922 | chargeVectorShuffle[2]->push_back(vEta); |
| 923 | chargeVectorShuffle[3]->push_back(vPhi); |
| 924 | chargeVectorShuffle[4]->push_back(vP[0]); |
| 925 | chargeVectorShuffle[5]->push_back(vP[1]); |
| 926 | chargeVectorShuffle[6]->push_back(vP[2]); |
| 927 | chargeVectorShuffle[7]->push_back(vPt); |
| 928 | chargeVectorShuffle[8]->push_back(vE); |
| b8057da6 |
929 | } |
| 930 | |
| 931 | gNumberOfAcceptedTracks += 1; |
| 932 | |
| 933 | } //track loop |
| 934 | }//Vz cut |
| 935 | }//Vy cut |
| 936 | }//Vx cut |
| 937 | }//proper vertex resolution |
| 938 | }//proper number of contributors |
| 939 | }//vertex object valid |
| 940 | }//triggered event |
| 941 | }//AOD analysis |
| 942 | |
| 943 | //MC-ESD analysis |
| 944 | if(gAnalysisLevel == "MCESD") { |
| 945 | AliMCEvent* mcEvent = MCEvent(); |
| 946 | if (!mcEvent) { |
| 947 | Printf("ERROR: mcEvent not available"); |
| 948 | return; |
| 949 | } |
| 950 | |
| 951 | AliESDEvent* gESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE |
| 952 | if (!gESD) { |
| 953 | Printf("ERROR: gESD not available"); |
| 954 | return; |
| 955 | } |
| 956 | |
| 957 | // store offline trigger bits |
| 958 | fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()); |
| 959 | |
| 960 | // event selection done in AliAnalysisTaskSE::Exec() --> this is not used |
| 961 | fHistEventStats->Fill(1); //all events |
| 962 | Bool_t isSelected = kTRUE; |
| 963 | if(fUseOfflineTrigger) |
| 964 | isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected(); |
| 965 | if(isSelected) { |
| 966 | fHistEventStats->Fill(2); //triggered events |
| 967 | |
| 968 | if(fUseCentrality) { |
| 969 | //Centrality stuff |
| 970 | AliCentrality *centrality = gESD->GetCentrality(); |
| 971 | |
| 972 | fCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data()); |
| 973 | |
| 974 | // take only events inside centrality class |
| 975 | if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin, |
| 976 | fCentralityPercentileMax, |
| 977 | fCentralityEstimator.Data())) |
| 978 | return; |
| 979 | |
| 980 | // centrality QA (V0M) |
| 981 | fHistV0M->Fill(gESD->GetVZEROData()->GetMTotV0A(), gESD->GetVZEROData()->GetMTotV0C()); |
| 982 | } |
| 983 | |
| 984 | const AliESDVertex *vertex = gESD->GetPrimaryVertex(); |
| 985 | if(vertex) { |
| 986 | if(vertex->GetNContributors() > 0) { |
| 987 | if(vertex->GetZRes() != 0) { |
| 988 | fHistEventStats->Fill(3); //events with a proper vertex |
| 989 | if(TMath::Abs(vertex->GetXv()) < fVxMax) { |
| 990 | if(TMath::Abs(vertex->GetYv()) < fVyMax) { |
| 991 | if(TMath::Abs(vertex->GetZv()) < fVzMax) { |
| 992 | fHistEventStats->Fill(4); //analayzed events |
| 993 | fHistVx->Fill(vertex->GetXv()); |
| 994 | fHistVy->Fill(vertex->GetYv()); |
| 995 | fHistVz->Fill(vertex->GetZv()); |
| 996 | |
| 997 | //Printf("There are %d tracks in this event", gESD->GetNumberOfTracks()); |
| 998 | for (Int_t iTracks = 0; iTracks < gESD->GetNumberOfTracks(); iTracks++) { |
| 999 | AliESDtrack* track = dynamic_cast<AliESDtrack *>(gESD->GetTrack(iTracks)); |
| 1000 | if (!track) { |
| 1001 | Printf("ERROR: Could not receive track %d", iTracks); |
| 1002 | continue; |
| 1003 | } |
| 1004 | |
| 1005 | Int_t label = TMath::Abs(track->GetLabel()); |
| 1006 | if(label > mcEvent->GetNumberOfTracks()) continue; |
| 1007 | if(label > mcEvent->GetNumberOfPrimaries()) continue; |
| 1008 | |
| 1009 | AliMCParticle* mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label)); |
| 1010 | if(!mcTrack) continue; |
| 1011 | |
| 1012 | // take only TPC only tracks |
| 9fd4b54e |
1013 | trackTPC = new AliESDtrack(); |
| 1014 | if(!track->FillTPCOnlyTrack(*trackTPC)) continue; |
| b8057da6 |
1015 | |
| 1016 | //ESD track cuts |
| 1017 | if(fESDtrackCuts) |
| 9fd4b54e |
1018 | if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue; |
| b8057da6 |
1019 | |
| 1020 | // fill QA histograms |
| 1021 | Float_t b[2]; |
| 1022 | Float_t bCov[3]; |
| 9fd4b54e |
1023 | trackTPC->GetImpactParameters(b,bCov); |
| b8057da6 |
1024 | if (bCov[0]<=0 || bCov[2]<=0) { |
| 1025 | AliDebug(1, "Estimated b resolution lower or equal zero!"); |
| 1026 | bCov[0]=0; bCov[2]=0; |
| 1027 | } |
| 1028 | |
| 1029 | Int_t nClustersTPC = -1; |
| 9fd4b54e |
1030 | nClustersTPC = trackTPC->GetTPCNclsIter1(); // TPC standalone |
| b8057da6 |
1031 | //nClustersTPC = track->GetTPCclusters(0); // global track |
| 1032 | Float_t chi2PerClusterTPC = -1; |
| 1033 | if (nClustersTPC!=0) { |
| 9fd4b54e |
1034 | chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC); // TPC standalone |
| b8057da6 |
1035 | //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC); // global track |
| 1036 | } |
| 1037 | |
| 9fd4b54e |
1038 | vCharge = trackTPC->Charge(); |
| 1039 | vY = trackTPC->Y(); |
| 1040 | vEta = trackTPC->Eta(); |
| 1041 | vPhi = trackTPC->Phi() * TMath::RadToDeg(); |
| 1042 | vE = trackTPC->E(); |
| 1043 | vPt = trackTPC->Pt(); |
| 1044 | trackTPC->PxPyPz(vP); |
| 1045 | |
| 1046 | fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC); |
| b8057da6 |
1047 | fHistDCA->Fill(b[1],b[0]); |
| 1048 | fHistChi2->Fill(chi2PerClusterTPC); |
| 9fd4b54e |
1049 | fHistPt->Fill(vPt); |
| 1050 | fHistEta->Fill(vEta); |
| 1051 | fHistPhi->Fill(vPhi); |
| 1052 | fHistRapidity->Fill(vY); |
| 1053 | if(vCharge > 0) fHistPhiPos->Fill(vPhi); |
| 1054 | else if(vCharge < 0) fHistPhiNeg->Fill(vPhi); |
| b8057da6 |
1055 | |
| 1056 | // fill charge vector |
| 9fd4b54e |
1057 | chargeVector[0]->push_back(vCharge); |
| 1058 | chargeVector[1]->push_back(vY); |
| 1059 | chargeVector[2]->push_back(vEta); |
| 1060 | chargeVector[3]->push_back(vPhi); |
| 1061 | chargeVector[4]->push_back(vP[0]); |
| 1062 | chargeVector[5]->push_back(vP[1]); |
| 1063 | chargeVector[6]->push_back(vP[2]); |
| 1064 | chargeVector[7]->push_back(vPt); |
| 1065 | chargeVector[8]->push_back(vE); |
| b8057da6 |
1066 | |
| 1067 | if(fRunShuffling) { |
| 9fd4b54e |
1068 | chargeVectorShuffle[0]->push_back(vCharge); |
| 1069 | chargeVectorShuffle[1]->push_back(vY); |
| 1070 | chargeVectorShuffle[2]->push_back(vEta); |
| 1071 | chargeVectorShuffle[3]->push_back(vPhi); |
| 1072 | chargeVectorShuffle[4]->push_back(vP[0]); |
| 1073 | chargeVectorShuffle[5]->push_back(vP[1]); |
| 1074 | chargeVectorShuffle[6]->push_back(vP[2]); |
| 1075 | chargeVectorShuffle[7]->push_back(vPt); |
| 1076 | chargeVectorShuffle[8]->push_back(vE); |
| b8057da6 |
1077 | } |
| 1078 | |
| 9fd4b54e |
1079 | delete trackTPC; |
| b8057da6 |
1080 | gNumberOfAcceptedTracks += 1; |
| 1081 | |
| 1082 | } //track loop |
| 1083 | }//Vz cut |
| 1084 | }//Vy cut |
| 1085 | }//Vx cut |
| 1086 | }//proper vertex resolution |
| 1087 | }//proper number of contributors |
| 1088 | }//vertex object valid |
| 1089 | }//triggered event |
| 1090 | }//MC-ESD analysis |
| 1091 | |
| 1092 | //MC analysis |
| 1093 | else if(gAnalysisLevel == "MC") { |
| 1094 | AliMCEvent* mcEvent = MCEvent(); |
| 1095 | if (!mcEvent) { |
| 1096 | Printf("ERROR: mcEvent not available"); |
| 1097 | return; |
| 1098 | } |
| 1099 | fHistEventStats->Fill(1); //total events |
| 1100 | fHistEventStats->Fill(2); //offline trigger |
| 1101 | |
| 1102 | Double_t gReactionPlane = 0., gImpactParameter = 0.; |
| 1103 | if(fUseCentrality) { |
| 1104 | //Get the MC header |
| 1105 | AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader()); |
| 1106 | if (headerH) { |
| 1107 | //Printf("====================================================="); |
| 1108 | //Printf("Reaction plane angle: %lf",headerH->ReactionPlaneAngle()); |
| 1109 | //Printf("Impact parameter: %lf",headerH->ImpactParameter()); |
| 1110 | //Printf("====================================================="); |
| 1111 | gReactionPlane = headerH->ReactionPlaneAngle(); |
| 1112 | gImpactParameter = headerH->ImpactParameter(); |
| 1113 | fCentrality = gImpactParameter; |
| 1114 | } |
| 1115 | fCentrality = gImpactParameter; |
| 1116 | |
| 1117 | // take only events inside centrality class (DIDN'T CHANGE THIS UP TO NOW) |
| 1118 | if((fImpactParameterMin > gImpactParameter) || (fImpactParameterMax < gImpactParameter)) |
| 1119 | return; |
| 1120 | } |
| 1121 | |
| 1122 | AliGenEventHeader *header = mcEvent->GenEventHeader(); |
| 1123 | if(!header) return; |
| 1124 | |
| 1125 | TArrayF gVertexArray; |
| 1126 | header->PrimaryVertex(gVertexArray); |
| 1127 | //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)", |
| 1128 | //gVertexArray.At(0), |
| 1129 | //gVertexArray.At(1), |
| 1130 | //gVertexArray.At(2)); |
| 1131 | fHistEventStats->Fill(3); //events with a proper vertex |
| 1132 | if(TMath::Abs(gVertexArray.At(0)) < fVxMax) { |
| 1133 | if(TMath::Abs(gVertexArray.At(1)) < fVyMax) { |
| 1134 | if(TMath::Abs(gVertexArray.At(2)) < fVzMax) { |
| 1135 | fHistEventStats->Fill(4); //analayzed events |
| 1136 | fHistVx->Fill(gVertexArray.At(0)); |
| 1137 | fHistVy->Fill(gVertexArray.At(1)); |
| 1138 | fHistVz->Fill(gVertexArray.At(2)); |
| 1139 | |
| 1140 | Printf("There are %d tracks in this event", mcEvent->GetNumberOfPrimaries()); |
| 1141 | for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfPrimaries(); iTracks++) { |
| 1142 | AliMCParticle* track = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(iTracks)); |
| 1143 | if (!track) { |
| 1144 | Printf("ERROR: Could not receive particle %d", iTracks); |
| 1145 | continue; |
| 1146 | } |
| 1147 | |
| 1148 | //exclude non stable particles |
| 1149 | if(!mcEvent->IsPhysicalPrimary(iTracks)) continue; |
| 1150 | |
| 9fd4b54e |
1151 | vEta = track->Eta(); |
| 1152 | vPt = track->Pt(); |
| 1153 | vY = track->Y(); |
| b8057da6 |
1154 | |
| 9fd4b54e |
1155 | if( vPt < fPtMin || vPt > fPtMax) |
| b8057da6 |
1156 | continue; |
| 1157 | if (!fUsePID) { |
| 9fd4b54e |
1158 | if( vEta < fEtaMin || vEta > fEtaMax) continue; |
| b8057da6 |
1159 | } |
| 1160 | else if (fUsePID){ |
| 9fd4b54e |
1161 | if( vY < fEtaMin || vY > fEtaMax) continue; |
| b8057da6 |
1162 | } |
| 1163 | |
| 1164 | //analyze one set of particles |
| 1165 | if(fUseMCPdgCode) { |
| 1166 | TParticle *particle = track->Particle(); |
| 1167 | if(!particle) continue; |
| 1168 | |
| 1169 | Int_t gPdgCode = particle->GetPdgCode(); |
| 1170 | if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode)) |
| 1171 | continue; |
| 1172 | } |
| 1173 | |
| 1174 | //Use the acceptance parameterization |
| 1175 | if(fAcceptanceParameterization) { |
| 1176 | Double_t gRandomNumber = gRandom->Rndm(); |
| 1177 | if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt())) |
| 1178 | continue; |
| 1179 | } |
| 1180 | |
| 1181 | //Exclude resonances |
| 1182 | if(fExcludeResonancesInMC) { |
| 1183 | TParticle *particle = track->Particle(); |
| 1184 | if(!particle) continue; |
| 1185 | |
| 1186 | Bool_t kExcludeParticle = kFALSE; |
| 1187 | Int_t gMotherIndex = particle->GetFirstMother(); |
| 1188 | if(gMotherIndex != -1) { |
| 1189 | AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(gMotherIndex)); |
| 1190 | if(motherTrack) { |
| 1191 | TParticle *motherParticle = motherTrack->Particle(); |
| 1192 | if(motherParticle) { |
| 1193 | Int_t pdgCodeOfMother = motherParticle->GetPdgCode(); |
| 1194 | //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) { |
| 1195 | if(pdgCodeOfMother == 113) { |
| 1196 | kExcludeParticle = kTRUE; |
| 1197 | } |
| 1198 | } |
| 1199 | } |
| 1200 | } |
| 1201 | |
| 1202 | //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi |
| 1203 | if(kExcludeParticle) continue; |
| 1204 | } |
| 1205 | |
| 9fd4b54e |
1206 | vCharge = track->Charge(); |
| 1207 | vPhi = track->Phi(); |
| 1208 | vE = track->E(); |
| 1209 | track->PxPyPz(vP); |
| 1210 | //Printf("phi (before): %lf",vPhi); |
| b8057da6 |
1211 | |
| 9fd4b54e |
1212 | fHistPt->Fill(vPt); |
| 1213 | fHistEta->Fill(vEta); |
| 1214 | fHistPhi->Fill(vPhi); |
| 1215 | fHistRapidity->Fill(vY); |
| 1216 | if(vCharge > 0) fHistPhiPos->Fill(vPhi); |
| 1217 | else if(vCharge < 0) fHistPhiNeg->Fill(vPhi); |
| b8057da6 |
1218 | |
| 1219 | //Flow after burner |
| 1220 | if(fUseFlowAfterBurner) { |
| 1221 | Double_t precisionPhi = 0.001; |
| 1222 | Int_t maxNumberOfIterations = 100; |
| 1223 | |
| 9fd4b54e |
1224 | Double_t phi0 = vPhi; |
| 1225 | Double_t gV2 = fDifferentialV2->Eval(vPt); |
| b8057da6 |
1226 | |
| 1227 | for (Int_t j = 0; j < maxNumberOfIterations; j++) { |
| 9fd4b54e |
1228 | Double_t phiprev = vPhi; |
| 1229 | Double_t fl = vPhi - phi0 + gV2*TMath::Sin(2.*(vPhi - gReactionPlane)); |
| 1230 | Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(vPhi - gReactionPlane)); |
| 1231 | vPhi -= fl/fp; |
| 1232 | if (TMath::AreEqualAbs(phiprev,vPhi,precisionPhi)) break; |
| b8057da6 |
1233 | } |
| 9fd4b54e |
1234 | //Printf("phi (after): %lf\n",vPhi); |
| 648f1a5a |
1235 | Double_t vDeltaphiBefore = phi0 - gReactionPlane; |
| 9fd4b54e |
1236 | if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi(); |
| 1237 | fHistPhiBefore->Fill(vDeltaphiBefore); |
| 1238 | |
| 1239 | Double_t vDeltaphiAfter = vPhi - gReactionPlane; |
| 1240 | if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi(); |
| 1241 | fHistPhiAfter->Fill(vDeltaphiAfter); |
| b8057da6 |
1242 | } |
| 1243 | |
| 9fd4b54e |
1244 | vPhi *= TMath::RadToDeg(); |
| b8057da6 |
1245 | |
| 1246 | // fill charge vector |
| 9fd4b54e |
1247 | chargeVector[0]->push_back(vCharge); |
| 1248 | chargeVector[1]->push_back(vY); |
| 1249 | chargeVector[2]->push_back(vEta); |
| 1250 | chargeVector[3]->push_back(vPhi); |
| 1251 | chargeVector[4]->push_back(vP[0]); |
| 1252 | chargeVector[5]->push_back(vP[1]); |
| 1253 | chargeVector[6]->push_back(vP[2]); |
| 1254 | chargeVector[7]->push_back(vPt); |
| 1255 | chargeVector[8]->push_back(vE); |
| b8057da6 |
1256 | |
| 1257 | if(fRunShuffling) { |
| 9fd4b54e |
1258 | chargeVectorShuffle[0]->push_back(vCharge); |
| 1259 | chargeVectorShuffle[1]->push_back(vY); |
| 1260 | chargeVectorShuffle[2]->push_back(vEta); |
| 1261 | chargeVectorShuffle[3]->push_back(vPhi); |
| 1262 | chargeVectorShuffle[4]->push_back(vP[0]); |
| 1263 | chargeVectorShuffle[5]->push_back(vP[1]); |
| 1264 | chargeVectorShuffle[6]->push_back(vP[2]); |
| 1265 | chargeVectorShuffle[7]->push_back(vPt); |
| 1266 | chargeVectorShuffle[8]->push_back(vE); |
| b8057da6 |
1267 | } |
| 1268 | gNumberOfAcceptedTracks += 1; |
| 1269 | |
| 1270 | } //track loop |
| 1271 | }//Vz cut |
| 1272 | }//Vy cut |
| 1273 | }//Vx cut |
| 1274 | }//MC analysis |
| 1275 | |
| 1276 | //multiplicity cut (used in pp) |
| 1277 | if(fUseMultiplicity) { |
| 1278 | if((gNumberOfAcceptedTracks < fNumberOfAcceptedTracksMin)||(gNumberOfAcceptedTracks > fNumberOfAcceptedTracksMax)) |
| 1279 | return; |
| 1280 | } |
| 1281 | fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks); |
| 1282 | |
| 1283 | // calculate balance function |
| 1284 | if(fUseMultiplicity) |
| 1285 | fBalance->CalculateBalance(gNumberOfAcceptedTracks,chargeVector,bSign); |
| 1286 | else |
| 1287 | fBalance->CalculateBalance(fCentrality,chargeVector,bSign); |
| 1288 | |
| 1289 | if(fRunShuffling) { |
| 1290 | // shuffle charges |
| 1291 | random_shuffle( chargeVectorShuffle[0]->begin(), chargeVectorShuffle[0]->end() ); |
| 1292 | if(fUseMultiplicity) |
| 1293 | fShuffledBalance->CalculateBalance(gNumberOfAcceptedTracks,chargeVectorShuffle,bSign); |
| 1294 | else |
| 1295 | fShuffledBalance->CalculateBalance(fCentrality,chargeVectorShuffle,bSign); |
| 1296 | } |
| 1297 | } |
| 1298 | |
| 1299 | //________________________________________________________________________ |
| 1300 | void AliAnalysisTaskBF::FinishTaskOutput(){ |
| 1301 | //Printf("END BF"); |
| 1302 | |
| 1303 | if (!fBalance) { |
| 1304 | Printf("ERROR: fBalance not available"); |
| 1305 | return; |
| 1306 | } |
| 1307 | if(fRunShuffling) { |
| 1308 | if (!fShuffledBalance) { |
| 1309 | Printf("ERROR: fShuffledBalance not available"); |
| 1310 | return; |
| 1311 | } |
| 1312 | } |
| 1313 | |
| 1314 | } |
| 1315 | |
| 1316 | //________________________________________________________________________ |
| 1317 | void AliAnalysisTaskBF::Terminate(Option_t *) { |
| 1318 | // Draw result to the screen |
| 1319 | // Called once at the end of the query |
| 1320 | |
| 1321 | // not implemented ... |
| 1322 | |
| 1323 | } |
git://git.uio.no / u / mrichter / AliRoot.git / blame