]>
Commit | Line | Data |
---|---|---|
0f67a57c | 1 | /* $Id$ */ |
2 | ||
3 | #include "AlidNdEtaCorrectionTask.h" | |
4 | ||
5 | #include <TCanvas.h> | |
6 | #include <TChain.h> | |
7 | #include <TFile.h> | |
8 | #include <TH1F.h> | |
54b096ef | 9 | #include <TH2F.h> |
1c15d51a | 10 | #include <TH3F.h> |
54b096ef | 11 | #include <TProfile.h> |
0f67a57c | 12 | #include <TParticle.h> |
1c15d51a | 13 | #include <TParticlePDG.h> |
c180f65d | 14 | #include <TDatabasePDG.h> |
81be4ee8 | 15 | #include <TRandom.h> |
0f67a57c | 16 | |
17 | #include <AliLog.h> | |
18 | #include <AliESDVertex.h> | |
19 | #include <AliESDEvent.h> | |
81be4ee8 | 20 | #include <AliESDRun.h> |
0f67a57c | 21 | #include <AliStack.h> |
22 | #include <AliHeader.h> | |
23 | #include <AliGenEventHeader.h> | |
24 | #include <AliMultiplicity.h> | |
25 | #include <AliAnalysisManager.h> | |
26 | #include <AliMCEventHandler.h> | |
27 | #include <AliMCEvent.h> | |
28 | #include <AliESDInputHandler.h> | |
29 | ||
745d6088 | 30 | #include "AliESDtrackCuts.h" |
0f67a57c | 31 | #include "AliPWG0Helper.h" |
0f67a57c | 32 | #include "dNdEta/dNdEtaAnalysis.h" |
33 | #include "dNdEta/AlidNdEtaCorrection.h" | |
70fdd197 | 34 | #include "AliTriggerAnalysis.h" |
81be4ee8 | 35 | #include "AliPhysicsSelection.h" |
0f67a57c | 36 | |
37 | ClassImp(AlidNdEtaCorrectionTask) | |
38 | ||
1c15d51a | 39 | AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask() : |
40 | AliAnalysisTask(), | |
41 | fESD(0), | |
42 | fOutput(0), | |
43 | fOption(), | |
a7f69e56 | 44 | fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn)), |
70fdd197 | 45 | fTrigger(AliTriggerAnalysis::kMB1), |
69b09e3b | 46 | fFillPhi(kFALSE), |
47 | fDeltaPhiCut(-1), | |
1d107532 | 48 | fSymmetrize(kFALSE), |
81be4ee8 | 49 | fMultAxisEta1(kFALSE), |
50 | fDiffTreatment(AliPWG0Helper::kMCFlags), | |
1c15d51a | 51 | fSignMode(0), |
52 | fOnlyPrimaries(kFALSE), | |
69b09e3b | 53 | fStatError(0), |
81be4ee8 | 54 | fSystSkipParticles(kFALSE), |
1c15d51a | 55 | fEsdTrackCuts(0), |
56 | fdNdEtaCorrection(0), | |
57 | fdNdEtaAnalysisMC(0), | |
58 | fdNdEtaAnalysisESD(0), | |
59 | fPIDParticles(0), | |
60 | fPIDTracks(0), | |
61 | fVertexCorrelation(0), | |
51f6de65 | 62 | fVertexCorrelationShift(0), |
1c15d51a | 63 | fVertexProfile(0), |
64 | fVertexShift(0), | |
65 | fVertexShiftNorm(0), | |
66 | fEtaCorrelation(0), | |
51f6de65 | 67 | fEtaCorrelationShift(0), |
1c15d51a | 68 | fEtaProfile(0), |
69 | fEtaResolution(0), | |
69b09e3b | 70 | fDeltaPhiCorrelation(0), |
1c15d51a | 71 | fpTResolution(0), |
72 | fEsdTrackCutsPrim(0), | |
73 | fEsdTrackCutsSec(0), | |
74 | fTemp1(0), | |
75 | fTemp2(0), | |
76 | fMultAll(0), | |
77 | fMultTr(0), | |
78 | fMultVtx(0), | |
68fa248f | 79 | fEventStats(0), |
80 | fEsdTrackCutsCheck(0), | |
81 | fEtaCorrelationAllESD(0), | |
82 | fpTCorrelation(0), | |
83 | fpTCorrelationShift(0), | |
84 | fpTCorrelationAllESD(0), | |
85 | fpTCorrelationShiftAllESD(0), | |
86 | fPtMin(0.15), | |
87 | fPtMC(0), | |
88 | fEtaMC(0), | |
89 | fPtESD(0), | |
90 | fEtaESD(0), | |
91 | fVtxMC(0), | |
92 | fNumberEventMC(0), | |
93 | fNumberEvent(0), | |
94 | fEventNumber(-1), | |
95 | fWeightSecondaries(kFALSE) | |
1c15d51a | 96 | { |
97 | // | |
98 | // Constructor. Initialization of pointers | |
99 | // | |
100 | ||
69b09e3b | 101 | for (Int_t i=0; i<4; i++) |
102 | fdNdEtaCorrectionSpecial[i] = 0; | |
1c15d51a | 103 | |
104 | for (Int_t i=0; i<8; i++) | |
105 | fDeltaPhi[i] = 0; | |
81be4ee8 | 106 | |
107 | AliLog::SetClassDebugLevel("AlidNdEtaCorrectionTask", AliLog::kWarning); | |
1c15d51a | 108 | } |
109 | ||
0f67a57c | 110 | AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask(const char* opt) : |
111 | AliAnalysisTask("AlidNdEtaCorrectionTask", ""), | |
112 | fESD(0), | |
113 | fOutput(0), | |
114 | fOption(opt), | |
a7f69e56 | 115 | fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn)), |
70fdd197 | 116 | fTrigger(AliTriggerAnalysis::kMB1), |
69b09e3b | 117 | fFillPhi(kFALSE), |
118 | fDeltaPhiCut(0), | |
1d107532 | 119 | fSymmetrize(kFALSE), |
81be4ee8 | 120 | fMultAxisEta1(kFALSE), |
121 | fDiffTreatment(AliPWG0Helper::kMCFlags), | |
0f67a57c | 122 | fSignMode(0), |
3d7758c1 | 123 | fOnlyPrimaries(kFALSE), |
69b09e3b | 124 | fStatError(0), |
81be4ee8 | 125 | fSystSkipParticles(kFALSE), |
0f67a57c | 126 | fEsdTrackCuts(0), |
127 | fdNdEtaCorrection(0), | |
128 | fdNdEtaAnalysisMC(0), | |
129 | fdNdEtaAnalysisESD(0), | |
130 | fPIDParticles(0), | |
96fcc8a7 | 131 | fPIDTracks(0), |
54b096ef | 132 | fVertexCorrelation(0), |
51f6de65 | 133 | fVertexCorrelationShift(0), |
54b096ef | 134 | fVertexProfile(0), |
1c15d51a | 135 | fVertexShift(0), |
54b096ef | 136 | fVertexShiftNorm(0), |
3d7758c1 | 137 | fEtaCorrelation(0), |
51f6de65 | 138 | fEtaCorrelationShift(0), |
3d7758c1 | 139 | fEtaProfile(0), |
1c15d51a | 140 | fEtaResolution(0), |
69b09e3b | 141 | fDeltaPhiCorrelation(0), |
1c15d51a | 142 | fpTResolution(0), |
143 | fEsdTrackCutsPrim(0), | |
144 | fEsdTrackCutsSec(0), | |
145 | fTemp1(0), | |
146 | fTemp2(0), | |
4ef701a0 | 147 | fMultAll(0), |
148 | fMultTr(0), | |
50ec344d | 149 | fMultVtx(0), |
68fa248f | 150 | fEventStats(0), |
151 | fEsdTrackCutsCheck(0), | |
152 | fEtaCorrelationAllESD(0), | |
153 | fpTCorrelation(0), | |
154 | fpTCorrelationShift(0), | |
155 | fpTCorrelationAllESD(0), | |
156 | fpTCorrelationShiftAllESD(0), | |
157 | fPtMin(0.15), | |
158 | fPtMC(0), | |
159 | fEtaMC(0), | |
160 | fPtESD(0), | |
161 | fEtaESD(0), | |
162 | fVtxMC(0), | |
163 | fNumberEventMC(0), | |
164 | fNumberEvent(0), | |
165 | fEventNumber(-1), | |
166 | fWeightSecondaries(kFALSE) | |
0f67a57c | 167 | { |
168 | // | |
169 | // Constructor. Initialization of pointers | |
170 | // | |
171 | ||
172 | // Define input and output slots here | |
173 | DefineInput(0, TChain::Class()); | |
174 | DefineOutput(0, TList::Class()); | |
96fcc8a7 | 175 | |
69b09e3b | 176 | for (Int_t i=0; i<4; i++) |
177 | fdNdEtaCorrectionSpecial[i] = 0; | |
4ef701a0 | 178 | |
3d7758c1 | 179 | for (Int_t i=0; i<8; i++) |
4ef701a0 | 180 | fDeltaPhi[i] = 0; |
81be4ee8 | 181 | |
182 | AliLog::SetClassDebugLevel("AlidNdEtaCorrectionTask", AliLog::kWarning); | |
0f67a57c | 183 | } |
184 | ||
185 | AlidNdEtaCorrectionTask::~AlidNdEtaCorrectionTask() | |
186 | { | |
187 | // | |
188 | // Destructor | |
189 | // | |
190 | ||
191 | // histograms are in the output list and deleted when the output | |
192 | // list is deleted by the TSelector dtor | |
df974e92 | 193 | |
194 | if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) { | |
0f67a57c | 195 | delete fOutput; |
196 | fOutput = 0; | |
197 | } | |
df974e92 | 198 | |
0f67a57c | 199 | } |
200 | ||
201 | //________________________________________________________________________ | |
202 | void AlidNdEtaCorrectionTask::ConnectInputData(Option_t *) | |
203 | { | |
204 | // Connect ESD | |
205 | // Called once | |
206 | ||
207 | Printf("AlidNdEtaCorrectionTask::ConnectInputData called"); | |
208 | ||
c17301f3 | 209 | AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); |
210 | ||
211 | if (!esdH) { | |
212 | Printf("ERROR: Could not get ESDInputHandler"); | |
0f67a57c | 213 | } else { |
c17301f3 | 214 | fESD = esdH->GetEvent(); |
0f67a57c | 215 | |
c17301f3 | 216 | // Enable only the needed branches |
217 | esdH->SetActiveBranches("AliESDHeader Vertex"); | |
0f67a57c | 218 | |
a7f69e56 | 219 | if (fAnalysisMode & AliPWG0Helper::kSPD) |
c17301f3 | 220 | esdH->SetActiveBranches("AliESDHeader Vertex AliMultiplicity"); |
0f67a57c | 221 | |
a7f69e56 | 222 | if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS) { |
c17301f3 | 223 | esdH->SetActiveBranches("AliESDHeader Vertex Tracks"); |
0f67a57c | 224 | } |
0f67a57c | 225 | } |
96fcc8a7 | 226 | |
227 | // disable info messages of AliMCEvent (per event) | |
228 | AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1); | |
0f67a57c | 229 | } |
230 | ||
231 | void AlidNdEtaCorrectionTask::CreateOutputObjects() | |
232 | { | |
233 | // create result objects and add to output list | |
234 | ||
68fa248f | 235 | AliDebug(2,Form("*********************************** fOption = %s", fOption.Data())); |
0f67a57c | 236 | if (fOption.Contains("only-positive")) |
237 | { | |
238 | Printf("INFO: Processing only positive particles."); | |
239 | fSignMode = 1; | |
240 | } | |
241 | else if (fOption.Contains("only-negative")) | |
242 | { | |
243 | Printf("INFO: Processing only negative particles."); | |
244 | fSignMode = -1; | |
245 | } | |
69b09e3b | 246 | |
247 | if (fOption.Contains("stat-error-1")) | |
248 | { | |
249 | Printf("INFO: Evaluation statistical errors. Mode: 1."); | |
250 | fStatError = 1; | |
251 | } | |
252 | else if (fOption.Contains("stat-error-2")) | |
253 | { | |
254 | Printf("INFO: Evaluation statistical errors. Mode: 2."); | |
255 | fStatError = 2; | |
256 | } | |
0f67a57c | 257 | |
0f67a57c | 258 | fOutput = new TList; |
259 | fOutput->SetOwner(); | |
260 | ||
770a1f1d | 261 | fdNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction", fAnalysisMode); |
0f67a57c | 262 | fOutput->Add(fdNdEtaCorrection); |
263 | ||
1c15d51a | 264 | fPIDParticles = new TH1F("fPIDParticles", "PID of generated primary particles", 10001, -5000.5, 5000.5); |
0f67a57c | 265 | fOutput->Add(fPIDParticles); |
266 | ||
1c15d51a | 267 | fPIDTracks = new TH1F("fPIDTracks", "MC PID of reconstructed tracks", 10001, -5000.5, 5000.5); |
0f67a57c | 268 | fOutput->Add(fPIDTracks); |
269 | ||
770a1f1d | 270 | fdNdEtaAnalysisMC = new dNdEtaAnalysis("dndetaMC", "dndetaMC", fAnalysisMode); |
0f67a57c | 271 | fOutput->Add(fdNdEtaAnalysisMC); |
272 | ||
770a1f1d | 273 | fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndetaESD", "dndetaESD", fAnalysisMode); |
0f67a57c | 274 | fOutput->Add(fdNdEtaAnalysisESD); |
96fcc8a7 | 275 | |
1c15d51a | 276 | if (fEsdTrackCuts) |
277 | { | |
278 | fEsdTrackCutsPrim = dynamic_cast<AliESDtrackCuts*> (fEsdTrackCuts->Clone("fEsdTrackCutsPrim")); | |
279 | fEsdTrackCutsSec = dynamic_cast<AliESDtrackCuts*> (fEsdTrackCuts->Clone("fEsdTrackCutsSec")); | |
68fa248f | 280 | fEsdTrackCutsCheck = dynamic_cast<AliESDtrackCuts*> (fEsdTrackCuts->Clone("fEsdTrackCutsCheck")); |
281 | fEsdTrackCutsCheck->SetPtRange(0.15); | |
282 | fEsdTrackCutsCheck->SetEtaRange(-1.,1.); | |
1c15d51a | 283 | fOutput->Add(fEsdTrackCutsPrim); |
284 | fOutput->Add(fEsdTrackCutsSec); | |
285 | } | |
286 | ||
96fcc8a7 | 287 | if (fOption.Contains("process-types")) { |
69b09e3b | 288 | fdNdEtaCorrectionSpecial[0] = new AlidNdEtaCorrection("dndeta_correction_ND", "dndeta_correction_ND", fAnalysisMode); |
289 | fdNdEtaCorrectionSpecial[1] = new AlidNdEtaCorrection("dndeta_correction_SD", "dndeta_correction_SD", fAnalysisMode); | |
290 | fdNdEtaCorrectionSpecial[2] = new AlidNdEtaCorrection("dndeta_correction_DD", "dndeta_correction_DD", fAnalysisMode); | |
96fcc8a7 | 291 | |
69b09e3b | 292 | fOutput->Add(fdNdEtaCorrectionSpecial[0]); |
293 | fOutput->Add(fdNdEtaCorrectionSpecial[1]); | |
294 | fOutput->Add(fdNdEtaCorrectionSpecial[2]); | |
295 | } | |
296 | ||
297 | if (fOption.Contains("particle-species")) { | |
298 | fdNdEtaCorrectionSpecial[0] = new AlidNdEtaCorrection("dndeta_correction_pi", "dndeta_correction_pi", fAnalysisMode); | |
299 | fdNdEtaCorrectionSpecial[1] = new AlidNdEtaCorrection("dndeta_correction_K", "dndeta_correction_K", fAnalysisMode); | |
300 | fdNdEtaCorrectionSpecial[2] = new AlidNdEtaCorrection("dndeta_correction_p", "dndeta_correction_p", fAnalysisMode); | |
301 | fdNdEtaCorrectionSpecial[3] = new AlidNdEtaCorrection("dndeta_correction_other", "dndeta_correction_other", fAnalysisMode); | |
302 | ||
303 | for (Int_t i=0; i<4; i++) | |
304 | fOutput->Add(fdNdEtaCorrectionSpecial[i]); | |
96fcc8a7 | 305 | } |
306 | ||
69b09e3b | 307 | |
81be4ee8 | 308 | //fTemp1 = new TH2F("fTemp1", "fTemp1", 4, 0.5, 4.5, 101, -1.5, 99.5); // nsd study |
309 | fTemp1 = new TH2F("fTemp1", "fTemp1", 300, -15, 15, 80, -2.0, 2.0); | |
1c15d51a | 310 | fOutput->Add(fTemp1); |
81be4ee8 | 311 | |
312 | fTemp2 = new TH2F("fTemp2", "fTemp2", 300, -15, 15, 80, -2.0, 2.0); | |
1c15d51a | 313 | fOutput->Add(fTemp2); |
54b096ef | 314 | |
0fc41645 | 315 | fVertexCorrelation = new TH2F("fVertexCorrelation", "fVertexCorrelation;MC z-vtx;ESD z-vtx", 120, -30, 30, 120, -30, 30); |
316 | fOutput->Add(fVertexCorrelation); | |
a7f69e56 | 317 | fVertexCorrelationShift = new TH3F("fVertexCorrelationShift", "fVertexCorrelationShift;MC z-vtx;MC z-vtx - ESD z-vtx;rec. tracks", 120, -30, 30, 100, -1, 1, 100, -0.5, 99.5); |
51f6de65 | 318 | fOutput->Add(fVertexCorrelationShift); |
54b096ef | 319 | fVertexProfile = new TProfile("fVertexProfile", "fVertexProfile;MC z-vtx;MC z-vtx - ESD z-vtx", 40, -20, 20); |
0fc41645 | 320 | fOutput->Add(fVertexProfile); |
1c15d51a | 321 | fVertexShift = new TH1F("fVertexShift", "fVertexShift;(MC z-vtx - ESD z-vtx);Entries", 201, -2, 2); |
322 | fOutput->Add(fVertexShift); | |
81be4ee8 | 323 | fVertexShiftNorm = new TH2F("fVertexShiftNorm", "fVertexShiftNorm;(MC z-vtx - ESD z-vtx);rec. tracks;Entries", 200, -100, 100, 100, -0.5, 99.5); |
0fc41645 | 324 | fOutput->Add(fVertexShiftNorm); |
325 | ||
3d7758c1 | 326 | fEtaCorrelation = new TH2F("fEtaCorrelation", "fEtaCorrelation;MC #eta;ESD #eta", 120, -3, 3, 120, -3, 3); |
0fc41645 | 327 | fOutput->Add(fEtaCorrelation); |
68fa248f | 328 | fEtaCorrelationAllESD = new TH2F("fEtaCorrelationAllESD", "fEtaCorrelationAllESD;MC #eta;ESD #eta", 120, -3, 3, 120, -3, 3); |
329 | fOutput->Add(fEtaCorrelationAllESD); | |
69b09e3b | 330 | fEtaCorrelationShift = new TH2F("fEtaCorrelationShift", "fEtaCorrelationShift;MC #eta;MC #eta - ESD #eta", 120, -3, 3, 100, -0.1, 0.1); |
51f6de65 | 331 | fOutput->Add(fEtaCorrelationShift); |
3d7758c1 | 332 | fEtaProfile = new TProfile("fEtaProfile", "fEtaProfile;MC #eta;MC #eta - ESD #eta", 120, -3, 3); |
0fc41645 | 333 | fOutput->Add(fEtaProfile); |
1c15d51a | 334 | fEtaResolution = new TH1F("fEtaResolution", "fEtaResolution;MC #eta - ESD #eta", 201, -0.2, 0.2); |
335 | fOutput->Add(fEtaResolution); | |
336 | ||
69b09e3b | 337 | fpTResolution = new TH2F("fpTResolution", ";MC p_{T} (GeV/c);(MC p_{T} - ESD p_{T}) / MC p_{T}", 160, 0, 20, 201, -0.2, 0.2); |
1c15d51a | 338 | fOutput->Add(fpTResolution); |
3d7758c1 | 339 | |
df974e92 | 340 | fpTCorrelation = new TH2F("fpTCorrelation", "fpTCorrelation;MC p_{T} (GeV/c);ESD p_{T}", 160, 0, 20, 160, 0, 20); |
341 | fOutput->Add(fpTCorrelation); | |
342 | fpTCorrelationShift = new TH2F("fpTCorrelationShift", "fpTCorrelationShift;MC p_{T} (GeV/c);MC p_{T} - ESD p_{T}", 160, 0, 20, 100, -1, 1); | |
343 | fOutput->Add(fpTCorrelationShift); | |
68fa248f | 344 | fpTCorrelationAllESD = new TH2F("fpTCorrelationAllESD", "fpTCorrelationAllESD;MC p_{T} (GeV/c);ESD p_{T}", 160, 0, 20, 160, 0, 20); |
345 | fOutput->Add(fpTCorrelationAllESD); | |
346 | fpTCorrelationShiftAllESD = new TH2F("fpTCorrelationShiftAllESD", "fpTCorrelationShiftAllESD;MC p_{T} (GeV/c);MC p_{T} - ESD p_{T}", 160, 0, 20, 100, -1, 1); | |
347 | fOutput->Add(fpTCorrelationShiftAllESD); | |
348 | ||
4ef701a0 | 349 | fMultAll = new TH1F("fMultAll", "fMultAll", 500, -0.5, 499.5); |
0fc41645 | 350 | fOutput->Add(fMultAll); |
4ef701a0 | 351 | fMultVtx = new TH1F("fMultVtx", "fMultVtx", 500, -0.5, 499.5); |
0fc41645 | 352 | fOutput->Add(fMultVtx); |
4ef701a0 | 353 | fMultTr = new TH1F("fMultTr", "fMultTr", 500, -0.5, 499.5); |
0fc41645 | 354 | fOutput->Add(fMultTr); |
4ef701a0 | 355 | |
3d7758c1 | 356 | for (Int_t i=0; i<8; i++) |
0fc41645 | 357 | { |
69b09e3b | 358 | fDeltaPhi[i] = new TH2F(Form("fDeltaPhi_%d", i), ";#Delta phi;#phi;Entries", 2000, -0.1, 0.1, 18 * 5, 0, TMath::TwoPi()); |
0fc41645 | 359 | fOutput->Add(fDeltaPhi[i]); |
360 | } | |
50ec344d | 361 | |
81be4ee8 | 362 | fEventStats = new TH2F("fEventStats", "fEventStats;event type;status;count", 109, -6.5, 102.5, 4, -0.5, 3.5); |
0fc41645 | 363 | fOutput->Add(fEventStats); |
69b09e3b | 364 | fEventStats->GetXaxis()->SetBinLabel(1, "INEL"); // x = -5 |
365 | fEventStats->GetXaxis()->SetBinLabel(2, "NSD"); // x = -4 | |
366 | fEventStats->GetXaxis()->SetBinLabel(3, "ND"); // x = -3 | |
367 | fEventStats->GetXaxis()->SetBinLabel(4, "SD"); // x = -2 | |
368 | fEventStats->GetXaxis()->SetBinLabel(5, "DD"); // x = -1 | |
81be4ee8 | 369 | |
370 | fEventStats->GetXaxis()->SetBinLabel(108, "INEL=0"); // x = -101 | |
371 | fEventStats->GetXaxis()->SetBinLabel(109, "INEL>0"); // x = -102 | |
69b09e3b | 372 | |
81be4ee8 | 373 | for (Int_t i=-1; i<100; i++) |
69b09e3b | 374 | fEventStats->GetXaxis()->SetBinLabel(7+i, Form("%d", i)); |
50ec344d | 375 | |
376 | fEventStats->GetYaxis()->SetBinLabel(1, "nothing"); | |
377 | fEventStats->GetYaxis()->SetBinLabel(2, "trg"); | |
378 | fEventStats->GetYaxis()->SetBinLabel(3, "vtx"); | |
379 | fEventStats->GetYaxis()->SetBinLabel(4, "trgvtx"); | |
1c15d51a | 380 | |
381 | if (fEsdTrackCuts) | |
382 | { | |
383 | fEsdTrackCuts->SetName("fEsdTrackCuts"); | |
81be4ee8 | 384 | // TODO like this we send an empty object back... |
385 | fOutput->Add(fEsdTrackCuts->Clone()); | |
1c15d51a | 386 | } |
68fa248f | 387 | fPtMC = new TH1F("fPtMC", "Pt from MC for selected tracks;MC p_{T} (GeV/c)", 160, 0, 20); |
388 | fOutput->Add(fPtMC); | |
389 | fEtaMC = new TH1F("fEtaMC", "Eta from MC for selected tracks;MC #eta", 120, -3, 3); | |
390 | fOutput->Add(fEtaMC); | |
391 | fPtESD = new TH1F("fPtESD", "Pt from ESD for selected tracks;ESD p_{T} (GeV/c)", 160, 0, 20); | |
392 | fOutput->Add(fPtESD); | |
393 | fEtaESD = new TH1F("fEtaESD", "Eta from ESD for selected tracks;ESD #eta", 120, -3, 3); | |
394 | fOutput->Add(fEtaESD); | |
395 | ||
396 | fVtxMC = new TH1F("fVtxMC", "Vtx,z from MC for all events;MC vtx_z (cm)", 100, -30, 30); | |
397 | fOutput->Add(fVtxMC); | |
398 | ||
399 | fNumberEventMC = new TH1F("fNumberEventMC","Number of event accepted at MC level",600000,-0.5,600000-0.5); | |
400 | fOutput->Add(fNumberEventMC); | |
401 | ||
402 | fNumberEvent = new TH1F("fNumberEvent","Number of event accepted at Reco level",600000,-0.5,600000-0.5); | |
403 | fOutput->Add(fNumberEvent); | |
0f67a57c | 404 | } |
405 | ||
406 | void AlidNdEtaCorrectionTask::Exec(Option_t*) | |
407 | { | |
408 | // process the event | |
409 | ||
68fa248f | 410 | fEventNumber++; |
0f67a57c | 411 | // Check prerequisites |
412 | if (!fESD) | |
413 | { | |
414 | AliDebug(AliLog::kError, "ESD branch not available"); | |
415 | return; | |
416 | } | |
417 | ||
3d7758c1 | 418 | if (fOnlyPrimaries) |
419 | Printf("WARNING: Processing only primaries. For systematical studies only!"); | |
69b09e3b | 420 | |
421 | if (fStatError > 0) | |
422 | Printf("WARNING: Statistical error evaluation active!"); | |
a7f69e56 | 423 | |
81be4ee8 | 424 | AliInputEventHandler* inputHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); |
425 | if (!inputHandler) | |
426 | { | |
427 | Printf("ERROR: Could not receive input handler"); | |
428 | return; | |
429 | } | |
430 | ||
431 | Bool_t eventTriggered = inputHandler->IsEventSelected(); | |
54b096ef | 432 | |
81be4ee8 | 433 | static AliTriggerAnalysis* triggerAnalysis = 0; |
434 | if (!triggerAnalysis) | |
435 | { | |
436 | AliPhysicsSelection* physicsSelection = dynamic_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection()); | |
437 | if (physicsSelection) | |
438 | triggerAnalysis = physicsSelection->GetTriggerAnalysis(); | |
439 | } | |
440 | ||
441 | if (eventTriggered) | |
442 | eventTriggered = triggerAnalysis->IsTriggerFired(fESD, fTrigger); | |
443 | ||
3d7758c1 | 444 | if (!eventTriggered) |
445 | Printf("No trigger"); | |
446 | ||
54b096ef | 447 | // post the data already here |
448 | PostData(0, fOutput); | |
69b09e3b | 449 | |
54b096ef | 450 | // MC info |
451 | AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); | |
452 | if (!eventHandler) { | |
453 | Printf("ERROR: Could not retrieve MC event handler"); | |
454 | return; | |
455 | } | |
456 | ||
457 | AliMCEvent* mcEvent = eventHandler->MCEvent(); | |
458 | if (!mcEvent) { | |
459 | Printf("ERROR: Could not retrieve MC event"); | |
460 | return; | |
461 | } | |
462 | ||
463 | AliStack* stack = mcEvent->Stack(); | |
464 | if (!stack) | |
465 | { | |
466 | AliDebug(AliLog::kError, "Stack not available"); | |
467 | return; | |
468 | } | |
469 | ||
470 | AliHeader* header = mcEvent->Header(); | |
471 | if (!header) | |
472 | { | |
473 | AliDebug(AliLog::kError, "Header not available"); | |
474 | return; | |
475 | } | |
476 | ||
81be4ee8 | 477 | // get process type |
478 | Int_t processType = AliPWG0Helper::GetEventProcessType(fESD, header, stack, fDiffTreatment); | |
479 | ||
6b62a9c7 | 480 | AliDebug(AliLog::kDebug+1, Form("Found process type %d", processType)); |
54b096ef | 481 | |
6b62a9c7 | 482 | if (processType == AliPWG0Helper::kInvalidProcess) |
81be4ee8 | 483 | { |
484 | AliDebug(AliLog::kWarning, "Unknown process type. Setting to ND"); | |
485 | processType = AliPWG0Helper::kND; | |
486 | } | |
54b096ef | 487 | |
488 | // get the MC vertex | |
489 | AliGenEventHeader* genHeader = header->GenEventHeader(); | |
490 | TArrayF vtxMC(3); | |
491 | genHeader->PrimaryVertex(vtxMC); | |
68fa248f | 492 | fVtxMC->Fill(vtxMC[2]); |
493 | AliDebug(2,Form("MC vtx: x = %.6f, y = %.6f, z = %.6f",vtxMC[0],vtxMC[1],vtxMC[2])); | |
54b096ef | 494 | |
495 | // get the ESD vertex | |
770a1f1d | 496 | Bool_t eventVertex = kFALSE; |
a7f69e56 | 497 | Double_t vtx[3]; |
498 | const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode); | |
499 | if (vtxESD && AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode)) | |
54b096ef | 500 | { |
501 | vtxESD->GetXYZ(vtx); | |
81be4ee8 | 502 | eventVertex = kTRUE; |
503 | ||
504 | // remove vertices outside +- 15 cm | |
505 | if (TMath::Abs(vtx[2]) > 15) | |
506 | { | |
507 | eventVertex = kFALSE; | |
508 | vtxESD = 0; | |
509 | } | |
567160d6 | 510 | } |
a7f69e56 | 511 | else |
512 | vtxESD = 0; | |
513 | ||
0f67a57c | 514 | // create list of (label, eta, pt) tuples |
515 | Int_t inputCount = 0; | |
516 | Int_t* labelArr = 0; | |
3d7758c1 | 517 | Int_t* labelArr2 = 0; // only for case of SPD |
0f67a57c | 518 | Float_t* etaArr = 0; |
69b09e3b | 519 | Float_t* thirdDimArr = 0; |
4ef701a0 | 520 | Float_t* deltaPhiArr = 0; |
a7f69e56 | 521 | if (fAnalysisMode & AliPWG0Helper::kSPD) |
0f67a57c | 522 | { |
81be4ee8 | 523 | if (vtxESD) |
0f67a57c | 524 | { |
81be4ee8 | 525 | // get tracklets |
526 | const AliMultiplicity* mult = fESD->GetMultiplicity(); | |
527 | if (!mult) | |
528 | { | |
529 | AliDebug(AliLog::kError, "AliMultiplicity not available"); | |
530 | return; | |
531 | } | |
532 | ||
533 | labelArr = new Int_t[mult->GetNumberOfTracklets()]; | |
534 | labelArr2 = new Int_t[mult->GetNumberOfTracklets()]; | |
535 | etaArr = new Float_t[mult->GetNumberOfTracklets()]; | |
536 | thirdDimArr = new Float_t[mult->GetNumberOfTracklets()]; | |
537 | deltaPhiArr = new Float_t[mult->GetNumberOfTracklets()]; | |
538 | ||
539 | Bool_t foundInEta10 = kFALSE; | |
540 | ||
541 | // get multiplicity from SPD tracklets | |
542 | for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i) | |
543 | { | |
544 | //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i)); | |
545 | ||
546 | Float_t phi = mult->GetPhi(i); | |
547 | if (phi < 0) | |
548 | phi += TMath::Pi() * 2; | |
549 | Float_t deltaPhi = mult->GetDeltaPhi(i); | |
550 | ||
551 | if (TMath::Abs(deltaPhi) > 1) | |
552 | printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi); | |
553 | ||
554 | if (fOnlyPrimaries) | |
555 | if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0))) | |
556 | continue; | |
557 | ||
558 | if (fDeltaPhiCut > 0 && (TMath::Abs(deltaPhi) > fDeltaPhiCut || TMath::Abs(mult->GetDeltaTheta(i)) > fDeltaPhiCut / 0.08 * 0.025)) | |
3d7758c1 | 559 | continue; |
81be4ee8 | 560 | |
561 | if (fSystSkipParticles && gRandom->Uniform() < 0.0153) | |
562 | { | |
563 | Printf("Skipped tracklet!"); | |
564 | continue; | |
565 | } | |
566 | ||
567 | // TEST exclude potentially inefficient phi region | |
568 | //if (phi > 5.70 || phi < 0.06) | |
569 | // continue; | |
570 | ||
571 | // we have to repeat the trigger here, because the tracklet might have been kicked out fSystSkipParticles | |
572 | if (TMath::Abs(mult->GetEta(i)) < 1) | |
573 | foundInEta10 = kTRUE; | |
574 | ||
575 | etaArr[inputCount] = mult->GetEta(i); | |
576 | if (fSymmetrize) | |
577 | etaArr[inputCount] = TMath::Abs(etaArr[inputCount]); | |
578 | labelArr[inputCount] = mult->GetLabel(i, 0); | |
579 | labelArr2[inputCount] = mult->GetLabel(i, 1); | |
580 | thirdDimArr[inputCount] = phi; | |
581 | deltaPhiArr[inputCount] = deltaPhi; | |
582 | ++inputCount; | |
583 | } | |
69b09e3b | 584 | |
81be4ee8 | 585 | /* |
586 | for (Int_t i=0; i<mult->GetNumberOfSingleClusters(); ++i) | |
587 | { | |
588 | if (TMath::Abs(TMath::Log(TMath::Tan(mult->GetThetaSingle(i)/2.))) < 1); | |
589 | { | |
590 | foundInEta10 = kTRUE; | |
591 | break; | |
592 | } | |
593 | } | |
594 | */ | |
595 | ||
596 | if (fSystSkipParticles && (fTrigger & AliTriggerAnalysis::kOneParticle) && !foundInEta10) | |
597 | eventTriggered = kFALSE; | |
0f67a57c | 598 | } |
599 | } | |
a7f69e56 | 600 | else if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS) |
0f67a57c | 601 | { |
602 | if (!fEsdTrackCuts) | |
603 | { | |
604 | AliDebug(AliLog::kError, "fESDTrackCuts not available"); | |
605 | return; | |
606 | } | |
a7f69e56 | 607 | |
81be4ee8 | 608 | Bool_t foundInEta10 = kFALSE; |
609 | ||
69b09e3b | 610 | if (vtxESD) |
1c15d51a | 611 | { |
68fa248f | 612 | // control histograms on pT |
613 | Int_t nfromstack = stack->GetNtrack(); | |
614 | AliDebug(3,Form(" from stack we have %d tracks\n",nfromstack)); | |
615 | for (Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){ | |
616 | AliESDtrack* esdTrackcheck = dynamic_cast<AliESDtrack*> (fESD->GetTrack(itrack)); | |
617 | if (!esdTrackcheck){ | |
618 | AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", itrack)); | |
619 | continue; | |
620 | } | |
621 | if (fOnlyPrimaries){ | |
622 | Int_t label = TMath::Abs(esdTrackcheck->GetLabel()); | |
623 | AliDebug(4,Form("label = %d\n",label)); | |
624 | if (label == 0 || label > nfromstack) continue; | |
625 | if (stack->IsPhysicalPrimary(label) == kFALSE) continue; | |
626 | } | |
627 | ||
628 | Int_t label = TMath::Abs(esdTrackcheck->GetLabel()); | |
629 | if (label == 0 || label > nfromstack) continue; | |
630 | if (!stack->Particle(label)){ | |
631 | AliDebug(4,Form("WARNING: No particle for %d", label)); | |
632 | } | |
633 | else{ | |
634 | if (!fEsdTrackCuts->AcceptTrack(esdTrackcheck)){ | |
635 | TParticle* particle = stack->Particle(label); | |
636 | if (fEsdTrackCutsCheck->AcceptTrack(esdTrackcheck)){ | |
637 | //if (TMath::Abs(particle->Eta() < 0.8) && particle->Pt() > fPtMin){ | |
638 | Float_t ptMC = particle->Pt(); | |
639 | Float_t etaMC = particle->Eta(); | |
640 | Float_t ptESD = esdTrackcheck->Pt(); | |
641 | Float_t etaESD = esdTrackcheck->Eta(); | |
642 | fEtaCorrelationAllESD->Fill(etaMC,etaESD); | |
643 | fpTCorrelationAllESD->Fill(ptMC,ptESD); | |
644 | fpTCorrelationShiftAllESD->Fill(ptMC,ptMC-ptESD); | |
645 | } | |
646 | } | |
647 | } | |
648 | } // end loop over all ESDs | |
649 | ||
69b09e3b | 650 | // get multiplicity from ESD tracks |
a7f69e56 | 651 | TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, (fAnalysisMode & AliPWG0Helper::kTPC)); |
69b09e3b | 652 | Int_t nGoodTracks = list->GetEntries(); |
653 | ||
654 | Printf("Accepted %d tracks", nGoodTracks); | |
655 | ||
656 | labelArr = new Int_t[nGoodTracks]; | |
657 | labelArr2 = new Int_t[nGoodTracks]; | |
658 | etaArr = new Float_t[nGoodTracks]; | |
659 | thirdDimArr = new Float_t[nGoodTracks]; | |
660 | deltaPhiArr = new Float_t[nGoodTracks]; | |
661 | ||
662 | // loop over esd tracks | |
663 | for (Int_t i=0; i<nGoodTracks; i++) | |
1c15d51a | 664 | { |
69b09e3b | 665 | AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i)); |
666 | if (!esdTrack) | |
1c15d51a | 667 | { |
69b09e3b | 668 | AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i)); |
1c15d51a | 669 | continue; |
670 | } | |
69b09e3b | 671 | |
68fa248f | 672 | AliDebug(3,Form("particle %d: pt = %.6f, eta = %.6f",i,esdTrack->Pt(), esdTrack->Eta())); |
df974e92 | 673 | // 2 Options for INEL>0 trigger - choose one |
674 | // 1. HL | |
675 | //if (esdTrack->Pt() < 0.15) | |
676 | // foundInEta10 = kTRUE; | |
677 | // 2. MB Working Group definition | |
68fa248f | 678 | if (esdTrack->Pt() < fPtMin) |
81be4ee8 | 679 | continue; |
a7f69e56 | 680 | |
681 | if (fOnlyPrimaries) | |
682 | { | |
683 | Int_t label = TMath::Abs(esdTrack->GetLabel()); | |
684 | if (label == 0) | |
685 | continue; | |
686 | ||
687 | if (stack->IsPhysicalPrimary(label) == kFALSE) | |
688 | continue; | |
81be4ee8 | 689 | } |
690 | ||
68fa248f | 691 | Int_t label = TMath::Abs(esdTrack->GetLabel()); |
692 | if (!stack->Particle(label)){ | |
693 | AliDebug(3,Form("WARNING: No particle for %d", label)); | |
694 | } | |
695 | else{ | |
696 | TParticle* particle = stack->Particle(label); | |
697 | Float_t ptMC = particle->Pt(); | |
698 | Float_t etaMC = particle->Eta(); | |
699 | fPtMC->Fill(ptMC); | |
700 | fEtaMC->Fill(etaMC); | |
701 | fPtESD->Fill(esdTrack->Pt()); | |
702 | fEtaESD->Fill(esdTrack->Eta()); | |
703 | } | |
704 | ||
705 | // 2 Options for INEL>0 trigger - choose one | |
706 | // 1. HL | |
707 | //if (TMath::Abs(esdTrack->Eta()) < 1 && esdTrack->Pt() > 0.15) | |
708 | // foundInEta10 = kTRUE; | |
709 | // 2. MB Working Group definition | |
710 | if (TMath::Abs(esdTrack->Eta()) < 0.8 && esdTrack->Pt() > fPtMin) | |
81be4ee8 | 711 | foundInEta10 = kTRUE; |
68fa248f | 712 | |
69b09e3b | 713 | etaArr[inputCount] = esdTrack->Eta(); |
1d107532 | 714 | if (fSymmetrize) |
715 | etaArr[inputCount] = TMath::Abs(etaArr[inputCount]); | |
69b09e3b | 716 | labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel()); |
717 | labelArr2[inputCount] = labelArr[inputCount]; // no second label for tracks | |
718 | thirdDimArr[inputCount] = esdTrack->Pt(); | |
719 | deltaPhiArr[inputCount] = 0; // no delta phi for tracks | |
720 | ++inputCount; | |
721 | } | |
1c15d51a | 722 | |
69b09e3b | 723 | delete list; |
1c15d51a | 724 | |
81be4ee8 | 725 | // TODO this code crashes for TPCITS because particles are requested from the stack for some labels that are out of bound |
726 | if (0 && eventTriggered) | |
69b09e3b | 727 | { |
728 | // collect values for primaries and secondaries | |
729 | for (Int_t iTrack = 0; iTrack < fESD->GetNumberOfTracks(); iTrack++) | |
1c15d51a | 730 | { |
69b09e3b | 731 | AliESDtrack* track = 0; |
732 | ||
a7f69e56 | 733 | if (fAnalysisMode & AliPWG0Helper::kTPC) |
69b09e3b | 734 | track = AliESDtrackCuts::GetTPCOnlyTrack(fESD, iTrack); |
a7f69e56 | 735 | else if (fAnalysisMode & AliPWG0Helper::kTPCITS) |
69b09e3b | 736 | track = fESD->GetTrack(iTrack); |
737 | ||
738 | if (!track) | |
739 | continue; | |
740 | ||
741 | Int_t label = TMath::Abs(track->GetLabel()); | |
742 | if (!stack->Particle(label) || !stack->Particle(label)->GetPDG()) | |
743 | { | |
744 | Printf("WARNING: No particle for %d", label); | |
745 | if (stack->Particle(label)) | |
746 | stack->Particle(label)->Print(); | |
747 | continue; | |
748 | } | |
749 | ||
750 | if (stack->Particle(label)->GetPDG()->Charge() == 0) | |
751 | continue; | |
752 | ||
a7f69e56 | 753 | if (TMath::Abs(track->Eta()) < 0.8 && track->Pt() > 0.15) |
69b09e3b | 754 | { |
755 | if (stack->IsPhysicalPrimary(label)) | |
1c15d51a | 756 | { |
69b09e3b | 757 | // primary |
758 | if (fEsdTrackCutsPrim->AcceptTrack(track)) | |
759 | { | |
a7f69e56 | 760 | // if (AliESDtrackCuts::GetSigmaToVertex(track) > 900) |
761 | // { | |
762 | // Printf("Track %d has nsigma of %f. Printing track and vertex...", iTrack, AliESDtrackCuts::GetSigmaToVertex(track)); | |
763 | // Float_t b[2]; | |
764 | // Float_t r[3]; | |
765 | // track->GetImpactParameters(b, r); | |
766 | // Printf("Impact parameter %f %f and resolution: %f %f %f", b[0], b[1], r[0], r[1], r[2]); | |
767 | // track->Print(""); | |
768 | // if (vtxESD) | |
769 | // vtxESD->Print(); | |
770 | // } | |
69b09e3b | 771 | } |
1c15d51a | 772 | } |
69b09e3b | 773 | else |
774 | { | |
775 | // secondary | |
776 | fEsdTrackCutsSec->AcceptTrack(track); | |
1c15d51a | 777 | } |
778 | } | |
69b09e3b | 779 | |
780 | // TODO mem leak in the continue statements above | |
a7f69e56 | 781 | if (fAnalysisMode & AliPWG0Helper::kTPC) |
69b09e3b | 782 | delete track; |
1c15d51a | 783 | } |
1c15d51a | 784 | } |
785 | } | |
81be4ee8 | 786 | |
787 | if (!foundInEta10) | |
788 | eventTriggered = kFALSE; | |
68fa248f | 789 | else{ |
790 | //Printf("The event triggered. Its number in file is %d",fESD->GetEventNumberInFile()); | |
791 | fNumberEvent->Fill(fESD->GetEventNumberInFile()); | |
792 | } | |
0f67a57c | 793 | } |
54b096ef | 794 | else |
0f67a57c | 795 | return; |
0f67a57c | 796 | |
81be4ee8 | 797 | // fill process type |
798 | Int_t biny = (Int_t) eventTriggered + 2 * (Int_t) eventVertex; | |
799 | // INEL | |
800 | fEventStats->Fill(-6, biny); | |
801 | // NSD | |
802 | if (processType != AliPWG0Helper::kSD) | |
803 | fEventStats->Fill(-5, biny); | |
804 | // SD, ND, DD | |
805 | if (processType == AliPWG0Helper::kND) | |
806 | fEventStats->Fill(-4, biny); | |
807 | if (processType == AliPWG0Helper::kSD) | |
808 | fEventStats->Fill(-3, biny); | |
809 | if (processType == AliPWG0Helper::kDD) | |
810 | fEventStats->Fill(-2, biny); | |
811 | ||
0f67a57c | 812 | // loop over mc particles |
813 | Int_t nPrim = stack->GetNprimary(); | |
4ef701a0 | 814 | Int_t nAccepted = 0; |
0f67a57c | 815 | |
81be4ee8 | 816 | Bool_t oneParticleEvent = kFALSE; |
817 | for (Int_t iMc = 0; iMc < nPrim; ++iMc) | |
818 | { | |
819 | //Printf("Getting particle %d", iMc); | |
820 | TParticle* particle = stack->Particle(iMc); | |
821 | ||
822 | if (!particle) | |
823 | continue; | |
824 | ||
825 | if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE) | |
826 | continue; | |
827 | ||
df974e92 | 828 | // for INEL > 0, MB Working Group definition use the second option |
829 | // 1. standard | |
830 | //if (TMath::Abs(particle->Eta()) < 1.0) | |
831 | // 2. MB Working Group definition | |
832 | if (TMath::Abs(particle->Eta()) < 0.8 && particle->Pt() > fPtMin) | |
81be4ee8 | 833 | { |
834 | oneParticleEvent = kTRUE; | |
df974e92 | 835 | fNumberEventMC->Fill(fESD->GetEventNumberInFile()); |
81be4ee8 | 836 | break; |
837 | } | |
838 | } | |
839 | ||
840 | if (TMath::Abs(vtxMC[2]) < 5.5) | |
841 | { | |
842 | if (oneParticleEvent) | |
843 | fEventStats->Fill(102, biny); | |
844 | else | |
845 | fEventStats->Fill(101, biny); | |
846 | } | |
847 | ||
0f67a57c | 848 | for (Int_t iMc = 0; iMc < nPrim; ++iMc) |
849 | { | |
850 | //Printf("Getting particle %d", iMc); | |
851 | TParticle* particle = stack->Particle(iMc); | |
852 | ||
853 | if (!particle) | |
854 | continue; | |
855 | ||
856 | if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE) | |
1c15d51a | 857 | { |
858 | //if (TMath::Abs(particle->GetPdgCode()) > 3000 && TMath::Abs(particle->Eta()) < 1.0) | |
859 | // fPIDParticles->Fill(particle->GetPdgCode()); | |
0f67a57c | 860 | continue; |
1c15d51a | 861 | } |
0f67a57c | 862 | |
863 | if (SignOK(particle->GetPDG()) == kFALSE) | |
864 | continue; | |
81be4ee8 | 865 | |
68fa248f | 866 | // for INEL > 0, MB Working Group definition use the second option |
867 | // 1. standard | |
868 | //if (fPIDParticles && TMath::Abs(particle->Eta()) < 1.0) | |
869 | // 2. MB Working Group definition | |
870 | if (fPIDParticles && TMath::Abs(particle->Eta()) < 0.8 && particle->Pt() > fPtMin) | |
1c15d51a | 871 | fPIDParticles->Fill(particle->GetPdgCode()); |
872 | ||
0f67a57c | 873 | Float_t eta = particle->Eta(); |
1d107532 | 874 | if (fSymmetrize) |
875 | eta = TMath::Abs(eta); | |
69b09e3b | 876 | |
877 | Float_t thirdDim = -1; | |
a7f69e56 | 878 | if (fAnalysisMode & AliPWG0Helper::kSPD) |
69b09e3b | 879 | { |
880 | if (fFillPhi) | |
881 | { | |
882 | thirdDim = particle->Phi(); | |
883 | } | |
884 | else | |
885 | thirdDim = inputCount; | |
886 | } | |
887 | else | |
888 | thirdDim = particle->Pt(); | |
889 | ||
890 | // calculate y | |
891 | //Float_t y = 0.5 * TMath::Log((particle->Energy() + particle->Pz()) / (particle->Energy() - particle->Pz())); | |
892 | //fTemp1->Fill(eta); | |
893 | //fTemp2->Fill(y); | |
0f67a57c | 894 | |
81be4ee8 | 895 | Int_t processType2 = processType; |
896 | if (oneParticleEvent) | |
897 | processType2 |= AliPWG0Helper::kOnePart; | |
df974e92 | 898 | |
81be4ee8 | 899 | fdNdEtaCorrection->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2); |
0f67a57c | 900 | |
69b09e3b | 901 | if (fOption.Contains("process-types")) |
96fcc8a7 | 902 | { |
903 | // non diffractive | |
6b62a9c7 | 904 | if (processType==AliPWG0Helper::kND) |
81be4ee8 | 905 | fdNdEtaCorrectionSpecial[0]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2); |
96fcc8a7 | 906 | |
907 | // single diffractive | |
6b62a9c7 | 908 | if (processType==AliPWG0Helper::kSD) |
81be4ee8 | 909 | fdNdEtaCorrectionSpecial[1]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2); |
96fcc8a7 | 910 | |
911 | // double diffractive | |
6b62a9c7 | 912 | if (processType==AliPWG0Helper::kDD) |
81be4ee8 | 913 | fdNdEtaCorrectionSpecial[2]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2); |
69b09e3b | 914 | } |
915 | ||
916 | if (fOption.Contains("particle-species")) | |
917 | { | |
918 | Int_t id = -1; | |
919 | switch (TMath::Abs(particle->GetPdgCode())) | |
920 | { | |
921 | case 211: id = 0; break; | |
922 | case 321: id = 1; break; | |
923 | case 2212: id = 2; break; | |
924 | default: id = 3; break; | |
925 | } | |
81be4ee8 | 926 | fdNdEtaCorrectionSpecial[id]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType2); |
96fcc8a7 | 927 | } |
928 | ||
0f67a57c | 929 | if (eventTriggered) |
930 | if (eventVertex) | |
51f6de65 | 931 | fdNdEtaAnalysisMC->FillTrack(vtxMC[2], eta, thirdDim); |
4ef701a0 | 932 | |
0fc41645 | 933 | // TODO this value might be needed lower for the SPD study (only used for control histograms anyway) |
ea441adf | 934 | if (TMath::Abs(eta) < 1.5) // && pt > 0.2) |
4ef701a0 | 935 | nAccepted++; |
936 | } | |
937 | ||
81be4ee8 | 938 | if (AliPWG0Helper::GetLastProcessType() >= -1) |
939 | fEventStats->Fill(AliPWG0Helper::GetLastProcessType(), biny); | |
ea441adf | 940 | |
4ef701a0 | 941 | fMultAll->Fill(nAccepted); |
942 | if (eventTriggered) { | |
943 | fMultTr->Fill(nAccepted); | |
944 | if (eventVertex) | |
945 | fMultVtx->Fill(nAccepted); | |
0f67a57c | 946 | } |
947 | ||
3d7758c1 | 948 | Int_t processed = 0; |
69b09e3b | 949 | |
950 | Bool_t* primCount = 0; | |
951 | if (fStatError > 0) | |
952 | { | |
953 | primCount = new Bool_t[nPrim]; | |
954 | for (Int_t i=0; i<nPrim; ++i) | |
955 | primCount[i] = kFALSE; | |
956 | } | |
3d7758c1 | 957 | |
70fdd197 | 958 | Int_t nEta05 = 0; |
81be4ee8 | 959 | Int_t nEta10 = 0; |
0f67a57c | 960 | for (Int_t i=0; i<inputCount; ++i) |
961 | { | |
70fdd197 | 962 | if (TMath::Abs(etaArr[i]) < 0.5) |
963 | nEta05++; | |
81be4ee8 | 964 | if (TMath::Abs(etaArr[i]) < 1.0) |
965 | nEta10++; | |
966 | } | |
70fdd197 | 967 | |
81be4ee8 | 968 | for (Int_t i=0; i<inputCount; ++i) |
969 | { | |
0f67a57c | 970 | Int_t label = labelArr[i]; |
3d7758c1 | 971 | Int_t label2 = labelArr2[i]; |
0f67a57c | 972 | |
973 | if (label < 0) | |
974 | { | |
0fc41645 | 975 | Printf("WARNING: cannot find corresponding mc particle for track(let) %d with label %d.", i, label); |
0f67a57c | 976 | continue; |
977 | } | |
978 | ||
979 | TParticle* particle = stack->Particle(label); | |
980 | if (!particle) | |
981 | { | |
982 | AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label)); | |
983 | continue; | |
984 | } | |
1c15d51a | 985 | |
68fa248f | 986 | // for INEL > 0, MB Working Group definition use the second option |
987 | // 1. standard | |
988 | //if (TMath::Abs(particle->Eta()) < 1.0) | |
989 | // 2. INEL > 0 MB Working Group definition | |
990 | if (TMath::Abs(particle->Eta()) < 0.8 && particle->Pt()>fPtMin) | |
81be4ee8 | 991 | fPIDTracks->Fill(particle->GetPdgCode()); |
3d7758c1 | 992 | |
993 | // find particle that is filled in the correction map | |
994 | // this should be the particle that has been reconstructed | |
995 | // for TPC this is clear and always identified by the track's label | |
996 | // for SPD the labels might not be identical. In this case the particle closest to the beam line that is a primary is taken: | |
997 | // L1 L2 (P = primary, S = secondary) | |
998 | // P P' : --> P | |
999 | // P S : --> P | |
1000 | // S P : --> P | |
1001 | // S S' : --> S | |
1002 | ||
1003 | if (label != label2 && label2 >= 0) | |
1004 | { | |
1005 | if (stack->IsPhysicalPrimary(label) == kFALSE && stack->IsPhysicalPrimary(label2)) | |
1006 | { | |
1007 | particle = stack->Particle(label2); | |
1008 | if (!particle) | |
1009 | { | |
1010 | AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label2)); | |
1011 | continue; | |
1012 | } | |
1013 | } | |
1014 | } | |
0f67a57c | 1015 | |
1016 | if (eventTriggered && eventVertex) | |
1017 | { | |
3d7758c1 | 1018 | if (SignOK(particle->GetPDG()) == kFALSE) |
1019 | continue; | |
1020 | ||
1021 | processed++; | |
1022 | ||
1c15d51a | 1023 | // resolutions |
1d107532 | 1024 | if (fSymmetrize) |
1025 | fEtaResolution->Fill(TMath::Abs(particle->Eta()) - etaArr[i]); | |
1026 | else | |
1027 | fEtaResolution->Fill(particle->Eta() - etaArr[i]); | |
69b09e3b | 1028 | |
68fa248f | 1029 | if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS){ |
1030 | // for INEL > 0, MB Working Group definition use the second option | |
1031 | // 1. standard | |
1032 | //if (TMath::Abs(particle->Eta() < 0.9) && particle->Pt() > 0) | |
1033 | // 2. INEL > 0 MB WOrking Group definition | |
1034 | if (TMath::Abs(particle->Eta() < 0.8) && particle->Pt() > 0){ | |
df974e92 | 1035 | fpTResolution->Fill(particle->Pt(), (particle->Pt() - thirdDimArr[i]) / particle->Pt()); |
68fa248f | 1036 | fpTCorrelation->Fill(particle->Pt(),thirdDimArr[i]); |
1037 | fpTCorrelationShift->Fill(particle->Pt(),particle->Pt()-thirdDimArr[i]); | |
1038 | } | |
1039 | } | |
1c15d51a | 1040 | |
51f6de65 | 1041 | Float_t eta = -999; |
1042 | Float_t thirdDim = -1; | |
1043 | ||
3d7758c1 | 1044 | Bool_t firstIsPrim = stack->IsPhysicalPrimary(label); |
a7f69e56 | 1045 | // in case of same label the MC values are filled, otherwise (background) the reconstructed values |
1046 | if (label == label2) | |
3d7758c1 | 1047 | { |
51f6de65 | 1048 | eta = particle->Eta(); |
1d107532 | 1049 | if (fSymmetrize) |
1050 | eta = TMath::Abs(eta); | |
69b09e3b | 1051 | |
a7f69e56 | 1052 | if (fAnalysisMode & AliPWG0Helper::kSPD) |
69b09e3b | 1053 | { |
1054 | if (fFillPhi) | |
1055 | { | |
1056 | thirdDim = particle->Phi(); | |
1057 | } | |
1058 | else | |
1059 | thirdDim = inputCount; | |
1060 | } | |
1061 | else | |
1062 | thirdDim = particle->Pt(); | |
3d7758c1 | 1063 | } |
1064 | else | |
1065 | { | |
a7f69e56 | 1066 | if (fAnalysisMode & AliPWG0Helper::kSPD && !fFillPhi) |
69b09e3b | 1067 | { |
81be4ee8 | 1068 | thirdDim = (fMultAxisEta1) ? nEta10 : inputCount; |
69b09e3b | 1069 | } |
1070 | else | |
1071 | thirdDim = thirdDimArr[i]; | |
1072 | ||
51f6de65 | 1073 | eta = etaArr[i]; |
3d7758c1 | 1074 | } |
745d6088 | 1075 | |
69b09e3b | 1076 | Bool_t fillTrack = kTRUE; |
51f6de65 | 1077 | |
69b09e3b | 1078 | // statistical error evaluation active? |
1079 | if (fStatError > 0) | |
1080 | { | |
1081 | Bool_t statErrorDecision = kFALSE; | |
1082 | ||
1083 | // primary and not yet counted | |
1084 | if (label == label2 && firstIsPrim && !primCount[label]) | |
1085 | { | |
1086 | statErrorDecision = kTRUE; | |
1087 | primCount[label] = kTRUE; | |
1088 | } | |
1089 | ||
1090 | // in case of 1 we count only unique primaries, in case of 2 all the rest | |
1091 | if (fStatError == 1) | |
1092 | { | |
1093 | fillTrack = statErrorDecision; | |
1094 | } | |
1095 | else if (fStatError == 2) | |
1096 | fillTrack = !statErrorDecision; | |
1097 | } | |
1098 | ||
1099 | if (fillTrack) | |
81be4ee8 | 1100 | { |
68fa248f | 1101 | Double_t weight = 1.; |
1102 | if (fWeightSecondaries){ | |
1103 | if (!firstIsPrim){ | |
1104 | weight = GetSecondaryCorrection(thirdDim); | |
1105 | } | |
1106 | } | |
1107 | fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], eta, thirdDim, weight); | |
81be4ee8 | 1108 | fTemp2->Fill(vtxMC[2], eta); |
1109 | } | |
1110 | ||
51f6de65 | 1111 | // eta comparison for tracklets with the same label (others are background) |
1112 | if (label == label2) | |
1113 | { | |
1d107532 | 1114 | Float_t eta2 = particle->Eta(); |
1115 | if (fSymmetrize) | |
1116 | eta2 = TMath::Abs(eta2); | |
1117 | ||
1118 | fEtaProfile->Fill(eta2, eta2 - etaArr[i]); | |
68fa248f | 1119 | fEtaCorrelation->Fill(eta2, etaArr[i]); |
1d107532 | 1120 | fEtaCorrelationShift->Fill(eta2, eta2 - etaArr[i]); |
51f6de65 | 1121 | } |
3d7758c1 | 1122 | |
1d107532 | 1123 | if (fSymmetrize) |
1124 | fdNdEtaAnalysisESD->FillTrack(vtxMC[2], TMath::Abs(particle->Eta()), thirdDim); | |
1125 | else | |
1126 | fdNdEtaAnalysisESD->FillTrack(vtxMC[2], particle->Eta(), thirdDim); | |
3d7758c1 | 1127 | |
69b09e3b | 1128 | if (fOption.Contains("process-types")) |
96fcc8a7 | 1129 | { |
1130 | // non diffractive | |
6b62a9c7 | 1131 | if (processType == AliPWG0Helper::kND) |
69b09e3b | 1132 | fdNdEtaCorrectionSpecial[0]->FillTrackedParticle(vtxMC[2], eta, thirdDim); |
96fcc8a7 | 1133 | |
1134 | // single diffractive | |
1c15d51a | 1135 | if (processType == AliPWG0Helper::kSD) |
69b09e3b | 1136 | fdNdEtaCorrectionSpecial[1]->FillTrackedParticle(vtxMC[2], eta, thirdDim); |
96fcc8a7 | 1137 | |
1138 | // double diffractive | |
6b62a9c7 | 1139 | if (processType == AliPWG0Helper::kDD) |
69b09e3b | 1140 | fdNdEtaCorrectionSpecial[2]->FillTrackedParticle(vtxMC[2], eta, thirdDim); |
1141 | } | |
1142 | ||
1143 | if (fOption.Contains("particle-species")) | |
1144 | { | |
1145 | // find mother first | |
1146 | TParticle* mother = AliPWG0Helper::FindPrimaryMother(stack, label); | |
1147 | ||
1148 | Int_t id = -1; | |
1149 | switch (TMath::Abs(mother->GetPdgCode())) | |
1150 | { | |
1151 | case 211: id = 0; break; | |
1152 | case 321: id = 1; break; | |
1153 | case 2212: id = 2; break; | |
1154 | default: id = 3; break; | |
1155 | } | |
1156 | fdNdEtaCorrectionSpecial[id]->FillTrackedParticle(vtxMC[2], eta, thirdDim); | |
96fcc8a7 | 1157 | } |
4ef701a0 | 1158 | |
745d6088 | 1159 | // control histograms |
3d7758c1 | 1160 | Int_t hist = -1; |
1161 | if (label == label2) | |
4ef701a0 | 1162 | { |
69b09e3b | 1163 | if (firstIsPrim) |
3d7758c1 | 1164 | { |
69b09e3b | 1165 | hist = 0; |
3d7758c1 | 1166 | } |
1167 | else | |
1168 | hist = 1; | |
1169 | } | |
1170 | else if (label2 >= 0) | |
1171 | { | |
1172 | Bool_t secondIsPrim = stack->IsPhysicalPrimary(label2); | |
1173 | if (firstIsPrim && secondIsPrim) | |
1174 | { | |
1175 | hist = 2; | |
1176 | } | |
1177 | else if (firstIsPrim && !secondIsPrim) | |
1178 | { | |
1179 | hist = 3; | |
69b09e3b | 1180 | |
3d7758c1 | 1181 | // check if secondary is caused by the primary or it is a fake combination |
1182 | //Printf("PS case --> Label 1 is %d, label 2 is %d", label, label2); | |
1183 | Int_t mother = label2; | |
1184 | while (stack->Particle(mother) && stack->Particle(mother)->GetMother(0) >= 0) | |
1185 | { | |
1186 | //Printf(" %d created by %d, %d", mother, stack->Particle(mother)->GetMother(0), stack->Particle(mother)->GetMother(1)); | |
1187 | if (stack->Particle(mother)->GetMother(0) == label) | |
1188 | { | |
1189 | /*Printf("The untraceable decay was:"); | |
1190 | Printf(" from:"); | |
1191 | particle->Print(); | |
1192 | Printf(" to (status code %d):", stack->Particle(mother)->GetStatusCode()); | |
1193 | stack->Particle(mother)->Print();*/ | |
1194 | hist = 4; | |
1195 | } | |
1196 | mother = stack->Particle(mother)->GetMother(0); | |
1197 | } | |
1198 | } | |
1199 | else if (!firstIsPrim && secondIsPrim) | |
1200 | { | |
1201 | hist = 5; | |
1202 | } | |
1203 | else if (!firstIsPrim && !secondIsPrim) | |
1204 | { | |
1205 | hist = 6; | |
1206 | } | |
1207 | ||
4ef701a0 | 1208 | } |
1209 | else | |
3d7758c1 | 1210 | hist = 7; |
1211 | ||
69b09e3b | 1212 | fDeltaPhi[hist]->Fill(deltaPhiArr[i], thirdDimArr[i]); |
0f67a57c | 1213 | } |
1214 | } | |
69b09e3b | 1215 | |
1216 | if (primCount) | |
1217 | { | |
1218 | delete[] primCount; | |
1219 | primCount = 0; | |
1220 | } | |
0f67a57c | 1221 | |
3d7758c1 | 1222 | if (processed < inputCount) |
1223 | Printf("Only %d out of %d track(let)s were used", processed, inputCount); | |
1224 | ||
a7f69e56 | 1225 | if (eventTriggered && eventVertex) |
1226 | { | |
1227 | Double_t diff = vtxMC[2] - vtx[2]; | |
1228 | fVertexShift->Fill(diff); | |
81be4ee8 | 1229 | |
a7f69e56 | 1230 | fVertexCorrelation->Fill(vtxMC[2], vtx[2]); |
1231 | fVertexCorrelationShift->Fill(vtxMC[2], vtxMC[2] - vtx[2], inputCount); | |
1232 | fVertexProfile->Fill(vtxMC[2], vtxMC[2] - vtx[2]); | |
81be4ee8 | 1233 | |
1234 | if (vtxESD->IsFromVertexerZ() && inputCount > 0) | |
1235 | fVertexShiftNorm->Fill(diff, vtxESD->GetNContributors()); | |
a7f69e56 | 1236 | } |
1237 | ||
81be4ee8 | 1238 | Int_t multAxis = inputCount; |
1239 | if (fMultAxisEta1) | |
1240 | multAxis = nEta10; | |
1241 | ||
0f67a57c | 1242 | if (eventTriggered && eventVertex) |
96fcc8a7 | 1243 | { |
81be4ee8 | 1244 | fdNdEtaAnalysisMC->FillEvent(vtxMC[2], multAxis); |
1245 | fdNdEtaAnalysisESD->FillEvent(vtxMC[2], multAxis); | |
96fcc8a7 | 1246 | } |
0f67a57c | 1247 | |
81be4ee8 | 1248 | Int_t processType2 = processType; |
1249 | if (oneParticleEvent) | |
1250 | processType2 |= AliPWG0Helper::kOnePart; | |
1251 | ||
1252 | // stuff regarding the vertex reco correction and trigger bias correction | |
1253 | fdNdEtaCorrection->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2); | |
96fcc8a7 | 1254 | |
69b09e3b | 1255 | if (fOption.Contains("process-types")) |
96fcc8a7 | 1256 | { |
1257 | // non diffractive | |
a7f69e56 | 1258 | if (processType == AliPWG0Helper::kND) |
81be4ee8 | 1259 | fdNdEtaCorrectionSpecial[0]->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2); |
96fcc8a7 | 1260 | |
1261 | // single diffractive | |
6b62a9c7 | 1262 | if (processType == AliPWG0Helper::kSD) |
81be4ee8 | 1263 | fdNdEtaCorrectionSpecial[1]->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2); |
96fcc8a7 | 1264 | |
1265 | // double diffractive | |
6b62a9c7 | 1266 | if (processType == AliPWG0Helper::kDD) |
81be4ee8 | 1267 | fdNdEtaCorrectionSpecial[2]->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2); |
0f67a57c | 1268 | } |
69b09e3b | 1269 | |
1270 | if (fOption.Contains("particle-species")) | |
1271 | for (Int_t id=0; id<4; id++) | |
81be4ee8 | 1272 | fdNdEtaCorrectionSpecial[id]->FillEvent(vtxMC[2], multAxis, eventTriggered, eventVertex, processType2); |
69b09e3b | 1273 | |
1274 | if (etaArr) | |
1275 | delete[] etaArr; | |
1276 | if (labelArr) | |
1277 | delete[] labelArr; | |
1278 | if (labelArr2) | |
1279 | delete[] labelArr2; | |
1280 | if (thirdDimArr) | |
1281 | delete[] thirdDimArr; | |
1282 | if (deltaPhiArr) | |
1283 | delete[] deltaPhiArr; | |
0f67a57c | 1284 | } |
1285 | ||
1286 | void AlidNdEtaCorrectionTask::Terminate(Option_t *) | |
1287 | { | |
1288 | // The Terminate() function is the last function to be called during | |
1289 | // a query. It always runs on the client, it can be used to present | |
1290 | // the results graphically or save the results to file. | |
1291 | ||
1292 | fOutput = dynamic_cast<TList*> (GetOutputData(0)); | |
1293 | if (!fOutput) { | |
1294 | Printf("ERROR: fOutput not available"); | |
1295 | return; | |
1296 | } | |
1297 | ||
1298 | fdNdEtaCorrection = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction")); | |
1299 | fdNdEtaAnalysisMC = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaMC")); | |
1300 | fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaESD")); | |
1301 | if (!fdNdEtaCorrection || !fdNdEtaAnalysisMC || !fdNdEtaAnalysisESD) | |
1302 | { | |
1303 | AliDebug(AliLog::kError, "Could not read object from output list"); | |
1304 | return; | |
1305 | } | |
1306 | ||
1307 | fdNdEtaCorrection->Finish(); | |
1308 | ||
1309 | TString fileName; | |
1310 | fileName.Form("correction_map%s.root", fOption.Data()); | |
1311 | TFile* fout = new TFile(fileName, "RECREATE"); | |
1312 | ||
1c15d51a | 1313 | fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCuts")); |
0f67a57c | 1314 | if (fEsdTrackCuts) |
1315 | fEsdTrackCuts->SaveHistograms("esd_track_cuts"); | |
1c15d51a | 1316 | |
1317 | fEsdTrackCutsPrim = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCutsPrim")); | |
1318 | if (fEsdTrackCutsPrim) | |
1319 | fEsdTrackCutsPrim->SaveHistograms("esd_track_cuts_primaries"); | |
1320 | ||
1321 | fEsdTrackCutsSec = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCutsSec")); | |
1322 | if (fEsdTrackCutsSec) | |
1323 | fEsdTrackCutsSec->SaveHistograms("esd_track_cuts_secondaries"); | |
1324 | ||
0f67a57c | 1325 | fdNdEtaCorrection->SaveHistograms(); |
1326 | fdNdEtaAnalysisMC->SaveHistograms(); | |
1327 | fdNdEtaAnalysisESD->SaveHistograms(); | |
1328 | ||
69b09e3b | 1329 | if (fOutput->FindObject("dndeta_correction_ND")) |
1330 | { | |
1331 | fdNdEtaCorrectionSpecial[0] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_ND")); | |
1332 | fdNdEtaCorrectionSpecial[1] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_SD")); | |
1333 | fdNdEtaCorrectionSpecial[2] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_DD")); | |
1334 | } | |
1335 | else | |
1336 | { | |
1337 | fdNdEtaCorrectionSpecial[0] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_pi")); | |
1338 | fdNdEtaCorrectionSpecial[1] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_K")); | |
1339 | fdNdEtaCorrectionSpecial[2] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_p")); | |
1340 | fdNdEtaCorrectionSpecial[3] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_other")); | |
1341 | } | |
1342 | for (Int_t i=0; i<4; ++i) | |
1343 | if (fdNdEtaCorrectionSpecial[i]) | |
1344 | fdNdEtaCorrectionSpecial[i]->SaveHistograms(); | |
96fcc8a7 | 1345 | |
1c15d51a | 1346 | fTemp1 = dynamic_cast<TH1*> (fOutput->FindObject("fTemp1")); |
1347 | if (fTemp1) | |
1348 | fTemp1->Write(); | |
1349 | fTemp2 = dynamic_cast<TH1*> (fOutput->FindObject("fTemp2")); | |
1350 | if (fTemp2) | |
1351 | fTemp2->Write(); | |
96fcc8a7 | 1352 | |
0fc41645 | 1353 | fVertexCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fVertexCorrelation")); |
54b096ef | 1354 | if (fVertexCorrelation) |
1355 | fVertexCorrelation->Write(); | |
a7f69e56 | 1356 | fVertexCorrelationShift = dynamic_cast<TH3F*> (fOutput->FindObject("fVertexCorrelationShift")); |
51f6de65 | 1357 | if (fVertexCorrelationShift) |
a7f69e56 | 1358 | { |
1359 | ((TH2*) fVertexCorrelationShift->Project3D("yx"))->FitSlicesY(); | |
51f6de65 | 1360 | fVertexCorrelationShift->Write(); |
a7f69e56 | 1361 | } |
0fc41645 | 1362 | fVertexProfile = dynamic_cast<TProfile*> (fOutput->FindObject("fVertexProfile")); |
54b096ef | 1363 | if (fVertexProfile) |
1364 | fVertexProfile->Write(); | |
1c15d51a | 1365 | fVertexShift = dynamic_cast<TH1F*> (fOutput->FindObject("fVertexShift")); |
1366 | if (fVertexShift) | |
1367 | fVertexShift->Write(); | |
a7f69e56 | 1368 | fVertexShiftNorm = dynamic_cast<TH2F*> (fOutput->FindObject("fVertexShiftNorm")); |
54b096ef | 1369 | if (fVertexShiftNorm) |
a7f69e56 | 1370 | { |
1371 | fVertexShiftNorm->ProjectionX(); | |
54b096ef | 1372 | fVertexShiftNorm->Write(); |
a7f69e56 | 1373 | } |
3d7758c1 | 1374 | |
0fc41645 | 1375 | fEtaCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelation")); |
3d7758c1 | 1376 | if (fEtaCorrelation) |
1377 | fEtaCorrelation->Write(); | |
68fa248f | 1378 | fEtaCorrelationAllESD = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelationAllESD")); |
1379 | if (fEtaCorrelationAllESD) | |
1380 | fEtaCorrelationAllESD->Write(); | |
51f6de65 | 1381 | fEtaCorrelationShift = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelationShift")); |
1382 | if (fEtaCorrelationShift) | |
a7f69e56 | 1383 | { |
1384 | fEtaCorrelationShift->FitSlicesY(); | |
51f6de65 | 1385 | fEtaCorrelationShift->Write(); |
a7f69e56 | 1386 | } |
0fc41645 | 1387 | fEtaProfile = dynamic_cast<TProfile*> (fOutput->FindObject("fEtaProfile")); |
3d7758c1 | 1388 | if (fEtaProfile) |
1389 | fEtaProfile->Write(); | |
1c15d51a | 1390 | fEtaResolution = dynamic_cast<TH1F*> (fOutput->FindObject("fEtaResolution")); |
1391 | if (fEtaResolution) | |
1392 | fEtaResolution->Write(); | |
68fa248f | 1393 | fpTCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fpTCorrelation")); |
1394 | if (fpTCorrelation) | |
1395 | fpTCorrelation->Write(); | |
1396 | fpTCorrelationShift = dynamic_cast<TH2F*> (fOutput->FindObject("fpTCorrelationShift")); | |
1397 | if (fpTCorrelationShift) | |
1398 | { | |
1399 | fpTCorrelationShift->FitSlicesY(); | |
1400 | fpTCorrelationShift->Write(); | |
1401 | } | |
69b09e3b | 1402 | fpTResolution = dynamic_cast<TH2F*> (fOutput->FindObject("fpTResolution")); |
1c15d51a | 1403 | if (fpTResolution) |
69b09e3b | 1404 | { |
1405 | fpTResolution->FitSlicesY(); | |
1c15d51a | 1406 | fpTResolution->Write(); |
69b09e3b | 1407 | } |
3d7758c1 | 1408 | |
68fa248f | 1409 | fpTCorrelationAllESD = dynamic_cast<TH2F*> (fOutput->FindObject("fpTCorrelationAllESD")); |
1410 | if (fpTCorrelationAllESD) | |
1411 | fpTCorrelationAllESD->Write(); | |
1412 | fpTCorrelationShiftAllESD = dynamic_cast<TH2F*> (fOutput->FindObject("fpTCorrelationShiftAllESD")); | |
1413 | if (fpTCorrelationShiftAllESD) | |
1414 | { | |
1415 | fpTCorrelationShiftAllESD->FitSlicesY(); | |
1416 | fpTCorrelationShiftAllESD->Write(); | |
1417 | } | |
0fc41645 | 1418 | fMultAll = dynamic_cast<TH1F*> (fOutput->FindObject("fMultAll")); |
4ef701a0 | 1419 | if (fMultAll) |
1420 | fMultAll->Write(); | |
1421 | ||
0fc41645 | 1422 | fMultTr = dynamic_cast<TH1F*> (fOutput->FindObject("fMultTr")); |
4ef701a0 | 1423 | if (fMultTr) |
1424 | fMultTr->Write(); | |
0fc41645 | 1425 | |
1426 | fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx")); | |
4ef701a0 | 1427 | if (fMultVtx) |
1428 | fMultVtx->Write(); | |
1429 | ||
3d7758c1 | 1430 | for (Int_t i=0; i<8; ++i) |
0fc41645 | 1431 | { |
69b09e3b | 1432 | fDeltaPhi[i] = dynamic_cast<TH2*> (fOutput->FindObject(Form("fDeltaPhi_%d", i))); |
4ef701a0 | 1433 | if (fDeltaPhi[i]) |
1434 | fDeltaPhi[i]->Write(); | |
0fc41645 | 1435 | } |
54b096ef | 1436 | |
0fc41645 | 1437 | fEventStats = dynamic_cast<TH2F*> (fOutput->FindObject("fEventStats")); |
50ec344d | 1438 | if (fEventStats) |
1439 | fEventStats->Write(); | |
1440 | ||
1c15d51a | 1441 | fPIDParticles = dynamic_cast<TH1F*> (fOutput->FindObject("fPIDParticles")); |
1442 | if (fPIDParticles) | |
1443 | fPIDParticles->Write(); | |
1444 | ||
1445 | fPIDTracks = dynamic_cast<TH1F*> (fOutput->FindObject("fPIDTracks")); | |
1446 | if (fPIDTracks) | |
1447 | fPIDTracks->Write(); | |
0f67a57c | 1448 | |
68fa248f | 1449 | fPtMC = dynamic_cast<TH1F*> (fOutput->FindObject("fPtMC")); |
1450 | if (fPtMC) | |
1451 | fPtMC->Write(); | |
1452 | fEtaMC = dynamic_cast<TH1F*> (fOutput->FindObject("fEtaMC")); | |
1453 | if (fEtaMC) | |
1454 | fEtaMC->Write(); | |
1455 | fPtESD = dynamic_cast<TH1F*> (fOutput->FindObject("fPtESD")); | |
1456 | if (fPtESD) | |
1457 | fPtESD->Write(); | |
1458 | fEtaESD = dynamic_cast<TH1F*> (fOutput->FindObject("fEtaESD")); | |
1459 | if (fEtaESD) | |
1460 | fEtaESD->Write(); | |
1461 | fVtxMC = dynamic_cast<TH1F*> (fOutput->FindObject("fVtxMC")); | |
1462 | if (fVtxMC) | |
1463 | fVtxMC->Write(); | |
1464 | ||
1465 | fNumberEventMC = dynamic_cast<TH1F*> (fOutput->FindObject("fNumberEventMC")); | |
1466 | if (fNumberEventMC) | |
1467 | fNumberEventMC->Write(); | |
1468 | fNumberEvent = dynamic_cast<TH1F*> (fOutput->FindObject("fNumberEvent")); | |
1469 | if (fNumberEvent) | |
1470 | fNumberEvent->Write(); | |
1471 | ||
1c15d51a | 1472 | //fdNdEtaCorrection->DrawHistograms(); |
0f67a57c | 1473 | |
68fa248f | 1474 | Printf("Writing result to %s", fileName.Data()); |
96fcc8a7 | 1475 | |
0f67a57c | 1476 | if (fPIDParticles && fPIDTracks) |
1477 | { | |
0f67a57c | 1478 | TDatabasePDG* pdgDB = new TDatabasePDG; |
1479 | ||
1480 | for (Int_t i=0; i <= fPIDParticles->GetNbinsX()+1; ++i) | |
1481 | if (fPIDParticles->GetBinContent(i) > 0) | |
1c15d51a | 1482 | { |
1483 | TObject* pdgParticle = pdgDB->GetParticle((Int_t) fPIDParticles->GetBinCenter(i)); | |
1484 | printf("PDG = %d (%s): generated: %d, reconstructed: %d, ratio: %f\n", (Int_t) fPIDParticles->GetBinCenter(i), (pdgParticle) ? pdgParticle->GetName() : "not found", (Int_t) fPIDParticles->GetBinContent(i), (Int_t) fPIDTracks->GetBinContent(i), ((fPIDTracks->GetBinContent(i) > 0) ? fPIDParticles->GetBinContent(i) / fPIDTracks->GetBinContent(i) : -1)); | |
1485 | } | |
0f67a57c | 1486 | |
1487 | delete pdgDB; | |
1488 | pdgDB = 0; | |
1489 | } | |
1c15d51a | 1490 | |
1491 | fout->Write(); | |
1492 | fout->Close(); | |
0f67a57c | 1493 | } |
1494 | ||
1495 | Bool_t AlidNdEtaCorrectionTask::SignOK(TParticlePDG* particle) | |
1496 | { | |
1497 | // returns if a particle with this sign should be counted | |
1498 | // this is determined by the value of fSignMode, which should have the same sign | |
1499 | // as the charge | |
1500 | // if fSignMode is 0 all particles are counted | |
1501 | ||
1502 | if (fSignMode == 0) | |
1503 | return kTRUE; | |
1504 | ||
1505 | if (!particle) | |
1506 | { | |
1507 | printf("WARNING: not counting a particle that does not have a pdg particle\n"); | |
1508 | return kFALSE; | |
1509 | } | |
1510 | ||
1511 | Double_t charge = particle->Charge(); | |
1512 | ||
1513 | if (fSignMode > 0) | |
1514 | if (charge < 0) | |
1515 | return kFALSE; | |
1516 | ||
1517 | if (fSignMode < 0) | |
1518 | if (charge > 0) | |
1519 | return kFALSE; | |
1520 | ||
1521 | return kTRUE; | |
1522 | } | |
96fcc8a7 | 1523 | |
68fa248f | 1524 | Double_t AlidNdEtaCorrectionTask::GetSecondaryCorrection(Double_t pt){ |
1525 | ||
1526 | // getting the data driven correction factor to correct for | |
1527 | // the underestimate of secondaries in Pythia | |
1528 | // (A. Dainese + J. Otwinowski | |
1529 | ||
1530 | if (pt <= 0.17) return 1.0; | |
1531 | if (pt <= 0.4) return GetLinearInterpolationValue(0.17,1.0,0.4,1.07, pt); | |
1532 | if (pt <= 0.6) return GetLinearInterpolationValue(0.4,1.07,0.6,1.25, pt); | |
1533 | if (pt <= 1.2) return GetLinearInterpolationValue(0.6,1.25,1.2,1.5, pt); | |
1534 | return 1.5; | |
1535 | ||
1536 | } | |
1537 | ||
1538 | Double_t AlidNdEtaCorrectionTask::GetLinearInterpolationValue(Double_t const x1, Double_t const y1, Double_t const x2, Double_t const y2, const Double_t pt) | |
1539 | { | |
1540 | ||
1541 | // | |
1542 | // linear interpolation to be used to get the secondary correction (see AlidNdEtaCorrectionTask::GetSecondaryCorrection) | |
1543 | // | |
1544 | ||
1545 | return ((y2-y1)/(x2-x1))*pt+(y2-(((y2-y1)/(x2-x1))*x2)); | |
1546 | } | |
1547 |