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