adding flag for non-field data
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AlidNdEtaCorrectionTask.cxx
CommitLineData
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>
0f67a57c 15
16#include <AliLog.h>
17#include <AliESDVertex.h>
18#include <AliESDEvent.h>
19#include <AliStack.h>
20#include <AliHeader.h>
21#include <AliGenEventHeader.h>
22#include <AliMultiplicity.h>
23#include <AliAnalysisManager.h>
24#include <AliMCEventHandler.h>
25#include <AliMCEvent.h>
26#include <AliESDInputHandler.h>
27
745d6088 28#include "AliESDtrackCuts.h"
0f67a57c 29#include "AliPWG0Helper.h"
0f67a57c 30#include "dNdEta/dNdEtaAnalysis.h"
31#include "dNdEta/AlidNdEtaCorrection.h"
32
33ClassImp(AlidNdEtaCorrectionTask)
34
1c15d51a 35AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask() :
36 AliAnalysisTask(),
37 fESD(0),
38 fOutput(0),
39 fOption(),
a7f69e56 40 fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn)),
1c15d51a 41 fTrigger(AliPWG0Helper::kMB1),
69b09e3b 42 fFillPhi(kFALSE),
43 fDeltaPhiCut(-1),
1c15d51a 44 fSignMode(0),
45 fOnlyPrimaries(kFALSE),
69b09e3b 46 fStatError(0),
1c15d51a 47 fEsdTrackCuts(0),
48 fdNdEtaCorrection(0),
49 fdNdEtaAnalysisMC(0),
50 fdNdEtaAnalysisESD(0),
51 fPIDParticles(0),
52 fPIDTracks(0),
53 fVertexCorrelation(0),
51f6de65 54 fVertexCorrelationShift(0),
1c15d51a 55 fVertexProfile(0),
56 fVertexShift(0),
57 fVertexShiftNorm(0),
58 fEtaCorrelation(0),
51f6de65 59 fEtaCorrelationShift(0),
1c15d51a 60 fEtaProfile(0),
61 fEtaResolution(0),
69b09e3b 62 fDeltaPhiCorrelation(0),
1c15d51a 63 fpTResolution(0),
64 fEsdTrackCutsPrim(0),
65 fEsdTrackCutsSec(0),
66 fTemp1(0),
67 fTemp2(0),
68 fMultAll(0),
69 fMultTr(0),
70 fMultVtx(0),
71 fEventStats(0)
72{
73 //
74 // Constructor. Initialization of pointers
75 //
76
69b09e3b 77 for (Int_t i=0; i<4; i++)
78 fdNdEtaCorrectionSpecial[i] = 0;
1c15d51a 79
80 for (Int_t i=0; i<8; i++)
81 fDeltaPhi[i] = 0;
82}
83
0f67a57c 84AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask(const char* opt) :
85 AliAnalysisTask("AlidNdEtaCorrectionTask", ""),
86 fESD(0),
87 fOutput(0),
88 fOption(opt),
a7f69e56 89 fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn)),
0fc41645 90 fTrigger(AliPWG0Helper::kMB1),
69b09e3b 91 fFillPhi(kFALSE),
92 fDeltaPhiCut(0),
0f67a57c 93 fSignMode(0),
3d7758c1 94 fOnlyPrimaries(kFALSE),
69b09e3b 95 fStatError(0),
0f67a57c 96 fEsdTrackCuts(0),
97 fdNdEtaCorrection(0),
98 fdNdEtaAnalysisMC(0),
99 fdNdEtaAnalysisESD(0),
100 fPIDParticles(0),
96fcc8a7 101 fPIDTracks(0),
54b096ef 102 fVertexCorrelation(0),
51f6de65 103 fVertexCorrelationShift(0),
54b096ef 104 fVertexProfile(0),
1c15d51a 105 fVertexShift(0),
54b096ef 106 fVertexShiftNorm(0),
3d7758c1 107 fEtaCorrelation(0),
51f6de65 108 fEtaCorrelationShift(0),
3d7758c1 109 fEtaProfile(0),
1c15d51a 110 fEtaResolution(0),
69b09e3b 111 fDeltaPhiCorrelation(0),
1c15d51a 112 fpTResolution(0),
113 fEsdTrackCutsPrim(0),
114 fEsdTrackCutsSec(0),
115 fTemp1(0),
116 fTemp2(0),
4ef701a0 117 fMultAll(0),
118 fMultTr(0),
50ec344d 119 fMultVtx(0),
120 fEventStats(0)
0f67a57c 121{
122 //
123 // Constructor. Initialization of pointers
124 //
125
126 // Define input and output slots here
127 DefineInput(0, TChain::Class());
128 DefineOutput(0, TList::Class());
96fcc8a7 129
69b09e3b 130 for (Int_t i=0; i<4; i++)
131 fdNdEtaCorrectionSpecial[i] = 0;
4ef701a0 132
3d7758c1 133 for (Int_t i=0; i<8; i++)
4ef701a0 134 fDeltaPhi[i] = 0;
0f67a57c 135}
136
137AlidNdEtaCorrectionTask::~AlidNdEtaCorrectionTask()
138{
139 //
140 // Destructor
141 //
142
143 // histograms are in the output list and deleted when the output
144 // list is deleted by the TSelector dtor
145
146 if (fOutput) {
147 delete fOutput;
148 fOutput = 0;
149 }
150}
151
152//________________________________________________________________________
153void AlidNdEtaCorrectionTask::ConnectInputData(Option_t *)
154{
155 // Connect ESD
156 // Called once
157
158 Printf("AlidNdEtaCorrectionTask::ConnectInputData called");
159
c17301f3 160 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
161
162 if (!esdH) {
163 Printf("ERROR: Could not get ESDInputHandler");
0f67a57c 164 } else {
c17301f3 165 fESD = esdH->GetEvent();
0f67a57c 166
c17301f3 167 // Enable only the needed branches
168 esdH->SetActiveBranches("AliESDHeader Vertex");
0f67a57c 169
a7f69e56 170 if (fAnalysisMode & AliPWG0Helper::kSPD)
c17301f3 171 esdH->SetActiveBranches("AliESDHeader Vertex AliMultiplicity");
0f67a57c 172
a7f69e56 173 if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS) {
c17301f3 174 esdH->SetActiveBranches("AliESDHeader Vertex Tracks");
0f67a57c 175 }
0f67a57c 176 }
96fcc8a7 177
178 // disable info messages of AliMCEvent (per event)
179 AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
0f67a57c 180}
181
182void AlidNdEtaCorrectionTask::CreateOutputObjects()
183{
184 // create result objects and add to output list
185
186 if (fOption.Contains("only-positive"))
187 {
188 Printf("INFO: Processing only positive particles.");
189 fSignMode = 1;
190 }
191 else if (fOption.Contains("only-negative"))
192 {
193 Printf("INFO: Processing only negative particles.");
194 fSignMode = -1;
195 }
69b09e3b 196
197 if (fOption.Contains("stat-error-1"))
198 {
199 Printf("INFO: Evaluation statistical errors. Mode: 1.");
200 fStatError = 1;
201 }
202 else if (fOption.Contains("stat-error-2"))
203 {
204 Printf("INFO: Evaluation statistical errors. Mode: 2.");
205 fStatError = 2;
206 }
0f67a57c 207
0f67a57c 208 fOutput = new TList;
209 fOutput->SetOwner();
210
770a1f1d 211 fdNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction", fAnalysisMode);
0f67a57c 212 fOutput->Add(fdNdEtaCorrection);
213
1c15d51a 214 fPIDParticles = new TH1F("fPIDParticles", "PID of generated primary particles", 10001, -5000.5, 5000.5);
0f67a57c 215 fOutput->Add(fPIDParticles);
216
1c15d51a 217 fPIDTracks = new TH1F("fPIDTracks", "MC PID of reconstructed tracks", 10001, -5000.5, 5000.5);
0f67a57c 218 fOutput->Add(fPIDTracks);
219
770a1f1d 220 fdNdEtaAnalysisMC = new dNdEtaAnalysis("dndetaMC", "dndetaMC", fAnalysisMode);
0f67a57c 221 fOutput->Add(fdNdEtaAnalysisMC);
222
770a1f1d 223 fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndetaESD", "dndetaESD", fAnalysisMode);
0f67a57c 224 fOutput->Add(fdNdEtaAnalysisESD);
96fcc8a7 225
1c15d51a 226 if (fEsdTrackCuts)
227 {
228 fEsdTrackCutsPrim = dynamic_cast<AliESDtrackCuts*> (fEsdTrackCuts->Clone("fEsdTrackCutsPrim"));
229 fEsdTrackCutsSec = dynamic_cast<AliESDtrackCuts*> (fEsdTrackCuts->Clone("fEsdTrackCutsSec"));
230 fOutput->Add(fEsdTrackCutsPrim);
231 fOutput->Add(fEsdTrackCutsSec);
232 }
233
96fcc8a7 234 if (fOption.Contains("process-types")) {
69b09e3b 235 fdNdEtaCorrectionSpecial[0] = new AlidNdEtaCorrection("dndeta_correction_ND", "dndeta_correction_ND", fAnalysisMode);
236 fdNdEtaCorrectionSpecial[1] = new AlidNdEtaCorrection("dndeta_correction_SD", "dndeta_correction_SD", fAnalysisMode);
237 fdNdEtaCorrectionSpecial[2] = new AlidNdEtaCorrection("dndeta_correction_DD", "dndeta_correction_DD", fAnalysisMode);
96fcc8a7 238
69b09e3b 239 fOutput->Add(fdNdEtaCorrectionSpecial[0]);
240 fOutput->Add(fdNdEtaCorrectionSpecial[1]);
241 fOutput->Add(fdNdEtaCorrectionSpecial[2]);
242 }
243
244 if (fOption.Contains("particle-species")) {
245 fdNdEtaCorrectionSpecial[0] = new AlidNdEtaCorrection("dndeta_correction_pi", "dndeta_correction_pi", fAnalysisMode);
246 fdNdEtaCorrectionSpecial[1] = new AlidNdEtaCorrection("dndeta_correction_K", "dndeta_correction_K", fAnalysisMode);
247 fdNdEtaCorrectionSpecial[2] = new AlidNdEtaCorrection("dndeta_correction_p", "dndeta_correction_p", fAnalysisMode);
248 fdNdEtaCorrectionSpecial[3] = new AlidNdEtaCorrection("dndeta_correction_other", "dndeta_correction_other", fAnalysisMode);
249
250 for (Int_t i=0; i<4; i++)
251 fOutput->Add(fdNdEtaCorrectionSpecial[i]);
96fcc8a7 252 }
253
69b09e3b 254
1c15d51a 255 /*
69b09e3b 256 fTemp1 = new TH2F("fTemp1", "fTemp1", 200, -0.08, 0.08, 200, -0.08, 0.08);
1c15d51a 257 fOutput->Add(fTemp1);
69b09e3b 258 fTemp2 = new TH1F("fTemp2", "fTemp2", 2000, -5, 5);
1c15d51a 259 fOutput->Add(fTemp2);
260 */
54b096ef 261
0fc41645 262 fVertexCorrelation = new TH2F("fVertexCorrelation", "fVertexCorrelation;MC z-vtx;ESD z-vtx", 120, -30, 30, 120, -30, 30);
263 fOutput->Add(fVertexCorrelation);
a7f69e56 264 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 265 fOutput->Add(fVertexCorrelationShift);
54b096ef 266 fVertexProfile = new TProfile("fVertexProfile", "fVertexProfile;MC z-vtx;MC z-vtx - ESD z-vtx", 40, -20, 20);
0fc41645 267 fOutput->Add(fVertexProfile);
1c15d51a 268 fVertexShift = new TH1F("fVertexShift", "fVertexShift;(MC z-vtx - ESD z-vtx);Entries", 201, -2, 2);
269 fOutput->Add(fVertexShift);
a7f69e56 270 fVertexShiftNorm = new TH2F("fVertexShiftNorm", "fVertexShiftNorm;(MC z-vtx - ESD z-vtx) / #sigma_{ESD z-vtx};rec. tracks;Entries", 200, -100, 100, 100, -0.5, 99.5);
0fc41645 271 fOutput->Add(fVertexShiftNorm);
272
3d7758c1 273 fEtaCorrelation = new TH2F("fEtaCorrelation", "fEtaCorrelation;MC #eta;ESD #eta", 120, -3, 3, 120, -3, 3);
0fc41645 274 fOutput->Add(fEtaCorrelation);
69b09e3b 275 fEtaCorrelationShift = new TH2F("fEtaCorrelationShift", "fEtaCorrelationShift;MC #eta;MC #eta - ESD #eta", 120, -3, 3, 100, -0.1, 0.1);
51f6de65 276 fOutput->Add(fEtaCorrelationShift);
3d7758c1 277 fEtaProfile = new TProfile("fEtaProfile", "fEtaProfile;MC #eta;MC #eta - ESD #eta", 120, -3, 3);
0fc41645 278 fOutput->Add(fEtaProfile);
1c15d51a 279 fEtaResolution = new TH1F("fEtaResolution", "fEtaResolution;MC #eta - ESD #eta", 201, -0.2, 0.2);
280 fOutput->Add(fEtaResolution);
281
69b09e3b 282 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 283 fOutput->Add(fpTResolution);
3d7758c1 284
4ef701a0 285 fMultAll = new TH1F("fMultAll", "fMultAll", 500, -0.5, 499.5);
0fc41645 286 fOutput->Add(fMultAll);
4ef701a0 287 fMultVtx = new TH1F("fMultVtx", "fMultVtx", 500, -0.5, 499.5);
0fc41645 288 fOutput->Add(fMultVtx);
4ef701a0 289 fMultTr = new TH1F("fMultTr", "fMultTr", 500, -0.5, 499.5);
0fc41645 290 fOutput->Add(fMultTr);
4ef701a0 291
3d7758c1 292 for (Int_t i=0; i<8; i++)
0fc41645 293 {
69b09e3b 294 fDeltaPhi[i] = new TH2F(Form("fDeltaPhi_%d", i), ";#Delta phi;#phi;Entries", 2000, -0.1, 0.1, 18 * 5, 0, TMath::TwoPi());
0fc41645 295 fOutput->Add(fDeltaPhi[i]);
296 }
50ec344d 297
69b09e3b 298 fEventStats = new TH2F("fEventStats", "fEventStats;event type;status;count", 106, -5.5, 100.5, 4, -0.5, 3.5);
0fc41645 299 fOutput->Add(fEventStats);
69b09e3b 300 fEventStats->GetXaxis()->SetBinLabel(1, "INEL"); // x = -5
301 fEventStats->GetXaxis()->SetBinLabel(2, "NSD"); // x = -4
302 fEventStats->GetXaxis()->SetBinLabel(3, "ND"); // x = -3
303 fEventStats->GetXaxis()->SetBinLabel(4, "SD"); // x = -2
304 fEventStats->GetXaxis()->SetBinLabel(5, "DD"); // x = -1
305
ea441adf 306 for (Int_t i=0; i<100; i++)
69b09e3b 307 fEventStats->GetXaxis()->SetBinLabel(7+i, Form("%d", i));
50ec344d 308
309 fEventStats->GetYaxis()->SetBinLabel(1, "nothing");
310 fEventStats->GetYaxis()->SetBinLabel(2, "trg");
311 fEventStats->GetYaxis()->SetBinLabel(3, "vtx");
312 fEventStats->GetYaxis()->SetBinLabel(4, "trgvtx");
1c15d51a 313
314 if (fEsdTrackCuts)
315 {
316 fEsdTrackCuts->SetName("fEsdTrackCuts");
317 fOutput->Add(fEsdTrackCuts);
318 }
0f67a57c 319}
320
321void AlidNdEtaCorrectionTask::Exec(Option_t*)
322{
323 // process the event
324
325 // Check prerequisites
326 if (!fESD)
327 {
328 AliDebug(AliLog::kError, "ESD branch not available");
329 return;
330 }
331
3d7758c1 332 if (fOnlyPrimaries)
333 Printf("WARNING: Processing only primaries. For systematical studies only!");
69b09e3b 334
335 if (fStatError > 0)
336 Printf("WARNING: Statistical error evaluation active!");
a7f69e56 337
0f67a57c 338 // trigger definition
69b09e3b 339 Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD, fTrigger);
a7f69e56 340 //Printf("Trigger mask: %lld", fESD->GetTriggerMask());
54b096ef 341
3d7758c1 342 if (!eventTriggered)
343 Printf("No trigger");
344
54b096ef 345 // post the data already here
346 PostData(0, fOutput);
69b09e3b 347
54b096ef 348 // MC info
349 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
350 if (!eventHandler) {
351 Printf("ERROR: Could not retrieve MC event handler");
352 return;
353 }
354
355 AliMCEvent* mcEvent = eventHandler->MCEvent();
356 if (!mcEvent) {
357 Printf("ERROR: Could not retrieve MC event");
358 return;
359 }
360
361 AliStack* stack = mcEvent->Stack();
362 if (!stack)
363 {
364 AliDebug(AliLog::kError, "Stack not available");
365 return;
366 }
367
368 AliHeader* header = mcEvent->Header();
369 if (!header)
370 {
371 AliDebug(AliLog::kError, "Header not available");
372 return;
373 }
374
6b62a9c7 375 // get process type;
376 Int_t processType = AliPWG0Helper::GetEventProcessType(header);
377 AliDebug(AliLog::kDebug+1, Form("Found process type %d", processType));
54b096ef 378
6b62a9c7 379 if (processType == AliPWG0Helper::kInvalidProcess)
ea441adf 380 AliDebug(AliLog::kError, "Unknown process type.");
54b096ef 381
382 // get the MC vertex
383 AliGenEventHeader* genHeader = header->GenEventHeader();
384 TArrayF vtxMC(3);
385 genHeader->PrimaryVertex(vtxMC);
386
387 // get the ESD vertex
770a1f1d 388 Bool_t eventVertex = kFALSE;
a7f69e56 389 Double_t vtx[3];
390 const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
391 if (vtxESD && AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
54b096ef 392 {
a7f69e56 393 eventVertex = kTRUE;
54b096ef 394 vtxESD->GetXYZ(vtx);
567160d6 395 }
a7f69e56 396 else
397 vtxESD = 0;
398
50ec344d 399 // fill process type
400 Int_t biny = (Int_t) eventTriggered + 2 * (Int_t) eventVertex;
401 // INEL
69b09e3b 402 fEventStats->Fill(-5, biny);
567160d6 403 // NSD
ea441adf 404 if (processType != AliPWG0Helper::kSD)
69b09e3b 405 fEventStats->Fill(-4, biny);
406 // SD, ND, DD
407 if (processType == AliPWG0Helper::kND)
408 fEventStats->Fill(-3, biny);
409 if (processType == AliPWG0Helper::kSD)
ea441adf 410 fEventStats->Fill(-2, biny);
69b09e3b 411 if (processType == AliPWG0Helper::kDD)
412 fEventStats->Fill(-1, biny);
413
0f67a57c 414 // create list of (label, eta, pt) tuples
415 Int_t inputCount = 0;
416 Int_t* labelArr = 0;
3d7758c1 417 Int_t* labelArr2 = 0; // only for case of SPD
0f67a57c 418 Float_t* etaArr = 0;
69b09e3b 419 Float_t* thirdDimArr = 0;
4ef701a0 420 Float_t* deltaPhiArr = 0;
a7f69e56 421 if (fAnalysisMode & AliPWG0Helper::kSPD)
0f67a57c 422 {
423 // get tracklets
424 const AliMultiplicity* mult = fESD->GetMultiplicity();
425 if (!mult)
426 {
427 AliDebug(AliLog::kError, "AliMultiplicity not available");
428 return;
429 }
430
431 labelArr = new Int_t[mult->GetNumberOfTracklets()];
3d7758c1 432 labelArr2 = new Int_t[mult->GetNumberOfTracklets()];
0f67a57c 433 etaArr = new Float_t[mult->GetNumberOfTracklets()];
69b09e3b 434 thirdDimArr = new Float_t[mult->GetNumberOfTracklets()];
4ef701a0 435 deltaPhiArr = new Float_t[mult->GetNumberOfTracklets()];
0f67a57c 436
1c15d51a 437 // get multiplicity from SPD tracklets
0f67a57c 438 for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
439 {
440 //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
441
69b09e3b 442 Float_t phi = mult->GetPhi(i);
443 if (phi < 0)
444 phi += TMath::Pi() * 2;
4ef701a0 445 Float_t deltaPhi = mult->GetDeltaPhi(i);
4ef701a0 446
447 if (TMath::Abs(deltaPhi) > 1)
448 printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
0f67a57c 449
3d7758c1 450 if (fOnlyPrimaries)
451 if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
452 continue;
453
69b09e3b 454 if (fDeltaPhiCut > 0 && TMath::Abs(deltaPhi) > fDeltaPhiCut)
455 continue;
456
0f67a57c 457 etaArr[inputCount] = mult->GetEta(i);
3d7758c1 458 labelArr[inputCount] = mult->GetLabel(i, 0);
459 labelArr2[inputCount] = mult->GetLabel(i, 1);
69b09e3b 460 thirdDimArr[inputCount] = phi;
4ef701a0 461 deltaPhiArr[inputCount] = deltaPhi;
0f67a57c 462 ++inputCount;
463 }
464 }
a7f69e56 465 else if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
0f67a57c 466 {
467 if (!fEsdTrackCuts)
468 {
469 AliDebug(AliLog::kError, "fESDTrackCuts not available");
470 return;
471 }
a7f69e56 472
69b09e3b 473 if (vtxESD)
1c15d51a 474 {
69b09e3b 475 // get multiplicity from ESD tracks
a7f69e56 476 TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, (fAnalysisMode & AliPWG0Helper::kTPC));
69b09e3b 477 Int_t nGoodTracks = list->GetEntries();
478
479 Printf("Accepted %d tracks", nGoodTracks);
480
481 labelArr = new Int_t[nGoodTracks];
482 labelArr2 = new Int_t[nGoodTracks];
483 etaArr = new Float_t[nGoodTracks];
484 thirdDimArr = new Float_t[nGoodTracks];
485 deltaPhiArr = new Float_t[nGoodTracks];
486
487 // loop over esd tracks
488 for (Int_t i=0; i<nGoodTracks; i++)
1c15d51a 489 {
69b09e3b 490 AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
491 if (!esdTrack)
1c15d51a 492 {
69b09e3b 493 AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
1c15d51a 494 continue;
495 }
69b09e3b 496
a7f69e56 497 //Printf("status is: %u", esdTrack->GetStatus());
498
499 if (fOnlyPrimaries)
500 {
501 Int_t label = TMath::Abs(esdTrack->GetLabel());
502 if (label == 0)
503 continue;
504
505 if (stack->IsPhysicalPrimary(label) == kFALSE)
506 continue;
507 }
69b09e3b 508
509 etaArr[inputCount] = esdTrack->Eta();
510 labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
511 labelArr2[inputCount] = labelArr[inputCount]; // no second label for tracks
512 thirdDimArr[inputCount] = esdTrack->Pt();
513 deltaPhiArr[inputCount] = 0; // no delta phi for tracks
514 ++inputCount;
515 }
1c15d51a 516
69b09e3b 517 delete list;
1c15d51a 518
69b09e3b 519 if (eventTriggered)
520 {
521 // collect values for primaries and secondaries
522 for (Int_t iTrack = 0; iTrack < fESD->GetNumberOfTracks(); iTrack++)
1c15d51a 523 {
69b09e3b 524 AliESDtrack* track = 0;
525
a7f69e56 526 if (fAnalysisMode & AliPWG0Helper::kTPC)
69b09e3b 527 track = AliESDtrackCuts::GetTPCOnlyTrack(fESD, iTrack);
a7f69e56 528 else if (fAnalysisMode & AliPWG0Helper::kTPCITS)
69b09e3b 529 track = fESD->GetTrack(iTrack);
530
531 if (!track)
532 continue;
533
534 Int_t label = TMath::Abs(track->GetLabel());
535 if (!stack->Particle(label) || !stack->Particle(label)->GetPDG())
536 {
537 Printf("WARNING: No particle for %d", label);
538 if (stack->Particle(label))
539 stack->Particle(label)->Print();
540 continue;
541 }
542
543 if (stack->Particle(label)->GetPDG()->Charge() == 0)
544 continue;
545
a7f69e56 546 if (TMath::Abs(track->Eta()) < 0.8 && track->Pt() > 0.15)
69b09e3b 547 {
548 if (stack->IsPhysicalPrimary(label))
1c15d51a 549 {
69b09e3b 550 // primary
551 if (fEsdTrackCutsPrim->AcceptTrack(track))
552 {
a7f69e56 553// if (AliESDtrackCuts::GetSigmaToVertex(track) > 900)
554// {
555// Printf("Track %d has nsigma of %f. Printing track and vertex...", iTrack, AliESDtrackCuts::GetSigmaToVertex(track));
556// Float_t b[2];
557// Float_t r[3];
558// track->GetImpactParameters(b, r);
559// Printf("Impact parameter %f %f and resolution: %f %f %f", b[0], b[1], r[0], r[1], r[2]);
560// track->Print("");
561// if (vtxESD)
562// vtxESD->Print();
563// }
69b09e3b 564 }
1c15d51a 565 }
69b09e3b 566 else
567 {
568 // secondary
569 fEsdTrackCutsSec->AcceptTrack(track);
1c15d51a 570 }
571 }
69b09e3b 572
573 // TODO mem leak in the continue statements above
a7f69e56 574 if (fAnalysisMode & AliPWG0Helper::kTPC)
69b09e3b 575 delete track;
1c15d51a 576 }
1c15d51a 577 }
578 }
0f67a57c 579 }
54b096ef 580 else
0f67a57c 581 return;
0f67a57c 582
0f67a57c 583 // loop over mc particles
584 Int_t nPrim = stack->GetNprimary();
4ef701a0 585 Int_t nAccepted = 0;
0f67a57c 586
587 for (Int_t iMc = 0; iMc < nPrim; ++iMc)
588 {
589 //Printf("Getting particle %d", iMc);
590 TParticle* particle = stack->Particle(iMc);
591
592 if (!particle)
593 continue;
594
595 if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
1c15d51a 596 {
597 //if (TMath::Abs(particle->GetPdgCode()) > 3000 && TMath::Abs(particle->Eta()) < 1.0)
598 // fPIDParticles->Fill(particle->GetPdgCode());
0f67a57c 599 continue;
1c15d51a 600 }
0f67a57c 601
602 if (SignOK(particle->GetPDG()) == kFALSE)
603 continue;
604
1c15d51a 605 if (fPIDParticles && TMath::Abs(particle->Eta()) < 1.0)
606 fPIDParticles->Fill(particle->GetPdgCode());
607
0f67a57c 608 Float_t eta = particle->Eta();
69b09e3b 609
610 Float_t thirdDim = -1;
a7f69e56 611 if (fAnalysisMode & AliPWG0Helper::kSPD)
69b09e3b 612 {
613 if (fFillPhi)
614 {
615 thirdDim = particle->Phi();
616 }
617 else
618 thirdDim = inputCount;
619 }
620 else
621 thirdDim = particle->Pt();
622
623 // calculate y
624 //Float_t y = 0.5 * TMath::Log((particle->Energy() + particle->Pz()) / (particle->Energy() - particle->Pz()));
625 //fTemp1->Fill(eta);
626 //fTemp2->Fill(y);
0f67a57c 627
51f6de65 628 fdNdEtaCorrection->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType);
0f67a57c 629
69b09e3b 630 if (fOption.Contains("process-types"))
96fcc8a7 631 {
632 // non diffractive
6b62a9c7 633 if (processType==AliPWG0Helper::kND)
69b09e3b 634 fdNdEtaCorrectionSpecial[0]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType);
96fcc8a7 635
636 // single diffractive
6b62a9c7 637 if (processType==AliPWG0Helper::kSD)
69b09e3b 638 fdNdEtaCorrectionSpecial[1]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType);
96fcc8a7 639
640 // double diffractive
6b62a9c7 641 if (processType==AliPWG0Helper::kDD)
69b09e3b 642 fdNdEtaCorrectionSpecial[2]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType);
643 }
644
645 if (fOption.Contains("particle-species"))
646 {
647 Int_t id = -1;
648 switch (TMath::Abs(particle->GetPdgCode()))
649 {
650 case 211: id = 0; break;
651 case 321: id = 1; break;
652 case 2212: id = 2; break;
653 default: id = 3; break;
654 }
655 fdNdEtaCorrectionSpecial[id]->FillMCParticle(vtxMC[2], eta, thirdDim, eventTriggered, eventVertex, processType);
96fcc8a7 656 }
657
0f67a57c 658 if (eventTriggered)
659 if (eventVertex)
51f6de65 660 fdNdEtaAnalysisMC->FillTrack(vtxMC[2], eta, thirdDim);
4ef701a0 661
0fc41645 662 // TODO this value might be needed lower for the SPD study (only used for control histograms anyway)
ea441adf 663 if (TMath::Abs(eta) < 1.5) // && pt > 0.2)
4ef701a0 664 nAccepted++;
665 }
666
69b09e3b 667 fEventStats->Fill(AliPWG0Helper::GetLastProcessType(), biny);
ea441adf 668
4ef701a0 669 fMultAll->Fill(nAccepted);
670 if (eventTriggered) {
671 fMultTr->Fill(nAccepted);
672 if (eventVertex)
673 fMultVtx->Fill(nAccepted);
0f67a57c 674 }
675
3d7758c1 676 Int_t processed = 0;
69b09e3b 677
678 Bool_t* primCount = 0;
679 if (fStatError > 0)
680 {
681 primCount = new Bool_t[nPrim];
682 for (Int_t i=0; i<nPrim; ++i)
683 primCount[i] = kFALSE;
684 }
3d7758c1 685
0f67a57c 686 for (Int_t i=0; i<inputCount; ++i)
687 {
688 Int_t label = labelArr[i];
3d7758c1 689 Int_t label2 = labelArr2[i];
0f67a57c 690
691 if (label < 0)
692 {
0fc41645 693 Printf("WARNING: cannot find corresponding mc particle for track(let) %d with label %d.", i, label);
0f67a57c 694 continue;
695 }
696
697 TParticle* particle = stack->Particle(label);
698 if (!particle)
699 {
700 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
701 continue;
702 }
1c15d51a 703
704 fPIDTracks->Fill(particle->GetPdgCode());
3d7758c1 705
706 // find particle that is filled in the correction map
707 // this should be the particle that has been reconstructed
708 // for TPC this is clear and always identified by the track's label
709 // for SPD the labels might not be identical. In this case the particle closest to the beam line that is a primary is taken:
710 // L1 L2 (P = primary, S = secondary)
711 // P P' : --> P
712 // P S : --> P
713 // S P : --> P
714 // S S' : --> S
715
716 if (label != label2 && label2 >= 0)
717 {
718 if (stack->IsPhysicalPrimary(label) == kFALSE && stack->IsPhysicalPrimary(label2))
719 {
720 particle = stack->Particle(label2);
721 if (!particle)
722 {
723 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label2));
724 continue;
725 }
726 }
727 }
0f67a57c 728
729 if (eventTriggered && eventVertex)
730 {
3d7758c1 731 if (SignOK(particle->GetPDG()) == kFALSE)
732 continue;
733
734 processed++;
735
1c15d51a 736 // resolutions
737 fEtaResolution->Fill(particle->Eta() - etaArr[i]);
69b09e3b 738
a7f69e56 739 if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
69b09e3b 740 if (TMath::Abs(particle->Eta() < 0.9) && particle->Pt() > 0)
741 fpTResolution->Fill(particle->Pt(), (particle->Pt() - thirdDimArr[i]) / particle->Pt());
1c15d51a 742
51f6de65 743 Float_t eta = -999;
744 Float_t thirdDim = -1;
745
3d7758c1 746 Bool_t firstIsPrim = stack->IsPhysicalPrimary(label);
a7f69e56 747 // in case of same label the MC values are filled, otherwise (background) the reconstructed values
748 if (label == label2)
3d7758c1 749 {
51f6de65 750 eta = particle->Eta();
69b09e3b 751
a7f69e56 752 if (fAnalysisMode & AliPWG0Helper::kSPD)
69b09e3b 753 {
754 if (fFillPhi)
755 {
756 thirdDim = particle->Phi();
757 }
758 else
759 thirdDim = inputCount;
760 }
761 else
762 thirdDim = particle->Pt();
3d7758c1 763 }
764 else
765 {
a7f69e56 766 if (fAnalysisMode & AliPWG0Helper::kSPD && !fFillPhi)
69b09e3b 767 {
768 thirdDim = inputCount;
769 }
770 else
771 thirdDim = thirdDimArr[i];
772
51f6de65 773 eta = etaArr[i];
3d7758c1 774 }
745d6088 775
69b09e3b 776 Bool_t fillTrack = kTRUE;
51f6de65 777
69b09e3b 778 // statistical error evaluation active?
779 if (fStatError > 0)
780 {
781 Bool_t statErrorDecision = kFALSE;
782
783 // primary and not yet counted
784 if (label == label2 && firstIsPrim && !primCount[label])
785 {
786 statErrorDecision = kTRUE;
787 primCount[label] = kTRUE;
788 }
789
790 // in case of 1 we count only unique primaries, in case of 2 all the rest
791 if (fStatError == 1)
792 {
793 fillTrack = statErrorDecision;
794 }
795 else if (fStatError == 2)
796 fillTrack = !statErrorDecision;
797 }
798
799 if (fillTrack)
800 fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], eta, thirdDim);
51f6de65 801
802 // eta comparison for tracklets with the same label (others are background)
803 if (label == label2)
804 {
805 fEtaProfile->Fill(particle->Eta(), particle->Eta() - etaArr[i]);
806 fEtaCorrelation->Fill(etaArr[i], particle->Eta());
69b09e3b 807 fEtaCorrelationShift->Fill(particle->Eta(), particle->Eta() - etaArr[i]);
51f6de65 808 }
3d7758c1 809
51f6de65 810 fdNdEtaAnalysisESD->FillTrack(vtxMC[2], particle->Eta(), thirdDim);
3d7758c1 811
69b09e3b 812 if (fOption.Contains("process-types"))
96fcc8a7 813 {
814 // non diffractive
6b62a9c7 815 if (processType == AliPWG0Helper::kND)
69b09e3b 816 fdNdEtaCorrectionSpecial[0]->FillTrackedParticle(vtxMC[2], eta, thirdDim);
96fcc8a7 817
818 // single diffractive
1c15d51a 819 if (processType == AliPWG0Helper::kSD)
69b09e3b 820 fdNdEtaCorrectionSpecial[1]->FillTrackedParticle(vtxMC[2], eta, thirdDim);
96fcc8a7 821
822 // double diffractive
6b62a9c7 823 if (processType == AliPWG0Helper::kDD)
69b09e3b 824 fdNdEtaCorrectionSpecial[2]->FillTrackedParticle(vtxMC[2], eta, thirdDim);
825 }
826
827 if (fOption.Contains("particle-species"))
828 {
829 // find mother first
830 TParticle* mother = AliPWG0Helper::FindPrimaryMother(stack, label);
831
832 Int_t id = -1;
833 switch (TMath::Abs(mother->GetPdgCode()))
834 {
835 case 211: id = 0; break;
836 case 321: id = 1; break;
837 case 2212: id = 2; break;
838 default: id = 3; break;
839 }
840 fdNdEtaCorrectionSpecial[id]->FillTrackedParticle(vtxMC[2], eta, thirdDim);
96fcc8a7 841 }
4ef701a0 842
745d6088 843 // control histograms
3d7758c1 844 Int_t hist = -1;
845 if (label == label2)
4ef701a0 846 {
69b09e3b 847 if (firstIsPrim)
3d7758c1 848 {
69b09e3b 849 hist = 0;
3d7758c1 850 }
851 else
852 hist = 1;
853 }
854 else if (label2 >= 0)
855 {
856 Bool_t secondIsPrim = stack->IsPhysicalPrimary(label2);
857 if (firstIsPrim && secondIsPrim)
858 {
859 hist = 2;
860 }
861 else if (firstIsPrim && !secondIsPrim)
862 {
863 hist = 3;
69b09e3b 864
3d7758c1 865 // check if secondary is caused by the primary or it is a fake combination
866 //Printf("PS case --> Label 1 is %d, label 2 is %d", label, label2);
867 Int_t mother = label2;
868 while (stack->Particle(mother) && stack->Particle(mother)->GetMother(0) >= 0)
869 {
870 //Printf(" %d created by %d, %d", mother, stack->Particle(mother)->GetMother(0), stack->Particle(mother)->GetMother(1));
871 if (stack->Particle(mother)->GetMother(0) == label)
872 {
873 /*Printf("The untraceable decay was:");
874 Printf(" from:");
875 particle->Print();
876 Printf(" to (status code %d):", stack->Particle(mother)->GetStatusCode());
877 stack->Particle(mother)->Print();*/
878 hist = 4;
879 }
880 mother = stack->Particle(mother)->GetMother(0);
881 }
882 }
883 else if (!firstIsPrim && secondIsPrim)
884 {
885 hist = 5;
886 }
887 else if (!firstIsPrim && !secondIsPrim)
888 {
889 hist = 6;
890 }
891
4ef701a0 892 }
893 else
3d7758c1 894 hist = 7;
895
69b09e3b 896 fDeltaPhi[hist]->Fill(deltaPhiArr[i], thirdDimArr[i]);
0f67a57c 897 }
898 }
69b09e3b 899
900 if (primCount)
901 {
902 delete[] primCount;
903 primCount = 0;
904 }
0f67a57c 905
3d7758c1 906 if (processed < inputCount)
907 Printf("Only %d out of %d track(let)s were used", processed, inputCount);
908
0f67a57c 909 if (eventTriggered && eventVertex)
96fcc8a7 910 {
a7f69e56 911 Double_t diff = vtxMC[2] - vtx[2];
912 fVertexShift->Fill(diff);
913 if (vtxESD->GetZRes() > 0)
914 fVertexShiftNorm->Fill(diff / vtxESD->GetZRes(), inputCount);
915
916 fVertexCorrelation->Fill(vtxMC[2], vtx[2]);
917 fVertexCorrelationShift->Fill(vtxMC[2], vtxMC[2] - vtx[2], inputCount);
918 fVertexProfile->Fill(vtxMC[2], vtxMC[2] - vtx[2]);
919 }
920
921 if (eventTriggered && eventVertex)
922 {
0f67a57c 923 fdNdEtaAnalysisMC->FillEvent(vtxMC[2], inputCount);
96fcc8a7 924 fdNdEtaAnalysisESD->FillEvent(vtxMC[2], inputCount);
925 }
0f67a57c 926
96fcc8a7 927 // stuff regarding the vertex reco correction and trigger bias correction
0f67a57c 928 fdNdEtaCorrection->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
96fcc8a7 929
69b09e3b 930 if (fOption.Contains("process-types"))
96fcc8a7 931 {
932 // non diffractive
a7f69e56 933 if (processType == AliPWG0Helper::kND)
69b09e3b 934 fdNdEtaCorrectionSpecial[0]->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
96fcc8a7 935
936 // single diffractive
6b62a9c7 937 if (processType == AliPWG0Helper::kSD)
69b09e3b 938 fdNdEtaCorrectionSpecial[1]->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
96fcc8a7 939
940 // double diffractive
6b62a9c7 941 if (processType == AliPWG0Helper::kDD)
69b09e3b 942 fdNdEtaCorrectionSpecial[2]->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
0f67a57c 943 }
69b09e3b 944
945 if (fOption.Contains("particle-species"))
946 for (Int_t id=0; id<4; id++)
947 fdNdEtaCorrectionSpecial[id]->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
948
949 if (etaArr)
950 delete[] etaArr;
951 if (labelArr)
952 delete[] labelArr;
953 if (labelArr2)
954 delete[] labelArr2;
955 if (thirdDimArr)
956 delete[] thirdDimArr;
957 if (deltaPhiArr)
958 delete[] deltaPhiArr;
0f67a57c 959}
960
961void AlidNdEtaCorrectionTask::Terminate(Option_t *)
962{
963 // The Terminate() function is the last function to be called during
964 // a query. It always runs on the client, it can be used to present
965 // the results graphically or save the results to file.
966
967 fOutput = dynamic_cast<TList*> (GetOutputData(0));
968 if (!fOutput) {
969 Printf("ERROR: fOutput not available");
970 return;
971 }
972
973 fdNdEtaCorrection = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction"));
974 fdNdEtaAnalysisMC = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaMC"));
975 fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaESD"));
976 if (!fdNdEtaCorrection || !fdNdEtaAnalysisMC || !fdNdEtaAnalysisESD)
977 {
978 AliDebug(AliLog::kError, "Could not read object from output list");
979 return;
980 }
981
982 fdNdEtaCorrection->Finish();
983
984 TString fileName;
985 fileName.Form("correction_map%s.root", fOption.Data());
986 TFile* fout = new TFile(fileName, "RECREATE");
987
1c15d51a 988 fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCuts"));
0f67a57c 989 if (fEsdTrackCuts)
990 fEsdTrackCuts->SaveHistograms("esd_track_cuts");
1c15d51a 991
992 fEsdTrackCutsPrim = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCutsPrim"));
993 if (fEsdTrackCutsPrim)
994 fEsdTrackCutsPrim->SaveHistograms("esd_track_cuts_primaries");
995
996 fEsdTrackCutsSec = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCutsSec"));
997 if (fEsdTrackCutsSec)
998 fEsdTrackCutsSec->SaveHistograms("esd_track_cuts_secondaries");
999
0f67a57c 1000 fdNdEtaCorrection->SaveHistograms();
1001 fdNdEtaAnalysisMC->SaveHistograms();
1002 fdNdEtaAnalysisESD->SaveHistograms();
1003
69b09e3b 1004 if (fOutput->FindObject("dndeta_correction_ND"))
1005 {
1006 fdNdEtaCorrectionSpecial[0] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_ND"));
1007 fdNdEtaCorrectionSpecial[1] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_SD"));
1008 fdNdEtaCorrectionSpecial[2] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_DD"));
1009 }
1010 else
1011 {
1012 fdNdEtaCorrectionSpecial[0] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_pi"));
1013 fdNdEtaCorrectionSpecial[1] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_K"));
1014 fdNdEtaCorrectionSpecial[2] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_p"));
1015 fdNdEtaCorrectionSpecial[3] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_other"));
1016 }
1017 for (Int_t i=0; i<4; ++i)
1018 if (fdNdEtaCorrectionSpecial[i])
1019 fdNdEtaCorrectionSpecial[i]->SaveHistograms();
96fcc8a7 1020
1c15d51a 1021 fTemp1 = dynamic_cast<TH1*> (fOutput->FindObject("fTemp1"));
1022 if (fTemp1)
1023 fTemp1->Write();
1024 fTemp2 = dynamic_cast<TH1*> (fOutput->FindObject("fTemp2"));
1025 if (fTemp2)
1026 fTemp2->Write();
96fcc8a7 1027
0fc41645 1028 fVertexCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fVertexCorrelation"));
54b096ef 1029 if (fVertexCorrelation)
1030 fVertexCorrelation->Write();
a7f69e56 1031 fVertexCorrelationShift = dynamic_cast<TH3F*> (fOutput->FindObject("fVertexCorrelationShift"));
51f6de65 1032 if (fVertexCorrelationShift)
a7f69e56 1033 {
1034 ((TH2*) fVertexCorrelationShift->Project3D("yx"))->FitSlicesY();
51f6de65 1035 fVertexCorrelationShift->Write();
a7f69e56 1036 }
0fc41645 1037 fVertexProfile = dynamic_cast<TProfile*> (fOutput->FindObject("fVertexProfile"));
54b096ef 1038 if (fVertexProfile)
1039 fVertexProfile->Write();
1c15d51a 1040 fVertexShift = dynamic_cast<TH1F*> (fOutput->FindObject("fVertexShift"));
1041 if (fVertexShift)
1042 fVertexShift->Write();
a7f69e56 1043 fVertexShiftNorm = dynamic_cast<TH2F*> (fOutput->FindObject("fVertexShiftNorm"));
54b096ef 1044 if (fVertexShiftNorm)
a7f69e56 1045 {
1046 fVertexShiftNorm->ProjectionX();
54b096ef 1047 fVertexShiftNorm->Write();
a7f69e56 1048 }
3d7758c1 1049
0fc41645 1050 fEtaCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelation"));
3d7758c1 1051 if (fEtaCorrelation)
1052 fEtaCorrelation->Write();
51f6de65 1053 fEtaCorrelationShift = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelationShift"));
1054 if (fEtaCorrelationShift)
a7f69e56 1055 {
1056 fEtaCorrelationShift->FitSlicesY();
51f6de65 1057 fEtaCorrelationShift->Write();
a7f69e56 1058 }
0fc41645 1059 fEtaProfile = dynamic_cast<TProfile*> (fOutput->FindObject("fEtaProfile"));
3d7758c1 1060 if (fEtaProfile)
1061 fEtaProfile->Write();
1c15d51a 1062 fEtaResolution = dynamic_cast<TH1F*> (fOutput->FindObject("fEtaResolution"));
1063 if (fEtaResolution)
1064 fEtaResolution->Write();
69b09e3b 1065 fpTResolution = dynamic_cast<TH2F*> (fOutput->FindObject("fpTResolution"));
1c15d51a 1066 if (fpTResolution)
69b09e3b 1067 {
1068 fpTResolution->FitSlicesY();
1c15d51a 1069 fpTResolution->Write();
69b09e3b 1070 }
3d7758c1 1071
0fc41645 1072 fMultAll = dynamic_cast<TH1F*> (fOutput->FindObject("fMultAll"));
4ef701a0 1073 if (fMultAll)
1074 fMultAll->Write();
1075
0fc41645 1076 fMultTr = dynamic_cast<TH1F*> (fOutput->FindObject("fMultTr"));
4ef701a0 1077 if (fMultTr)
1078 fMultTr->Write();
0fc41645 1079
1080 fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx"));
4ef701a0 1081 if (fMultVtx)
1082 fMultVtx->Write();
1083
3d7758c1 1084 for (Int_t i=0; i<8; ++i)
0fc41645 1085 {
69b09e3b 1086 fDeltaPhi[i] = dynamic_cast<TH2*> (fOutput->FindObject(Form("fDeltaPhi_%d", i)));
4ef701a0 1087 if (fDeltaPhi[i])
1088 fDeltaPhi[i]->Write();
0fc41645 1089 }
54b096ef 1090
0fc41645 1091 fEventStats = dynamic_cast<TH2F*> (fOutput->FindObject("fEventStats"));
50ec344d 1092 if (fEventStats)
1093 fEventStats->Write();
1094
1c15d51a 1095 fPIDParticles = dynamic_cast<TH1F*> (fOutput->FindObject("fPIDParticles"));
1096 if (fPIDParticles)
1097 fPIDParticles->Write();
1098
1099 fPIDTracks = dynamic_cast<TH1F*> (fOutput->FindObject("fPIDTracks"));
1100 if (fPIDTracks)
1101 fPIDTracks->Write();
0f67a57c 1102
1c15d51a 1103 //fdNdEtaCorrection->DrawHistograms();
0f67a57c 1104
96fcc8a7 1105 Printf("Writting result to %s", fileName.Data());
1106
0f67a57c 1107 if (fPIDParticles && fPIDTracks)
1108 {
0f67a57c 1109 TDatabasePDG* pdgDB = new TDatabasePDG;
1110
1111 for (Int_t i=0; i <= fPIDParticles->GetNbinsX()+1; ++i)
1112 if (fPIDParticles->GetBinContent(i) > 0)
1c15d51a 1113 {
1114 TObject* pdgParticle = pdgDB->GetParticle((Int_t) fPIDParticles->GetBinCenter(i));
1115 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));
1116 }
0f67a57c 1117
1118 delete pdgDB;
1119 pdgDB = 0;
1120 }
1c15d51a 1121
1122 fout->Write();
1123 fout->Close();
0f67a57c 1124}
1125
1126Bool_t AlidNdEtaCorrectionTask::SignOK(TParticlePDG* particle)
1127{
1128 // returns if a particle with this sign should be counted
1129 // this is determined by the value of fSignMode, which should have the same sign
1130 // as the charge
1131 // if fSignMode is 0 all particles are counted
1132
1133 if (fSignMode == 0)
1134 return kTRUE;
1135
1136 if (!particle)
1137 {
1138 printf("WARNING: not counting a particle that does not have a pdg particle\n");
1139 return kFALSE;
1140 }
1141
1142 Double_t charge = particle->Charge();
1143
1144 if (fSignMode > 0)
1145 if (charge < 0)
1146 return kFALSE;
1147
1148 if (fSignMode < 0)
1149 if (charge > 0)
1150 return kFALSE;
1151
1152 return kTRUE;
1153}
96fcc8a7 1154