]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG0/dNdEta/AlidNdEtaTask.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
bug fix
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AlidNdEtaTask.cxx
CommitLineData
0f67a57c 1/* $Id$ */
2
3#include "AlidNdEtaTask.h"
4
5#include <TStyle.h>
6#include <TSystem.h>
7#include <TCanvas.h>
8#include <TVector3.h>
9#include <TChain.h>
10#include <TFile.h>
11#include <TH1F.h>
68fa248f 12#include <TH1D.h>
0f67a57c 13#include <TH2F.h>
14#include <TH3F.h>
15#include <TParticle.h>
16#include <TRandom.h>
17#include <TNtuple.h>
18#include <TObjString.h>
19#include <TF1.h>
69b09e3b 20#include <TGraph.h>
0f67a57c 21
22#include <AliLog.h>
23#include <AliESDVertex.h>
24#include <AliESDEvent.h>
25#include <AliStack.h>
26#include <AliHeader.h>
27#include <AliGenEventHeader.h>
28#include <AliMultiplicity.h>
29#include <AliAnalysisManager.h>
30#include <AliMCEventHandler.h>
31#include <AliMCEvent.h>
32#include <AliESDInputHandler.h>
1d107532 33#include <AliESDInputHandlerRP.h>
69b09e3b 34#include <AliESDHeader.h>
0f67a57c 35
745d6088 36#include "AliESDtrackCuts.h"
0f67a57c 37#include "AliPWG0Helper.h"
38#include "AliCorrection.h"
39#include "AliCorrectionMatrix3D.h"
40#include "dNdEta/dNdEtaAnalysis.h"
70fdd197 41#include "AliTriggerAnalysis.h"
12bb57f1 42#include "AliPhysicsSelection.h"
43
44//#define FULLALIROOT
45
46#ifdef FULLALIROOT
47 #include "../ITS/AliITSRecPoint.h"
48 #include "AliCDBManager.h"
49 #include "AliCDBEntry.h"
50 #include "AliGeomManager.h"
51 #include "TGeoManager.h"
52#endif
0f67a57c 53
54ClassImp(AlidNdEtaTask)
55
56AlidNdEtaTask::AlidNdEtaTask(const char* opt) :
76532b17 57 AliAnalysisTaskSE("AlidNdEtaTask"),
0f67a57c 58 fESD(0),
59 fOutput(0),
60 fOption(opt),
a7f69e56 61 fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn)),
70fdd197 62 fTrigger(AliTriggerAnalysis::kMB1),
69b09e3b 63 fFillPhi(kFALSE),
64 fDeltaPhiCut(-1),
0f67a57c 65 fReadMC(kFALSE),
54b096ef 66 fUseMCVertex(kFALSE),
c17301f3 67 fOnlyPrimaries(kFALSE),
3d7758c1 68 fUseMCKine(kFALSE),
1d107532 69 fSymmetrize(kFALSE),
81be4ee8 70 fMultAxisEta1(kFALSE),
71 fDiffTreatment(AliPWG0Helper::kMCFlags),
0f67a57c 72 fEsdTrackCuts(0),
73 fdNdEtaAnalysisESD(0),
74 fMult(0),
75 fMultVtx(0),
76 fEvents(0),
54b096ef 77 fVertexResolution(0),
0f67a57c 78 fdNdEtaAnalysis(0),
69b09e3b 79 fdNdEtaAnalysisND(0),
0fc41645 80 fdNdEtaAnalysisNSD(0),
81be4ee8 81 fdNdEtaAnalysisOnePart(0),
0f67a57c 82 fdNdEtaAnalysisTr(0),
83 fdNdEtaAnalysisTrVtx(0),
84 fdNdEtaAnalysisTracks(0),
54b096ef 85 fPartPt(0),
51f6de65 86 fVertex(0),
1d107532 87 fVertexVsMult(0),
ea441adf 88 fPhi(0),
69b09e3b 89 fRawPt(0),
ea441adf 90 fEtaPhi(0),
12bb57f1 91 fModuleMap(0),
69b09e3b 92 fDeltaPhi(0),
a7f69e56 93 fDeltaTheta(0),
69b09e3b 94 fFiredChips(0),
1d107532 95 fTrackletsVsClusters(0),
96 fTrackletsVsUnassigned(0),
1d107532 97 fStats(0),
68fa248f 98 fStats2(0),
99 fPtMin(0.15),
100 fEta(0x0),
101 fEtaMC(0x0),
102 fHistEvents(0),
103 fHistEventsMC(0),
104 fTrigEffNum(0),
105 fTrigEffDen(0),
106 fVtxEffNum(0),
107 fVtxEffDen(0),
108 fVtxTrigEffNum(0),
109 fVtxTrigEffDen(0)
0f67a57c 110{
111 //
112 // Constructor. Initialization of pointers
113 //
114
115 // Define input and output slots here
76532b17 116 //DefineInput(0, TChain::Class());
117 DefineOutput(1, TList::Class());
69b09e3b 118
119 fZPhi[0] = 0;
120 fZPhi[1] = 0;
567160d6 121
122 AliLog::SetClassDebugLevel("AlidNdEtaTask", AliLog::kWarning);
0f67a57c 123}
124
125AlidNdEtaTask::~AlidNdEtaTask()
126{
127 //
128 // Destructor
129 //
130
131 // histograms are in the output list and deleted when the output
132 // list is deleted by the TSelector dtor
133
1103419a 134 if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
0f67a57c 135 delete fOutput;
136 fOutput = 0;
137 }
138}
139
76532b17 140Bool_t AlidNdEtaTask::UserNotify()
567160d6 141{
142 static Int_t count = 0;
143 count++;
69b09e3b 144 Printf("Processing %d. file: %s", count, ((TTree*) GetInputData(0))->GetCurrentFile()->GetName());
567160d6 145 return kTRUE;
146}
147
0f67a57c 148//________________________________________________________________________
76532b17 149void AlidNdEtaTask::ConnectInputData(Option_t *opt)
0f67a57c 150{
151 // Connect ESD
152 // Called once
76532b17 153
154 AliAnalysisTaskSE::ConnectInputData(opt);
0f67a57c 155
156 Printf("AlidNdEtaTask::ConnectInputData called");
157
76532b17 158 /*
567160d6 159 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
160
161 if (!esdH) {
162 Printf("ERROR: Could not get ESDInputHandler");
0f67a57c 163 } else {
567160d6 164 fESD = esdH->GetEvent();
a7f69e56 165
166 TString branches("AliESDHeader Vertex ");
0f67a57c 167
70fdd197 168 if (fAnalysisMode & AliPWG0Helper::kSPD || fTrigger & AliTriggerAnalysis::kOfflineFlag)
a7f69e56 169 branches += "AliMultiplicity ";
170
171 if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
172 branches += "Tracks ";
173
c17301f3 174 // Enable only the needed branches
a7f69e56 175 esdH->SetActiveBranches(branches);
0f67a57c 176 }
76532b17 177 */
54b096ef 178
179 // disable info messages of AliMCEvent (per event)
180 AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
12bb57f1 181
182 #ifdef FULLALIROOT
02e1bcdb 183 AliCDBManager::Instance()->SetDefaultStorage("raw://");
184 //AliCDBManager::Instance()->SetDefaultStorage("MC", "Residual");
12bb57f1 185 AliCDBManager::Instance()->SetRun(0);
186
187 AliCDBManager* mgr = AliCDBManager::Instance();
188 AliCDBEntry* obj = mgr->Get(AliCDBPath("GRP", "Geometry", "Data"));
189 AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
190
191 AliGeomManager::GetNalignable("ITS");
192 AliGeomManager::ApplyAlignObjsFromCDB("ITS");
193 #endif
0f67a57c 194}
195
76532b17 196void AlidNdEtaTask::UserCreateOutputObjects()
0f67a57c 197{
198 // create result objects and add to output list
199
54b096ef 200 Printf("AlidNdEtaTask::CreateOutputObjects");
1103419a 201 //AliLog::SetClassDebugLevel("AliPhysicsSelection", AliLog::kDebug);
54b096ef 202
5a6310fe 203 if (fOnlyPrimaries)
204 Printf("WARNING: Processing only primaries (MC information used). This is for systematical checks only.");
205
206 if (fUseMCKine)
207 Printf("WARNING: Using MC kine information. This is for systematical checks only.");
208
209 if (fUseMCVertex)
210 Printf("WARNING: Replacing vertex by MC vertex. This is for systematical checks only.");
211
0f67a57c 212 fOutput = new TList;
213 fOutput->SetOwner();
214
770a1f1d 215 fdNdEtaAnalysisESD = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD", fAnalysisMode);
0f67a57c 216 fOutput->Add(fdNdEtaAnalysisESD);
217
218 fMult = new TH1F("fMult", "fMult;Ntracks;Count", 201, -0.5, 200.5);
219 fOutput->Add(fMult);
220
221 fMultVtx = new TH1F("fMultVtx", "fMultVtx;Ntracks;Count", 201, -0.5, 200.5);
222 fOutput->Add(fMultVtx);
223
224 for (Int_t i=0; i<3; ++i)
225 {
81be4ee8 226 fPartEta[i] = new TH1F(Form("dndeta_check_%d", i), Form("dndeta_check_%d", i), 60, -3, 3);
0f67a57c 227 fPartEta[i]->Sumw2();
228 fOutput->Add(fPartEta[i]);
229 }
230
51f6de65 231 fEvents = new TH1F("dndeta_check_vertex", "dndeta_check_vertex", 800, -40, 40);
0f67a57c 232 fOutput->Add(fEvents);
233
81be4ee8 234 fVertexResolution = new TH2F("dndeta_vertex_resolution_z", ";z resolution;#Delta #phi", 1000, 0, 2, 100, 0, 0.2);
54b096ef 235 fOutput->Add(fVertexResolution);
236
ea441adf 237 fPhi = new TH1F("fPhi", "fPhi;#phi in rad.;count", 720, 0, 2 * TMath::Pi());
238 fOutput->Add(fPhi);
239
81be4ee8 240 fEtaPhi = new TH2F("fEtaPhi", "fEtaPhi;#eta;#phi in rad.;count", 80, -4, 4, 18*20, 0, 2 * TMath::Pi());
ea441adf 241 fOutput->Add(fEtaPhi);
12bb57f1 242
1d107532 243 fStats = new TH1F("fStats", "fStats", 5, 0.5, 5.5);
69b09e3b 244 fStats->GetXaxis()->SetBinLabel(1, "vertexer 3d");
245 fStats->GetXaxis()->SetBinLabel(2, "vertexer z");
a7f69e56 246 fStats->GetXaxis()->SetBinLabel(3, "trigger");
1d107532 247 fStats->GetXaxis()->SetBinLabel(4, "physics events");
248 fStats->GetXaxis()->SetBinLabel(5, "physics events after veto");
69b09e3b 249 fOutput->Add(fStats);
1d107532 250
76532b17 251 fStats2 = new TH2F("fStats2", "fStats2", 7, -0.5, 6.5, 12, -0.5, 11.5);
12bb57f1 252
253 fStats2->GetXaxis()->SetBinLabel(1, "No trigger");
81be4ee8 254 fStats2->GetXaxis()->SetBinLabel(2, "No Vertex");
255 fStats2->GetXaxis()->SetBinLabel(3, "Vtx res. too bad");
256 fStats2->GetXaxis()->SetBinLabel(4, "|z-vtx| > 15");
257 fStats2->GetXaxis()->SetBinLabel(5, "|z-vtx| > 10");
258 fStats2->GetXaxis()->SetBinLabel(6, "0 tracklets");
12bb57f1 259 fStats2->GetXaxis()->SetBinLabel(7, "Selected");
260
1d107532 261 fStats2->GetYaxis()->SetBinLabel(1, "n/a");
12bb57f1 262 fStats2->GetYaxis()->SetBinLabel(2, "empty");
263 fStats2->GetYaxis()->SetBinLabel(3, "BBA");
264 fStats2->GetYaxis()->SetBinLabel(4, "BBC");
265 fStats2->GetYaxis()->SetBinLabel(5, "BBA BBC");
1d107532 266 fStats2->GetYaxis()->SetBinLabel(6, "BGA");
267 fStats2->GetYaxis()->SetBinLabel(7, "BGC");
12bb57f1 268 fStats2->GetYaxis()->SetBinLabel(8, "BGA BGC");
269 fStats2->GetYaxis()->SetBinLabel(9, "BBA BGC");
270 fStats2->GetYaxis()->SetBinLabel(10, "BGA BBC");
76532b17 271 fStats2->GetYaxis()->SetBinLabel(11, "GFO");
272 fStats2->GetYaxis()->SetBinLabel(12, "NO GFO");
1d107532 273 fOutput->Add(fStats2);
69b09e3b 274
12bb57f1 275 fTrackletsVsClusters = new TH2F("fTrackletsVsClusters", ";tracklets;clusters in ITS", 50, -0.5, 49.5, 1000, -0.5, 999.5);
276 fOutput->Add(fTrackletsVsClusters);
277
a7f69e56 278 if (fAnalysisMode & AliPWG0Helper::kSPD)
69b09e3b 279 {
81be4ee8 280 fDeltaPhi = new TH1F("fDeltaPhi", "fDeltaPhi;#Delta #phi;Entries", 500, -0.16, 0.16);
69b09e3b 281 fOutput->Add(fDeltaPhi);
81be4ee8 282 fDeltaTheta = new TH1F("fDeltaTheta", "fDeltaTheta;#Delta #theta;Entries", 500, -0.05, 0.05);
a7f69e56 283 fOutput->Add(fDeltaTheta);
69b09e3b 284 fFiredChips = new TH2F("fFiredChips", "fFiredChips;Chips L1 + L2;tracklets", 1201, -0.5, 1201, 50, -0.5, 49.5);
285 fOutput->Add(fFiredChips);
1d107532 286 fTrackletsVsUnassigned = new TH2F("fTrackletsVsUnassigned", ";tracklets;unassigned clusters in L0", 50, -0.5, 49.5, 200, -0.5, 199.5);
287 fOutput->Add(fTrackletsVsUnassigned);
69b09e3b 288 for (Int_t i=0; i<2; i++)
289 {
12bb57f1 290 fZPhi[i] = new TH2F(Form("fZPhi_%d", i), Form("fZPhi Layer %d;z (cm);#phi (rad.)", i), 200, -15, 15, 200, 0, TMath::Pi() * 2);
69b09e3b 291 fOutput->Add(fZPhi[i]);
292 }
12bb57f1 293
294 fModuleMap = new TH1F("fModuleMap", "fModuleMap;module number;cluster count", 240, -0.5, 239.5);
295 fOutput->Add(fModuleMap);
69b09e3b 296 }
297
a7f69e56 298 if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
69b09e3b 299 {
300 fRawPt = new TH1F("fRawPt", "raw pt;p_{T};Count", 2000, 0, 100);
301 fOutput->Add(fRawPt);
302 }
ea441adf 303
1d107532 304 fVertex = new TH3F("vertex_check", "vertex_check;x;y;z", 100, -1, 1, 100, -1, 1, 100, -30, 30);
51f6de65 305 fOutput->Add(fVertex);
1d107532 306
307 fVertexVsMult = new TH3F("fVertexVsMult", "fVertexVsMult;x;y;multiplicity", 100, -1, 1, 100, -1, 1, 100, -0.5, 99.5);
308 fOutput->Add(fVertexVsMult);
51f6de65 309
0f67a57c 310 if (fReadMC)
311 {
770a1f1d 312 fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta", fAnalysisMode);
0f67a57c 313 fOutput->Add(fdNdEtaAnalysis);
314
69b09e3b 315 fdNdEtaAnalysisND = new dNdEtaAnalysis("dndetaND", "dndetaND", fAnalysisMode);
316 fOutput->Add(fdNdEtaAnalysisND);
317
0fc41645 318 fdNdEtaAnalysisNSD = new dNdEtaAnalysis("dndetaNSD", "dndetaNSD", fAnalysisMode);
319 fOutput->Add(fdNdEtaAnalysisNSD);
320
81be4ee8 321 fdNdEtaAnalysisOnePart = new dNdEtaAnalysis("dndetaOnePart", "dndetaOnePart", fAnalysisMode);
322 fOutput->Add(fdNdEtaAnalysisOnePart);
323
770a1f1d 324 fdNdEtaAnalysisTr = new dNdEtaAnalysis("dndetaTr", "dndetaTr", fAnalysisMode);
0f67a57c 325 fOutput->Add(fdNdEtaAnalysisTr);
326
770a1f1d 327 fdNdEtaAnalysisTrVtx = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx", fAnalysisMode);
0f67a57c 328 fOutput->Add(fdNdEtaAnalysisTrVtx);
329
770a1f1d 330 fdNdEtaAnalysisTracks = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks", fAnalysisMode);
0f67a57c 331 fOutput->Add(fdNdEtaAnalysisTracks);
332
0f67a57c 333 fPartPt = new TH1F("dndeta_check_pt", "dndeta_check_pt", 1000, 0, 10);
334 fPartPt->Sumw2();
335 fOutput->Add(fPartPt);
336 }
1c15d51a 337
338 if (fEsdTrackCuts)
339 {
340 fEsdTrackCuts->SetName("fEsdTrackCuts");
341 fOutput->Add(fEsdTrackCuts);
342 }
68fa248f 343
344 fEta = new TH1D("fEta", "Eta;#eta;count", 80, -4, 4);
345 fOutput->Add(fEta);
346
347 fEtaMC = new TH1D("fEtaMC", "Eta, MC;#eta;count", 80, -4, 4);
348 fOutput->Add(fEtaMC);
349
350 fHistEvents = new TH1D("fHistEvents", "N. of Events;accepted;count", 2, 0, 2);
351 fOutput->Add(fHistEvents);
352
353 fHistEventsMC = new TH1D("fHistEventsMC", "N. of MC Events;accepted;count", 2, 0, 2);
354 fOutput->Add(fHistEventsMC);
355
356 fTrigEffNum = new TH1D("fTrigEffNum", "N. of triggered events", 100,0,100);
357 fOutput->Add(fTrigEffNum);
358 fTrigEffDen = new TH1D("fTrigEffDen", "N. of true events", 100,0,100);
359 fOutput->Add(fTrigEffDen);
360 fVtxEffNum = new TH1D("fVtxEffNum", "N. of events with vtx", 100,0,100);
361 fOutput->Add(fVtxEffNum);
362 fVtxEffDen = new TH1D("fVtxEffDen", "N. of true events", 100,0,100);
363 fOutput->Add(fVtxEffDen);
364 fVtxTrigEffNum = new TH1D("fVtxTrigEffNum", "N. of triggered events with vtx", 100,0,100);
365 fOutput->Add(fVtxTrigEffNum);
366 fVtxTrigEffDen = new TH1D("fVtxTrigEffDen", "N. of triggered true events", 100,0,100);
367 fOutput->Add(fVtxTrigEffDen);
12bb57f1 368
12bb57f1 369 //AliLog::SetClassDebugLevel("AliPhysicsSelection", AliLog::kDebug);
0f67a57c 370}
371
76532b17 372Bool_t AlidNdEtaTask::IsEventInBinZero()
373{
374 // checks if the event goes to the 0 bin
375
68fa248f 376 Bool_t isZeroBin = kTRUE;
76532b17 377 fESD = (AliESDEvent*) fInputEvent;
378
379 AliInputEventHandler* inputHandler = static_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
380 if (!inputHandler)
381 {
382 Printf("ERROR: Could not receive input handler");
383 return kFALSE;
384 }
385
386 static AliTriggerAnalysis* triggerAnalysis = 0;
387 if (!triggerAnalysis)
388 {
389 AliPhysicsSelection* physicsSelection = dynamic_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
390 if (physicsSelection)
391 triggerAnalysis = physicsSelection->GetTriggerAnalysis();
392 }
393
394 if (!triggerAnalysis)
395 {
396 Printf("ERROR: Could not receive trigger analysis object");
397 return kFALSE;
398 }
399
400 if (!triggerAnalysis->IsTriggerFired(fESD, fTrigger))
401 return kFALSE;
402
68fa248f 403 // 0 bin check - from Michele
76532b17 404
68fa248f 405 const AliMultiplicity* mult = fESD->GetMultiplicity();
406 if (!mult){
407 Printf("AlidNdEtaTask::IsBinZero: Can't get multiplicity object from ESDs");
408 return kFALSE;
409 }
410 Int_t ntracklet = mult->GetNumberOfTracklets();
411 const AliESDVertex * vtxESD = fESD->GetPrimaryVertexSPD();
412 if(vtxESD) {
413 // If there is a vertex from vertexer z with delta phi > 0.02 we
414 // don't consider it rec (we keep the event in bin0). If quality
415 // is good eneough we check the number of tracklets
416 // if the vertex is more than 15 cm away, this is autamatically bin0
417 if( TMath::Abs(vtxESD->GetZ()) <= 15 ) {
418 if (vtxESD->IsFromVertexerZ()) {
419 if (vtxESD->GetDispersion()<=0.02 ) {
420 if(ntracklet>0) isZeroBin = kFALSE;
421 }
422 } else if(ntracklet>0) isZeroBin = kFALSE; // if the event is not from Vz we chek the n of tracklets
423 }
424 }
425 return isZeroBin;
76532b17 426}
427
428void AlidNdEtaTask::UserExec(Option_t*)
0f67a57c 429{
430 // process the event
431
76532b17 432 fESD = (AliESDEvent*) fInputEvent;
433
1c15d51a 434 // these variables are also used in the MC section below; however, if ESD is off they stay with their default values
435 Bool_t eventTriggered = kFALSE;
436 const AliESDVertex* vtxESD = 0;
437
51f6de65 438 // post the data already here
76532b17 439 PostData(1, fOutput);
51f6de65 440
1c15d51a 441 // ESD analysis
442 if (fESD)
0f67a57c 443 {
69b09e3b 444 AliESDHeader* esdHeader = fESD->GetHeader();
445 if (!esdHeader)
446 {
447 Printf("ERROR: esdHeader could not be retrieved");
448 return;
449 }
ff8c4f30 450
81be4ee8 451 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
452 if (!inputHandler)
1d107532 453 {
81be4ee8 454 Printf("ERROR: Could not receive input handler");
455 return;
12bb57f1 456 }
457
1103419a 458 // TODO use flags here!
81be4ee8 459 eventTriggered = inputHandler->IsEventSelected();
460
76532b17 461 static AliTriggerAnalysis* triggerAnalysis = 0;
462 if (!triggerAnalysis)
81be4ee8 463 {
464 AliPhysicsSelection* physicsSelection = dynamic_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
465 if (physicsSelection)
76532b17 466 triggerAnalysis = physicsSelection->GetTriggerAnalysis();
81be4ee8 467 }
468
469 if (eventTriggered)
76532b17 470 eventTriggered = triggerAnalysis->IsTriggerFired(fESD, fTrigger);
12bb57f1 471
76532b17 472 AliTriggerAnalysis::V0Decision v0A = triggerAnalysis->V0Trigger(fESD, AliTriggerAnalysis::kASide, kFALSE);
473 AliTriggerAnalysis::V0Decision v0C = triggerAnalysis->V0Trigger(fESD, AliTriggerAnalysis::kCSide, kFALSE);
474 Bool_t fastORHW = (triggerAnalysis->SPDFiredChips(fESD, 1) > 0);
12bb57f1 475
476 Int_t vZero = 0;
477 if (v0A != AliTriggerAnalysis::kV0Invalid && v0C != AliTriggerAnalysis::kV0Invalid)
478 {
479 if (v0A == AliTriggerAnalysis::kV0Empty && v0C == AliTriggerAnalysis::kV0Empty) vZero = 1;
480 if (v0A == AliTriggerAnalysis::kV0BB && v0C == AliTriggerAnalysis::kV0Empty) vZero = 2;
481 if (v0A == AliTriggerAnalysis::kV0Empty && v0C == AliTriggerAnalysis::kV0BB) vZero = 3;
482 if (v0A == AliTriggerAnalysis::kV0BB && v0C == AliTriggerAnalysis::kV0BB) vZero = 4;
483 if (v0A == AliTriggerAnalysis::kV0BG && v0C == AliTriggerAnalysis::kV0Empty) vZero = 5;
484 if (v0A == AliTriggerAnalysis::kV0Empty && v0C == AliTriggerAnalysis::kV0BG) vZero = 6;
485 if (v0A == AliTriggerAnalysis::kV0BG && v0C == AliTriggerAnalysis::kV0BG) vZero = 7;
486 if (v0A == AliTriggerAnalysis::kV0BB && v0C == AliTriggerAnalysis::kV0BG) vZero = 8;
487 if (v0A == AliTriggerAnalysis::kV0BG && v0C == AliTriggerAnalysis::kV0BB) vZero = 9;
488 }
81be4ee8 489
490 if (vZero == 0)
491 Printf("VZERO: %d %d %d %d", vZero, eventTriggered, v0A, v0C);
1d107532 492
12bb57f1 493 Bool_t filled = kFALSE;
494
12bb57f1 495 if (!eventTriggered)
496 {
497 fStats2->Fill(0.0, vZero);
76532b17 498 fStats2->Fill(0.0, (fastORHW) ? 10 : 11);
12bb57f1 499 filled = kTRUE;
500 }
501
12bb57f1 502 if (eventTriggered)
12bb57f1 503 fStats->Fill(3);
12bb57f1 504
81be4ee8 505 /*
506 // ITS cluster tree
507 AliESDInputHandlerRP* handlerRP = dynamic_cast<AliESDInputHandlerRP*> (AliAnalysisManager::GetAnalys
508isManager()->GetInputEventHandler());
509 if (handlerRP)
12bb57f1 510 {
81be4ee8 511 TTree* itsClusterTree = handlerRP->GetTreeR("ITS");
512 if (!itsClusterTree)
513 return;
514
515 TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
516 TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
517
518 itsClusterBranch->SetAddress(&itsClusters);
519
520 Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
521
522 Int_t totalClusters = 0;
523
524 // loop over the its subdetectors
525 for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
526
527 if (!itsClusterTree->GetEvent(iIts))
528 continue;
529
530 Int_t nClusters = itsClusters->GetEntriesFast();
531 totalClusters += nClusters;
532
533 #ifdef FULLALIROOT
534 if (fAnalysisMode & AliPWG0Helper::kSPD)
12bb57f1 535 {
81be4ee8 536 // loop over clusters
537 while (nClusters--) {
538 AliITSRecPoint* cluster = (AliITSRecPoint*) itsClusters->UncheckedAt(nClusters);
539
540 Int_t layer = cluster->GetLayer();
541
542 if (layer > 1)
543 continue;
544
545 Float_t xyz[3] = {0., 0., 0.};
546 cluster->GetGlobalXYZ(xyz);
547
548 Float_t phi = TMath::Pi() + TMath::ATan2(-xyz[1], -xyz[0]);
549 Float_t z = xyz[2];
550
551 fZPhi[layer]->Fill(z, phi);
552 fModuleMap->Fill(layer * 80 + cluster->GetDetectorIndex());
12bb57f1 553 }
1d107532 554 }
81be4ee8 555 #endif
1d107532 556 }
12bb57f1 557 }
81be4ee8 558 */
559
12bb57f1 560 vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
561 if (!vtxESD)
562 {
563 if (!filled)
564 {
81be4ee8 565 fStats2->Fill(1, vZero);
76532b17 566 fStats2->Fill(1, (fastORHW) ? 10 : 11);
12bb57f1 567 filled = kTRUE;
568 }
569 }
570 else
571 {
81be4ee8 572 if (!AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
573 {
574 if (!filled)
575 {
576 fStats2->Fill(2, vZero);
76532b17 577 fStats2->Fill(2, (fastORHW) ? 10 : 11);
81be4ee8 578 filled = kTRUE;
579 }
580 }
581
12bb57f1 582 Double_t vtx[3];
583 vtxESD->GetXYZ(vtx);
1d107532 584
81be4ee8 585 // try to compare spd vertex and vertexer from tracks
586 // remove vertices outside +- 15 cm
587 if (TMath::Abs(vtx[2]) > 15)
12bb57f1 588 {
589 if (!filled)
590 {
591 fStats2->Fill(3, vZero);
76532b17 592 fStats2->Fill(3, (fastORHW) ? 10 : 11);
12bb57f1 593 filled = kTRUE;
594 }
595 }
81be4ee8 596
597 if (TMath::Abs(vtx[2]) > 10)
598 {
599 if (!filled)
600 {
601 fStats2->Fill(4, vZero);
76532b17 602 fStats2->Fill(4, (fastORHW) ? 10 : 11);
81be4ee8 603 filled = kTRUE;
604 }
605 }
1d107532 606
12bb57f1 607 const AliMultiplicity* mult = fESD->GetMultiplicity();
608 if (!mult)
1103419a 609 {
68fa248f 610 Printf("Returning, no Multiplicity found");
12bb57f1 611 return;
1103419a 612 }
1d107532 613
12bb57f1 614 if (mult->GetNumberOfTracklets() == 0)
1d107532 615 {
12bb57f1 616 if (!filled)
1d107532 617 {
81be4ee8 618 fStats2->Fill(5, vZero);
76532b17 619 fStats2->Fill(5, (fastORHW) ? 10 : 11);
12bb57f1 620 filled = kTRUE;
1d107532 621 }
622 }
12bb57f1 623 }
624
12bb57f1 625 if (!filled)
81be4ee8 626 {
12bb57f1 627 fStats2->Fill(6, vZero);
76532b17 628 fStats2->Fill(6, (fastORHW) ? 10 : 11);
81be4ee8 629 //Printf("File: %s, IEV: %d, TRG: ---, Orbit: 0x%x, Period: %d, BC: %d", ((TTree*) GetInputData(0))->GetCurrentFile()->GetName(), fESD->GetEventNumberInFile(), fESD->GetOrbitNumber(),fESD->GetPeriodNumber(),fESD->GetBunchCrossNumber());
630 }
1d107532 631
a7f69e56 632 if (eventTriggered)
633 fStats->Fill(3);
12bb57f1 634
635 fStats->Fill(5);
636
1c15d51a 637 // get the ESD vertex
638 vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
ff8c4f30 639
1c15d51a 640 Double_t vtx[3];
3d7758c1 641
1c15d51a 642 // fill z vertex resolution
643 if (vtxESD)
51f6de65 644 {
81be4ee8 645 if (strcmp(vtxESD->GetTitle(), "vertexer: Z") == 0)
646 fVertexResolution->Fill(vtxESD->GetZRes(), vtxESD->GetDispersion());
647 else
648 fVertexResolution->Fill(vtxESD->GetZRes(), 0);
649
1d107532 650 //if (strcmp(vtxESD->GetTitle(), "vertexer: 3D") == 0)
651 {
652 fVertex->Fill(vtxESD->GetXv(), vtxESD->GetYv(), vtxESD->GetZv());
653 }
654
51f6de65 655 if (AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
656 {
657 vtxESD->GetXYZ(vtx);
69b09e3b 658
659 // vertex stats
660 if (strcmp(vtxESD->GetTitle(), "vertexer: 3D") == 0)
69b09e3b 661 fStats->Fill(1);
81be4ee8 662
663 if (strcmp(vtxESD->GetTitle(), "vertexer: Z") == 0)
69b09e3b 664 fStats->Fill(2);
81be4ee8 665
666 // remove vertices outside +- 15 cm
667 if (TMath::Abs(vtx[2]) > 15)
668 vtxESD = 0;
51f6de65 669 }
670 else
671 vtxESD = 0;
672 }
3d7758c1 673
1c15d51a 674 // needed for syst. studies
675 AliStack* stack = 0;
676 TArrayF vtxMC(3);
677
c17301f3 678 if (fUseMCVertex || fUseMCKine || fOnlyPrimaries || fReadMC) {
1c15d51a 679 if (!fReadMC) {
c17301f3 680 Printf("ERROR: fUseMCVertex or fUseMCKine or fOnlyPrimaries set without fReadMC set!");
1c15d51a 681 return;
682 }
683
684 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
685 if (!eventHandler) {
686 Printf("ERROR: Could not retrieve MC event handler");
687 return;
688 }
689
690 AliMCEvent* mcEvent = eventHandler->MCEvent();
691 if (!mcEvent) {
692 Printf("ERROR: Could not retrieve MC event");
693 return;
694 }
695
696 AliHeader* header = mcEvent->Header();
697 if (!header)
698 {
699 AliDebug(AliLog::kError, "Header not available");
700 return;
701 }
3d7758c1 702
3d7758c1 703 // get the MC vertex
704 AliGenEventHeader* genHeader = header->GenEventHeader();
1c15d51a 705 if (!genHeader)
706 {
707 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
708 return;
709 }
3d7758c1 710 genHeader->PrimaryVertex(vtxMC);
711
1c15d51a 712 if (fUseMCVertex)
1c15d51a 713 vtx[2] = vtxMC[2];
3d7758c1 714
c2b08c81 715 stack = mcEvent->Stack();
716 if (!stack)
3d7758c1 717 {
c2b08c81 718 AliDebug(AliLog::kError, "Stack not available");
719 return;
3d7758c1 720 }
721 }
81be4ee8 722
1c15d51a 723 // create list of (label, eta, pt) tuples
724 Int_t inputCount = 0;
725 Int_t* labelArr = 0;
726 Float_t* etaArr = 0;
69b09e3b 727 Float_t* thirdDimArr = 0;
a7f69e56 728 if (fAnalysisMode & AliPWG0Helper::kSPD)
0f67a57c 729 {
81be4ee8 730 if (vtxESD)
1c15d51a 731 {
81be4ee8 732 // get tracklets
733 const AliMultiplicity* mult = fESD->GetMultiplicity();
734 if (!mult)
1d107532 735 {
81be4ee8 736 AliDebug(AliLog::kError, "AliMultiplicity not available");
737 return;
738 }
739
740 Int_t arrayLength = mult->GetNumberOfTracklets();
741 if (fAnalysisMode & AliPWG0Helper::kSPDOnlyL0)
742 arrayLength += mult->GetNumberOfSingleClusters();
743
744 labelArr = new Int_t[arrayLength];
745 etaArr = new Float_t[arrayLength];
746 thirdDimArr = new Float_t[arrayLength];
747
748 // get multiplicity from SPD tracklets
749 for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
750 {
751 //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
752
753 if (fOnlyPrimaries)
754 if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
755 continue;
756
757 Float_t deltaPhi = mult->GetDeltaPhi(i);
758 // prevent values to be shifted by 2 Pi()
759 if (deltaPhi < -TMath::Pi())
760 deltaPhi += TMath::Pi() * 2;
761 if (deltaPhi > TMath::Pi())
762 deltaPhi -= TMath::Pi() * 2;
763
764 if (TMath::Abs(deltaPhi) > 1)
765 printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
766
767 Int_t label = mult->GetLabel(i, 0);
768 Float_t eta = mult->GetEta(i);
769
770 // control histograms
771 Float_t phi = mult->GetPhi(i);
772 if (phi < 0)
773 phi += TMath::Pi() * 2;
774
775 // TEST exclude probably inefficient phi region
776 //if (phi > 5.70 || phi < 0.06)
777 // continue;
778
779 fPhi->Fill(phi);
780
781 if (vtxESD && TMath::Abs(vtx[2]) < 10)
782 {
783 fEtaPhi->Fill(eta, phi);
784 fDeltaPhi->Fill(deltaPhi);
785 fDeltaTheta->Fill(mult->GetDeltaTheta(i));
786 }
787
788 if (fDeltaPhiCut > 0 && (TMath::Abs(deltaPhi) > fDeltaPhiCut || TMath::Abs(mult->GetDeltaTheta(i)) > fDeltaPhiCut / 0.08 * 0.025))
789 continue;
790
791 if (fUseMCKine)
792 {
793 if (label > 0)
794 {
795 TParticle* particle = stack->Particle(label);
796 eta = particle->Eta();
797 phi = particle->Phi();
798 }
799 else
800 Printf("WARNING: fUseMCKine set without fOnlyPrimaries and no label found");
801 }
802
803 if (fSymmetrize)
804 eta = TMath::Abs(eta);
805
806 etaArr[inputCount] = eta;
807 labelArr[inputCount] = label;
808 thirdDimArr[inputCount] = phi;
809 ++inputCount;
1d107532 810 }
811
81be4ee8 812 if (fAnalysisMode & AliPWG0Helper::kSPDOnlyL0)
51f6de65 813 {
81be4ee8 814 // get additional clusters from L0
815 for (Int_t i=0; i<mult->GetNumberOfSingleClusters(); ++i)
51f6de65 816 {
81be4ee8 817 etaArr[inputCount] = -TMath::Log(TMath::Tan(mult->GetThetaSingle(i)/2.));
818 labelArr[inputCount] = -1;
819 thirdDimArr[inputCount] = mult->GetPhiSingle(i);
820
821 ++inputCount;
51f6de65 822 }
51f6de65 823 }
69b09e3b 824
81be4ee8 825 if (!fFillPhi)
826 {
827 // fill multiplicity in third axis
828 for (Int_t i=0; i<inputCount; ++i)
829 thirdDimArr[i] = inputCount;
830 }
831
832 Int_t firedChips = mult->GetNumberOfFiredChips(0) + mult->GetNumberOfFiredChips(1);
833 fFiredChips->Fill(firedChips, inputCount);
834 Printf("Accepted %d tracklets (%d fired chips)", inputCount, firedChips);
835
836 fTrackletsVsUnassigned->Fill(inputCount, mult->GetNumberOfSingleClusters());
69b09e3b 837 }
0f67a57c 838 }
a7f69e56 839 else if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
0f67a57c 840 {
1c15d51a 841 if (!fEsdTrackCuts)
842 {
843 AliDebug(AliLog::kError, "fESDTrackCuts not available");
844 return;
845 }
0f67a57c 846
81be4ee8 847 Bool_t foundInEta10 = kFALSE;
848
69b09e3b 849 if (vtxESD)
0f67a57c 850 {
69b09e3b 851 // get multiplicity from ESD tracks
a7f69e56 852 TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, fAnalysisMode & AliPWG0Helper::kTPC);
69b09e3b 853 Int_t nGoodTracks = list->GetEntries();
a7f69e56 854 Printf("Accepted %d tracks out of %d total ESD tracks", nGoodTracks, fESD->GetNumberOfTracks());
69b09e3b 855
856 labelArr = new Int_t[nGoodTracks];
857 etaArr = new Float_t[nGoodTracks];
858 thirdDimArr = new Float_t[nGoodTracks];
859
860 // loop over esd tracks
861 for (Int_t i=0; i<nGoodTracks; i++)
1c15d51a 862 {
69b09e3b 863 AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
864 if (!esdTrack)
c17301f3 865 {
69b09e3b 866 AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
867 continue;
c17301f3 868 }
69b09e3b 869
870 Float_t phi = esdTrack->Phi();
871 if (phi < 0)
872 phi += TMath::Pi() * 2;
873
874 Float_t eta = esdTrack->Eta();
875 Int_t label = TMath::Abs(esdTrack->GetLabel());
876 Float_t pT = esdTrack->Pt();
a7f69e56 877
878 // force pT to fixed value without B field
879 if (!(fAnalysisMode & AliPWG0Helper::kFieldOn))
880 pT = 1;
69b09e3b 881
882 fPhi->Fill(phi);
883 fEtaPhi->Fill(eta, phi);
884 if (eventTriggered && vtxESD)
885 fRawPt->Fill(pT);
886
68fa248f 887 if (esdTrack->Pt() < fPtMin) // even if the pt cut is already in the esdTrackCuts....
81be4ee8 888 continue;
889
a7f69e56 890 if (fOnlyPrimaries)
891 {
68fa248f 892 if (label == 0)
893 continue;
894
895 if (stack->IsPhysicalPrimary(label) == kFALSE)
896 continue;
a7f69e56 897 }
68fa248f 898
899 // 2 types of INEL>0 trigger - choose one
900
901 // 1. HL
902 // INEL>0 trigger
903 // if (TMath::Abs(esdTrack->Eta()) < 1 && esdTrack->Pt() > 0.15)
904 // foundInEta10 = kTRUE;
905
906 // 2. MB working group
907 if (TMath::Abs(esdTrack->Eta()) < 0.8 && esdTrack->Pt() > fPtMin){ // this should be in the trigger selection as well, so one particle should always be found
908 foundInEta10 = kTRUE;
909 }
69b09e3b 910
911 if (fUseMCKine)
912 {
68fa248f 913 if (label > 0)
69b09e3b 914 {
915 TParticle* particle = stack->Particle(label);
916 eta = particle->Eta();
917 pT = particle->Pt();
68fa248f 918 // check when using INEL>0, MB Working Group definition
919 if (TMath::Abs(eta) >=0.8 || pT <= fPtMin){
920 AliDebug(2,Form("WARNING ************* USING KINE: eta = %f, pt = %f",eta,pT));
921 }
69b09e3b 922 }
923 else
924 Printf("WARNING: fUseMCKine set without fOnlyPrimaries and no label found");
68fa248f 925 }
69b09e3b 926
1d107532 927 if (fSymmetrize)
928 eta = TMath::Abs(eta);
69b09e3b 929 etaArr[inputCount] = eta;
930 labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
931 thirdDimArr[inputCount] = pT;
932 ++inputCount;
c17301f3 933 }
69b09e3b 934
81be4ee8 935 //if (inputCount > 30)
936 // Printf("Event with %d accepted TPC tracks. Period number: %d Orbit number: %x Bunch crossing number: %d", inputCount, fESD->GetPeriodNumber(), fESD->GetOrbitNumber(), fESD->GetBunchCrossNumber());
12bb57f1 937
69b09e3b 938 // TODO restrict inputCount used as measure for the multiplicity to |eta| < 1
939
940 delete list;
1c15d51a 941 }
81be4ee8 942
68fa248f 943 if (!foundInEta10){
81be4ee8 944 eventTriggered = kFALSE;
68fa248f 945 fHistEvents->Fill(0.1);
946 AliDebug(3,"Event rejected");
947 }
948 else{
949 if (eventTriggered){
950 fHistEvents->Fill(1.1);
951 AliDebug(3,Form("Event Accepted, with inputcount = %d",inputCount));
952 }
953 else{
954 AliDebug(3,"Event has foundInEta10 but was not triggered");
955 }
956 }
1c15d51a 957 }
958 else
959 return;
7307d52c 960
1c15d51a 961 // Processing of ESD information (always)
962 if (eventTriggered)
dd367a14 963 {
964 // control hist
1c15d51a 965 fMult->Fill(inputCount);
966 fdNdEtaAnalysisESD->FillTriggeredEvent(inputCount);
0f67a57c 967
1c15d51a 968 if (vtxESD)
dd367a14 969 {
1c15d51a 970 // control hist
1d107532 971
972 if (strcmp(vtxESD->GetTitle(), "vertexer: 3D") == 0)
973 fVertexVsMult->Fill(vtxESD->GetXv(), vtxESD->GetYv(), inputCount);
974
81be4ee8 975 Int_t part05 = 0;
976 Int_t part10 = 0;
1c15d51a 977 for (Int_t i=0; i<inputCount; ++i)
dd367a14 978 {
1c15d51a 979 Float_t eta = etaArr[i];
69b09e3b 980 Float_t thirdDim = thirdDimArr[i];
0f67a57c 981
69b09e3b 982 fdNdEtaAnalysisESD->FillTrack(vtx[2], eta, thirdDim);
1c15d51a 983
69b09e3b 984 if (TMath::Abs(vtx[2]) < 10)
1c15d51a 985 {
986 fPartEta[0]->Fill(eta);
987
988 if (vtx[2] < 0)
989 fPartEta[1]->Fill(eta);
990 else
991 fPartEta[2]->Fill(eta);
992 }
81be4ee8 993
994 if (TMath::Abs(eta) < 0.5)
995 part05++;
68fa248f 996 if (TMath::Abs(eta) < 1.0) // in the INEL>0, MB WG definition, this is in any case equivalent to <0.8, due to the EsdTrackCuts
81be4ee8 997 part10++;
dd367a14 998 }
81be4ee8 999
1000 Int_t multAxis = inputCount;
1001 if (fMultAxisEta1)
1002 multAxis = part10;
1003
1004 fMultVtx->Fill(multAxis);
1005 //if (TMath::Abs(vtx[2]) < 10)
1006 // fMultVtx->Fill(part05);
0f67a57c 1007
1c15d51a 1008 // for event count per vertex
81be4ee8 1009 fdNdEtaAnalysisESD->FillEvent(vtx[2], multAxis);
0f67a57c 1010
1c15d51a 1011 // control hist
81be4ee8 1012 if (multAxis > 0)
12bb57f1 1013 fEvents->Fill(vtx[2]);
1c15d51a 1014
1015 if (fReadMC)
1016 {
1017 // from tracks is only done for triggered and vertex reconstructed events
1018 for (Int_t i=0; i<inputCount; ++i)
1019 {
1020 Int_t label = labelArr[i];
1021
1022 if (label < 0)
1023 continue;
1024
1025 //Printf("Getting particle of track %d", label);
1026 TParticle* particle = stack->Particle(label);
1027
1028 if (!particle)
1029 {
1030 AliDebug(AliLog::kError, Form("ERROR: Could not retrieve particle %d.", label));
1031 continue;
1032 }
1033
69b09e3b 1034 Float_t thirdDim = -1;
a7f69e56 1035 if (fAnalysisMode & AliPWG0Helper::kSPD)
69b09e3b 1036 {
1037 if (fFillPhi)
1038 {
1039 thirdDim = particle->Phi();
1040 }
1041 else
81be4ee8 1042 thirdDim = multAxis;
69b09e3b 1043 }
1044 else
1045 thirdDim = particle->Pt();
1046
1d107532 1047 Float_t eta = particle->Eta();
1048 if (fSymmetrize)
1049 eta = TMath::Abs(eta);
1050 fdNdEtaAnalysisTracks->FillTrack(vtxMC[2], eta, thirdDim);
1c15d51a 1051 } // end of track loop
1052
1053 // for event count per vertex
81be4ee8 1054 fdNdEtaAnalysisTracks->FillEvent(vtxMC[2], multAxis);
1c15d51a 1055 }
1056 }
dd367a14 1057 }
69b09e3b 1058 if (etaArr)
1059 delete[] etaArr;
1060 if (labelArr)
1061 delete[] labelArr;
1062 if (thirdDimArr)
1063 delete[] thirdDimArr;
0f67a57c 1064 }
1065
1066 if (fReadMC) // Processing of MC information (optional)
1067 {
1068 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
1069 if (!eventHandler) {
1070 Printf("ERROR: Could not retrieve MC event handler");
1071 return;
1072 }
1073
1074 AliMCEvent* mcEvent = eventHandler->MCEvent();
1075 if (!mcEvent) {
1076 Printf("ERROR: Could not retrieve MC event");
1077 return;
1078 }
1079
1c15d51a 1080 AliStack* stack = mcEvent->Stack();
0f67a57c 1081 if (!stack)
1082 {
1083 AliDebug(AliLog::kError, "Stack not available");
1084 return;
1085 }
1086
1087 AliHeader* header = mcEvent->Header();
1088 if (!header)
1089 {
1090 AliDebug(AliLog::kError, "Header not available");
1091 return;
1092 }
1093
1094 // get the MC vertex
1095 AliGenEventHeader* genHeader = header->GenEventHeader();
567160d6 1096 if (!genHeader)
1097 {
1098 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
1099 return;
1100 }
1101
0f67a57c 1102 TArrayF vtxMC(3);
1103 genHeader->PrimaryVertex(vtxMC);
81be4ee8 1104
6b62a9c7 1105 // get process type
81be4ee8 1106 Int_t processType = AliPWG0Helper::GetEventProcessType(fESD, header, stack, fDiffTreatment);
6b62a9c7 1107 AliDebug(AliLog::kDebug+1, Form("Found process type %d", processType));
0fc41645 1108
81be4ee8 1109 if (processType == AliPWG0Helper::kInvalidProcess)
6b62a9c7 1110 AliDebug(AliLog::kError, Form("Unknown process type %d.", processType));
0fc41645 1111
0f67a57c 1112 // loop over mc particles
1113 Int_t nPrim = stack->GetNprimary();
1114
1115 Int_t nAcceptedParticles = 0;
81be4ee8 1116 Bool_t oneParticleEvent = kFALSE;
0f67a57c 1117
68fa248f 1118 Int_t nMCparticlesinRange = 0; // number of true particles in my range of interest
1119
51f6de65 1120 // count particles first, then fill
0f67a57c 1121 for (Int_t iMc = 0; iMc < nPrim; ++iMc)
1122 {
1c15d51a 1123 if (!stack->IsPhysicalPrimary(iMc))
1124 continue;
1125
0f67a57c 1126 //Printf("Getting particle %d", iMc);
1127 TParticle* particle = stack->Particle(iMc);
1128
1129 if (!particle)
1130 continue;
1131
1c15d51a 1132 if (particle->GetPDG()->Charge() == 0)
0f67a57c 1133 continue;
1134
1c15d51a 1135 //AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID()));
0f67a57c 1136 Float_t eta = particle->Eta();
0f67a57c 1137
68fa248f 1138 AliDebug(2,Form("particle %d: eta = %f, pt = %f",iMc,particle->Eta(),particle->Pt()));
1139
1140 // INEL>0: choose one
1141 // 1.HL
1142 // if (TMath::Abs(eta) < 1.0){
1143 // oneParticleEvent = kTRUE;
1144 // nMCparticlesinRange++;
1145 //}
1146 // 2.MB Working Group definition
1147 if (TMath::Abs(eta) < 0.8 && particle->Pt() > fPtMin){
1148 oneParticleEvent = kTRUE;
1149 nMCparticlesinRange++;
1150 }
81be4ee8 1151
1152 // make a rough eta cut (so that nAcceptedParticles is not too far off the true measured value (NB: these histograms are only gathered for comparison))
ea441adf 1153 if (TMath::Abs(eta) < 1.5) // && pt > 0.3)
2b26296f 1154 nAcceptedParticles++;
51f6de65 1155 }
1156
68fa248f 1157 if (TMath::Abs(vtxMC[2]) < 10.){
1158 // if MC vertex is inside 10
1159 if (vtxESD){
1160 fVtxEffNum->Fill(nMCparticlesinRange);
1161 }
1162 fVtxEffDen->Fill(nMCparticlesinRange);
1163 if (eventTriggered){
1164 if (vtxESD){
1165 fVtxTrigEffNum->Fill(nMCparticlesinRange);
1166 }
1167 fVtxTrigEffDen->Fill(nMCparticlesinRange);
1168 fTrigEffNum->Fill(nMCparticlesinRange);
1169 }
1170 fTrigEffDen->Fill(nMCparticlesinRange);
1171 }
1172
1173 if (oneParticleEvent){
1174 fHistEventsMC->Fill(1.1);
1175 }
1176 else{
1177 fHistEventsMC->Fill(0.1);
1178 }
1179
51f6de65 1180 for (Int_t iMc = 0; iMc < nPrim; ++iMc)
1181 {
1182 if (!stack->IsPhysicalPrimary(iMc))
1183 continue;
2b26296f 1184
51f6de65 1185 //Printf("Getting particle %d", iMc);
1186 TParticle* particle = stack->Particle(iMc);
1187
1188 if (!particle)
1189 continue;
1190
1191 if (particle->GetPDG()->Charge() == 0)
1192 continue;
1193
1194 Float_t eta = particle->Eta();
1d107532 1195 if (fSymmetrize)
1196 eta = TMath::Abs(eta);
1197
69b09e3b 1198 Float_t thirdDim = -1;
1199
a7f69e56 1200 if (fAnalysisMode & AliPWG0Helper::kSPD)
69b09e3b 1201 {
1202 if (fFillPhi)
1203 {
1204 thirdDim = particle->Phi();
1205 }
1206 else
1207 thirdDim = nAcceptedParticles;
1208 }
1209 else
1210 thirdDim = particle->Pt();
51f6de65 1211
1212 fdNdEtaAnalysis->FillTrack(vtxMC[2], eta, thirdDim);
0f67a57c 1213
69b09e3b 1214 if (processType != AliPWG0Helper::kSD)
51f6de65 1215 fdNdEtaAnalysisNSD->FillTrack(vtxMC[2], eta, thirdDim);
0fc41645 1216
69b09e3b 1217 if (processType == AliPWG0Helper::kND)
1218 fdNdEtaAnalysisND->FillTrack(vtxMC[2], eta, thirdDim);
81be4ee8 1219
1220 if (oneParticleEvent)
68fa248f 1221 AliDebug(3,Form("filling dNdEtaAnalysis object:: vtx = %f, eta = %f, pt = %f",vtxMC[2], eta, thirdDim));
81be4ee8 1222 fdNdEtaAnalysisOnePart->FillTrack(vtxMC[2], eta, thirdDim);
69b09e3b 1223
0f67a57c 1224 if (eventTriggered)
1225 {
51f6de65 1226 fdNdEtaAnalysisTr->FillTrack(vtxMC[2], eta, thirdDim);
770a1f1d 1227 if (vtxESD)
51f6de65 1228 fdNdEtaAnalysisTrVtx->FillTrack(vtxMC[2], eta, thirdDim);
0f67a57c 1229 }
1230
ff8c4f30 1231 if (TMath::Abs(eta) < 1.0 && particle->Pt() > 0 && particle->P() > 0)
1232 {
1d107532 1233 //Float_t value = 1. / TMath::TwoPi() / particle->Pt() * particle->Energy() / particle->P();
1234 Float_t value = 1;
ff8c4f30 1235 fPartPt->Fill(particle->Pt(), value);
68fa248f 1236 if (TMath::Abs(eta) < 0.8 && particle->Pt() > fPtMin){
1237 fEtaMC->Fill(eta);
1238 }
ff8c4f30 1239 }
0f67a57c 1240 }
1241
1242 fdNdEtaAnalysis->FillEvent(vtxMC[2], nAcceptedParticles);
ea441adf 1243 if (processType != AliPWG0Helper::kSD)
0fc41645 1244 fdNdEtaAnalysisNSD->FillEvent(vtxMC[2], nAcceptedParticles);
69b09e3b 1245 if (processType == AliPWG0Helper::kND)
1246 fdNdEtaAnalysisND->FillEvent(vtxMC[2], nAcceptedParticles);
81be4ee8 1247 if (oneParticleEvent)
1248 fdNdEtaAnalysisOnePart->FillEvent(vtxMC[2], nAcceptedParticles);
0fc41645 1249
0f67a57c 1250 if (eventTriggered)
1251 {
1252 fdNdEtaAnalysisTr->FillEvent(vtxMC[2], nAcceptedParticles);
770a1f1d 1253 if (vtxESD)
0f67a57c 1254 fdNdEtaAnalysisTrVtx->FillEvent(vtxMC[2], nAcceptedParticles);
1255 }
0f67a57c 1256 }
0f67a57c 1257}
1258
1259void AlidNdEtaTask::Terminate(Option_t *)
1260{
1261 // The Terminate() function is the last function to be called during
1262 // a query. It always runs on the client, it can be used to present
1263 // the results graphically or save the results to file.
1264
76532b17 1265 fOutput = dynamic_cast<TList*> (GetOutputData(1));
1c15d51a 1266 if (!fOutput)
0f67a57c 1267 Printf("ERROR: fOutput not available");
0f67a57c 1268
68fa248f 1269 TH1D* dNdEta = new TH1D("dNdEta","#eta;counts",80, -4, 4);
1270 TH1D* dNdEtaMC = new TH1D("dNdEtaMC","#eta,MC;counts",80, -4, 4);
1c15d51a 1271 if (fOutput)
1272 {
1273 fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("fdNdEtaAnalysisESD"));
1274 fMult = dynamic_cast<TH1F*> (fOutput->FindObject("fMult"));
1275 fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx"));
1276 for (Int_t i=0; i<3; ++i)
1277 fPartEta[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("dndeta_check_%d", i)));
1278 fEvents = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_vertex"));
81be4ee8 1279 fVertexResolution = dynamic_cast<TH2F*> (fOutput->FindObject("dndeta_vertex_resolution_z"));
ea441adf 1280
1d107532 1281 fVertex = dynamic_cast<TH3F*> (fOutput->FindObject("vertex_check"));
1282 fVertexVsMult = dynamic_cast<TH3F*> (fOutput->FindObject("fVertexVsMult"));
1c15d51a 1283 fPhi = dynamic_cast<TH1F*> (fOutput->FindObject("fPhi"));
69b09e3b 1284 fRawPt = dynamic_cast<TH1F*> (fOutput->FindObject("fRawPt"));
1c15d51a 1285 fEtaPhi = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaPhi"));
69b09e3b 1286 for (Int_t i=0; i<2; ++i)
1287 fZPhi[i] = dynamic_cast<TH2F*> (fOutput->FindObject(Form("fZPhi_%d", i)));
12bb57f1 1288 fModuleMap = dynamic_cast<TH1F*> (fOutput->FindObject("fModuleMap"));
1c15d51a 1289 fDeltaPhi = dynamic_cast<TH1F*> (fOutput->FindObject("fDeltaPhi"));
a7f69e56 1290 fDeltaTheta = dynamic_cast<TH1F*> (fOutput->FindObject("fDeltaTheta"));
69b09e3b 1291 fFiredChips = dynamic_cast<TH2F*> (fOutput->FindObject("fFiredChips"));
1d107532 1292 fTrackletsVsClusters = dynamic_cast<TH2F*> (fOutput->FindObject("fTrackletsVsClusters"));
1293 fTrackletsVsUnassigned = dynamic_cast<TH2F*> (fOutput->FindObject("fTrackletsVsUnassigned"));
69b09e3b 1294 fStats = dynamic_cast<TH1F*> (fOutput->FindObject("fStats"));
1d107532 1295 fStats2 = dynamic_cast<TH2F*> (fOutput->FindObject("fStats2"));
1c15d51a 1296
1297 fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCuts"));
68fa248f 1298
1299 fEta = dynamic_cast<TH1D*>(fOutput->FindObject("fEta"));
1300 fEtaMC = dynamic_cast<TH1D*>(fOutput->FindObject("fEtaMC"));
1301 fHistEvents = dynamic_cast<TH1D*>(fOutput->FindObject("fHistEvents"));
1302 fHistEventsMC = dynamic_cast<TH1D*>(fOutput->FindObject("fHistEventsMC"));
1303 fVtxEffDen = dynamic_cast<TH1D*>(fOutput->FindObject("fVtxEffDen"));
1304 fVtxEffNum = dynamic_cast<TH1D*>(fOutput->FindObject("fVtxEffNum"));
1305 fTrigEffDen = dynamic_cast<TH1D*>(fOutput->FindObject("fTrigEffDen"));
1306 fTrigEffNum = dynamic_cast<TH1D*>(fOutput->FindObject("fTrigEffNum"));
1307 fVtxTrigEffDen = dynamic_cast<TH1D*>(fOutput->FindObject("fVtxTrigEffDen"));
1308 fVtxTrigEffNum = dynamic_cast<TH1D*>(fOutput->FindObject("fVtxTrigEffNum"));
1309 Printf("Events selected = %f", fHistEvents->GetBinContent(fHistEvents->GetXaxis()->FindBin(1.1)));
1310 Printf("Events selected frm MC = %f", fHistEventsMC->GetBinContent(fHistEventsMC->GetXaxis()->FindBin(1.1)));
1311 Float_t nevents = fHistEvents->GetBinContent(fHistEvents->GetXaxis()->FindBin(1.1));
1312 Float_t neventsMC = fHistEventsMC->GetBinContent(fHistEventsMC->GetXaxis()->FindBin(1.1));
1313 for (Int_t ibin = 1; ibin <= fEta->GetNbinsX(); ibin++){
1314 Float_t eta = fEta->GetBinContent(ibin)/nevents/fEta->GetBinWidth(ibin);
1315 Float_t etaerr = fEta->GetBinError(ibin)/nevents/fEta->GetBinWidth(ibin);
1316 Float_t etaMC = fEtaMC->GetBinContent(ibin)/neventsMC/fEtaMC->GetBinWidth(ibin);
1317 Float_t etaerrMC = fEtaMC->GetBinError(ibin)/neventsMC/fEtaMC->GetBinWidth(ibin);
1318 dNdEta->SetBinContent(ibin,eta);
1319 dNdEta->SetBinError(ibin,etaerr);
1320 dNdEtaMC->SetBinContent(ibin,etaMC);
1321 dNdEtaMC->SetBinError(ibin,etaerrMC);
1322 }
1323 new TCanvas("eta", " eta ",50, 50, 550, 550) ;
1324 fEta->Draw();
1325 new TCanvas("etaMC", " etaMC ",50, 50, 550, 550) ;
1326 fEtaMC->Draw();
1327 new TCanvas("dNdEta", "#eta;dNdEta ",50, 50, 550, 550) ;
1328 dNdEta->Draw();
1329 new TCanvas("dNdEtaMC", "#eta,MC;dNdEta ",50, 50, 550, 550) ;
1330 dNdEtaMC->Draw();
1331 new TCanvas("Events", "Events;Events ",50, 50, 550, 550) ;
1332 fHistEvents->Draw();
1333 new TCanvas("Events, MC", "Events, MC;Events ",50, 50, 550, 550) ;
1334 fHistEventsMC->Draw();
1335
1336 TFile* outputFileCheck = new TFile("histogramsCheck.root", "RECREATE");
1337 fEta->Write();
1338 fEtaMC->Write();
1339 dNdEta->Write();
1340 dNdEtaMC->Write();
1341 outputFileCheck->Write();
1342 outputFileCheck->Close();
1343
1c15d51a 1344 }
0f67a57c 1345
1346 if (!fdNdEtaAnalysisESD)
1347 {
1348 AliDebug(AliLog::kError, "ERROR: fdNdEtaAnalysisESD not available");
0f67a57c 1349 }
51f6de65 1350 else
0f67a57c 1351 {
51f6de65 1352 if (fMult && fMultVtx)
1353 {
1354 new TCanvas;
1355 fMult->Draw();
1356 fMultVtx->SetLineColor(2);
1357 fMultVtx->Draw("SAME");
1358 }
0f67a57c 1359
69b09e3b 1360 if (fFiredChips)
1361 {
1362 new TCanvas;
1363 fFiredChips->Draw("COLZ");
1364 }
1365
51f6de65 1366 if (fPartEta[0])
1367 {
81be4ee8 1368 Int_t events1 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(-9.9), fEvents->GetXaxis()->FindBin(-0.001));
1369 Int_t events2 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(0.001), fEvents->GetXaxis()->FindBin(9.9));
0f67a57c 1370
51f6de65 1371 Printf("%d events with vertex used", events1 + events2);
81be4ee8 1372 Printf("%d total events with vertex, %d total triggered events, %d triggered events with 0 multiplicity", (Int_t) fMultVtx->Integral(), (Int_t) fMult->Integral(), (Int_t) fMult->GetBinContent(1));
54b096ef 1373
51f6de65 1374 if (events1 > 0 && events2 > 0)
1375 {
1376 fPartEta[0]->Scale(1.0 / (events1 + events2));
1377 fPartEta[1]->Scale(1.0 / events1);
1378 fPartEta[2]->Scale(1.0 / events2);
0f67a57c 1379
51f6de65 1380 for (Int_t i=0; i<3; ++i)
1381 fPartEta[i]->Scale(1.0 / fPartEta[i]->GetBinWidth(1));
0f67a57c 1382
51f6de65 1383 new TCanvas("control", "control", 500, 500);
1384 for (Int_t i=0; i<3; ++i)
1385 {
1386 fPartEta[i]->SetLineColor(i+1);
1387 fPartEta[i]->Draw((i==0) ? "" : "SAME");
1388 }
1389 }
0f67a57c 1390 }
51f6de65 1391
1392 if (fEvents)
1393 {
1394 new TCanvas("control3", "control3", 500, 500);
1395 fEvents->Draw();
54b096ef 1396 }
12bb57f1 1397
1398 if (fStats2)
1399 {
1400 new TCanvas;
1401 fStats2->Draw("TEXT");
1402 }
0f67a57c 1403
51f6de65 1404 TFile* fout = new TFile("analysis_esd_raw.root", "RECREATE");
0f67a57c 1405
51f6de65 1406 if (fdNdEtaAnalysisESD)
1407 fdNdEtaAnalysisESD->SaveHistograms();
0f67a57c 1408
51f6de65 1409 if (fEsdTrackCuts)
a7f69e56 1410 fEsdTrackCuts->SaveHistograms("esd_track_cuts");
0f67a57c 1411
51f6de65 1412 if (fMult)
1413 fMult->Write();
0f67a57c 1414
51f6de65 1415 if (fMultVtx)
1416 fMultVtx->Write();
0f67a57c 1417
51f6de65 1418 for (Int_t i=0; i<3; ++i)
1419 if (fPartEta[i])
1420 fPartEta[i]->Write();
0f67a57c 1421
51f6de65 1422 if (fEvents)
1423 fEvents->Write();
0f67a57c 1424
51f6de65 1425 if (fVertexResolution)
81be4ee8 1426 {
51f6de65 1427 fVertexResolution->Write();
81be4ee8 1428 fVertexResolution->ProjectionX();
1429 fVertexResolution->ProjectionY();
1430 }
0f67a57c 1431
51f6de65 1432 if (fDeltaPhi)
1433 fDeltaPhi->Write();
54b096ef 1434
a7f69e56 1435 if (fDeltaTheta)
1436 fDeltaTheta->Write();
1437
51f6de65 1438 if (fPhi)
1439 fPhi->Write();
54b096ef 1440
69b09e3b 1441 if (fRawPt)
1442 fRawPt->Write();
1443
51f6de65 1444 if (fEtaPhi)
1445 fEtaPhi->Write();
ea441adf 1446
69b09e3b 1447 for (Int_t i=0; i<2; ++i)
1448 if (fZPhi[i])
1449 fZPhi[i]->Write();
1450
12bb57f1 1451 if (fModuleMap)
1452 fModuleMap->Write();
1453
69b09e3b 1454 if (fFiredChips)
1455 fFiredChips->Write();
1456
1d107532 1457 if (fTrackletsVsClusters)
1458 fTrackletsVsClusters->Write();
1459
1460 if (fTrackletsVsUnassigned)
1461 fTrackletsVsUnassigned->Write();
1462
69b09e3b 1463 if (fStats)
1464 fStats->Write();
1465
1d107532 1466 if (fStats2)
1467 fStats2->Write();
1468
51f6de65 1469 if (fVertex)
1470 fVertex->Write();
ea441adf 1471
1d107532 1472 if (fVertexVsMult)
1473 fVertexVsMult->Write();
1474
68fa248f 1475 if (fVtxEffDen) fVtxEffDen->Write();
1476 else Printf("fVtxEffDen not found");
1477
1478 if (fVtxEffNum) fVtxEffNum->Write();
1479 else Printf("fVtxEffNum not found");
1480
1481 if (fVtxTrigEffNum) fVtxTrigEffNum->Write();
1482 else Printf("fVtxTrigEffNum not found");
1483
1484 if (fVtxTrigEffDen) fVtxTrigEffDen->Write();
1485 else Printf("fVtxTrigEffDen not found");
1486
1487 if (fTrigEffNum) fTrigEffNum->Write();
1488 else Printf("fTrigEffNum not found");
1489
1490 if (fTrigEffDen) fTrigEffDen->Write();
1491 else Printf("fTrigEffDen not found");
1492
51f6de65 1493 fout->Write();
1494 fout->Close();
0f67a57c 1495
51f6de65 1496 Printf("Writting result to analysis_esd_raw.root");
1497 }
0f67a57c 1498
1499 if (fReadMC)
1500 {
1c15d51a 1501 if (fOutput)
1502 {
1503 fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
69b09e3b 1504 fdNdEtaAnalysisND = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaND"));
1c15d51a 1505 fdNdEtaAnalysisNSD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaNSD"));
81be4ee8 1506 fdNdEtaAnalysisOnePart = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaOnePart"));
1c15d51a 1507 fdNdEtaAnalysisTr = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTr"));
1508 fdNdEtaAnalysisTrVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTrVtx"));
1509 fdNdEtaAnalysisTracks = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTracks"));
1510 fPartPt = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_pt"));
1c15d51a 1511 }
0f67a57c 1512
1513 if (!fdNdEtaAnalysis || !fdNdEtaAnalysisTr || !fdNdEtaAnalysisTrVtx || !fPartPt)
1514 {
1515 AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p %p", (void*) fdNdEtaAnalysis, (void*) fPartPt));
1516 return;
1517 }
1518
1519 fdNdEtaAnalysis->Finish(0, -1, AlidNdEtaCorrection::kNone);
69b09e3b 1520 fdNdEtaAnalysisND->Finish(0, -1, AlidNdEtaCorrection::kNone);
0fc41645 1521 fdNdEtaAnalysisNSD->Finish(0, -1, AlidNdEtaCorrection::kNone);
81be4ee8 1522 fdNdEtaAnalysisOnePart->Finish(0, -1, AlidNdEtaCorrection::kNone);
0f67a57c 1523 fdNdEtaAnalysisTr->Finish(0, -1, AlidNdEtaCorrection::kNone);
1524 fdNdEtaAnalysisTrVtx->Finish(0, -1, AlidNdEtaCorrection::kNone);
1525 fdNdEtaAnalysisTracks->Finish(0, -1, AlidNdEtaCorrection::kNone);
1526
1527 Int_t events = (Int_t) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Integral();
1528 fPartPt->Scale(1.0/events);
1529 fPartPt->Scale(1.0/fPartPt->GetBinWidth(1));
1530
51f6de65 1531 TFile* fout = new TFile("analysis_mc.root","RECREATE");
0f67a57c 1532
1533 fdNdEtaAnalysis->SaveHistograms();
69b09e3b 1534 fdNdEtaAnalysisND->SaveHistograms();
0fc41645 1535 fdNdEtaAnalysisNSD->SaveHistograms();
81be4ee8 1536 fdNdEtaAnalysisOnePart->SaveHistograms();
0f67a57c 1537 fdNdEtaAnalysisTr->SaveHistograms();
1538 fdNdEtaAnalysisTrVtx->SaveHistograms();
1539 fdNdEtaAnalysisTracks->SaveHistograms();
ea441adf 1540
1541 if (fPartPt)
1542 fPartPt->Write();
1543
0f67a57c 1544 fout->Write();
1545 fout->Close();
1546
1547 Printf("Writting result to analysis_mc.root");
1548
1549 if (fPartPt)
1550 {
1551 new TCanvas("control2", "control2", 500, 500);
770a1f1d 1552 fPartPt->Draw();
0f67a57c 1553 }
1554 }
1555}