3 #include "AlidNdEtaTask.h"
14 #include <TParticle.h>
17 #include <TObjString.h>
21 #include <AliESDVertex.h>
22 #include <AliESDEvent.h>
24 #include <AliHeader.h>
25 #include <AliGenEventHeader.h>
26 #include <AliMultiplicity.h>
27 #include <AliAnalysisManager.h>
28 #include <AliMCEventHandler.h>
29 #include <AliMCEvent.h>
30 #include <AliESDInputHandler.h>
32 #include "AliESDtrackCuts.h"
33 #include "AliPWG0Helper.h"
34 #include "AliCorrection.h"
35 #include "AliCorrectionMatrix3D.h"
36 #include "dNdEta/dNdEtaAnalysis.h"
38 ClassImp(AlidNdEtaTask)
40 AlidNdEtaTask::AlidNdEtaTask(const char* opt) :
41 AliAnalysisTask("AlidNdEtaTask", ""),
45 fAnalysisMode(AliPWG0Helper::kTPC),
46 fTrigger(AliPWG0Helper::kMB1),
51 fdNdEtaAnalysisESD(0),
57 fdNdEtaAnalysisNSD(0),
59 fdNdEtaAnalysisTrVtx(0),
60 fdNdEtaAnalysisTracks(0),
68 // Constructor. Initialization of pointers
71 // Define input and output slots here
72 DefineInput(0, TChain::Class());
73 DefineOutput(0, TList::Class());
75 AliLog::SetClassDebugLevel("AlidNdEtaTask", AliLog::kWarning);
78 AlidNdEtaTask::~AlidNdEtaTask()
84 // histograms are in the output list and deleted when the output
85 // list is deleted by the TSelector dtor
93 Bool_t AlidNdEtaTask::Notify()
95 static Int_t count = 0;
97 Printf("Processing %d. file", count);
101 //________________________________________________________________________
102 void AlidNdEtaTask::ConnectInputData(Option_t *)
107 Printf("AlidNdEtaTask::ConnectInputData called");
109 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
112 Printf("ERROR: Could not get ESDInputHandler");
114 fESD = esdH->GetEvent();
116 TTree* tree = esdH->GetTree();
118 Printf("ERROR: Could not read tree");
120 // Disable all branches and enable only the needed ones
121 tree->SetBranchStatus("*", 0);
123 tree->SetBranchStatus("AliESDHeader*", 1);
124 tree->SetBranchStatus("*Vertex*", 1);
126 if (fAnalysisMode == AliPWG0Helper::kSPD) {
127 tree->SetBranchStatus("AliMultiplicity*", 1);
130 if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS) {
131 //AliESDtrackCuts::EnableNeededBranches(tree);
132 tree->SetBranchStatus("Tracks*", 1);
138 // disable info messages of AliMCEvent (per event)
139 AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
142 void AlidNdEtaTask::CreateOutputObjects()
144 // create result objects and add to output list
146 Printf("AlidNdEtaTask::CreateOutputObjects");
151 fdNdEtaAnalysisESD = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD", fAnalysisMode);
152 fOutput->Add(fdNdEtaAnalysisESD);
154 fMult = new TH1F("fMult", "fMult;Ntracks;Count", 201, -0.5, 200.5);
157 fMultVtx = new TH1F("fMultVtx", "fMultVtx;Ntracks;Count", 201, -0.5, 200.5);
158 fOutput->Add(fMultVtx);
160 for (Int_t i=0; i<3; ++i)
162 fPartEta[i] = new TH1F(Form("dndeta_check_%d", i), Form("dndeta_check_%d", i), 60, -6, 6);
163 fPartEta[i]->Sumw2();
164 fOutput->Add(fPartEta[i]);
167 fEvents = new TH1F("dndeta_check_vertex", "dndeta_check_vertex", 160, -40, 40);
168 fOutput->Add(fEvents);
170 fVertexResolution = new TH1F("dndeta_vertex_resolution_z", "dndeta_vertex_resolution_z", 1000, 0, 10);
171 fOutput->Add(fVertexResolution);
173 fPhi = new TH1F("fPhi", "fPhi;#phi in rad.;count", 720, 0, 2 * TMath::Pi());
176 fEtaPhi = new TH2F("fEtaPhi", "fEtaPhi;#eta;#phi in rad.;count", 80, -4, 4, 18*5, 0, 2 * TMath::Pi());
177 fOutput->Add(fEtaPhi);
179 if (fAnalysisMode == AliPWG0Helper::kSPD)
180 fDeltaPhi = new TH1F("fDeltaPhi", "fDeltaPhi;#Delta #phi;Entries", 18*50, -3.14, 3.14);
184 fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta", fAnalysisMode);
185 fOutput->Add(fdNdEtaAnalysis);
187 fdNdEtaAnalysisNSD = new dNdEtaAnalysis("dndetaNSD", "dndetaNSD", fAnalysisMode);
188 fOutput->Add(fdNdEtaAnalysisNSD);
190 fdNdEtaAnalysisTr = new dNdEtaAnalysis("dndetaTr", "dndetaTr", fAnalysisMode);
191 fOutput->Add(fdNdEtaAnalysisTr);
193 fdNdEtaAnalysisTrVtx = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx", fAnalysisMode);
194 fOutput->Add(fdNdEtaAnalysisTrVtx);
196 fdNdEtaAnalysisTracks = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks", fAnalysisMode);
197 fOutput->Add(fdNdEtaAnalysisTracks);
199 fVertex = new TH3F("vertex_check", "vertex_check", 50, -50, 50, 50, -50, 50, 50, -50, 50);
200 fOutput->Add(fVertex);
202 fPartPt = new TH1F("dndeta_check_pt", "dndeta_check_pt", 1000, 0, 10);
204 fOutput->Add(fPartPt);
209 fEsdTrackCuts->SetName("fEsdTrackCuts");
210 fOutput->Add(fEsdTrackCuts);
214 void AlidNdEtaTask::Exec(Option_t*)
218 // these variables are also used in the MC section below; however, if ESD is off they stay with their default values
219 Bool_t eventTriggered = kFALSE;
220 const AliESDVertex* vtxESD = 0;
225 // trigger definition
226 eventTriggered = AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask(), fTrigger);
228 // get the ESD vertex
229 vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
238 // fill z vertex resolution
240 fVertexResolution->Fill(vtxESD->GetZRes());
242 // post the data already here
243 PostData(0, fOutput);
245 // needed for syst. studies
249 if (fUseMCVertex || fUseMCKine || fReadMC) {
251 Printf("ERROR: fUseMCVertex or fUseMCKine set without fReadMC set!");
255 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
257 Printf("ERROR: Could not retrieve MC event handler");
261 AliMCEvent* mcEvent = eventHandler->MCEvent();
263 Printf("ERROR: Could not retrieve MC event");
267 AliHeader* header = mcEvent->Header();
270 AliDebug(AliLog::kError, "Header not available");
275 AliGenEventHeader* genHeader = header->GenEventHeader();
278 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
281 genHeader->PrimaryVertex(vtxMC);
285 Printf("WARNING: Replacing vertex by MC vertex. This is for systematical checks only.");
289 stack = mcEvent->Stack();
292 AliDebug(AliLog::kError, "Stack not available");
297 // create list of (label, eta, pt) tuples
298 Int_t inputCount = 0;
302 if (fAnalysisMode == AliPWG0Helper::kSPD)
305 const AliMultiplicity* mult = fESD->GetMultiplicity();
308 AliDebug(AliLog::kError, "AliMultiplicity not available");
312 labelArr = new Int_t[mult->GetNumberOfTracklets()];
313 etaArr = new Float_t[mult->GetNumberOfTracklets()];
314 ptArr = new Float_t[mult->GetNumberOfTracklets()];
317 Printf("Processing only primaries (MC information used). This is for systematical checks only.");
319 // get multiplicity from ITS tracklets
320 for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
322 //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
325 if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
328 Float_t deltaPhi = mult->GetDeltaPhi(i);
329 // prevent values to be shifted by 2 Pi()
330 if (deltaPhi < -TMath::Pi())
331 deltaPhi += TMath::Pi() * 2;
332 if (deltaPhi > TMath::Pi())
333 deltaPhi -= TMath::Pi() * 2;
335 if (TMath::Abs(deltaPhi) > 1)
336 printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
338 Float_t phi = mult->GetPhi(i);
340 phi += TMath::Pi() * 2;
342 fEtaPhi->Fill(mult->GetEta(i), phi);
344 fDeltaPhi->Fill(deltaPhi);
346 etaArr[inputCount] = mult->GetEta(i);
347 labelArr[inputCount] = mult->GetLabel(i, 0);
348 ptArr[inputCount] = 0; // no pt for tracklets
352 //Printf("Accepted %d tracklets", inputCount);
354 else if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
358 AliDebug(AliLog::kError, "fESDTrackCuts not available");
362 // get multiplicity from ESD tracks
363 TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
364 Int_t nGoodTracks = list->GetEntries();
366 labelArr = new Int_t[nGoodTracks];
367 etaArr = new Float_t[nGoodTracks];
368 ptArr = new Float_t[nGoodTracks];
370 // loop over esd tracks
371 for (Int_t i=0; i<nGoodTracks; i++)
373 AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
376 AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
380 Float_t phi = esdTrack->Phi();
382 phi += TMath::Pi() * 2;
384 fEtaPhi->Fill(esdTrack->Eta(), phi);
386 etaArr[inputCount] = esdTrack->Eta();
387 labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
388 ptArr[inputCount] = esdTrack->Pt();
392 Printf("Accepted %d tracks", nGoodTracks);
399 // Processing of ESD information (always)
403 fMult->Fill(inputCount);
404 fdNdEtaAnalysisESD->FillTriggeredEvent(inputCount);
409 fMultVtx->Fill(inputCount);
411 for (Int_t i=0; i<inputCount; ++i)
413 Float_t eta = etaArr[i];
414 Float_t pt = ptArr[i];
416 fdNdEtaAnalysisESD->FillTrack(vtx[2], eta, pt);
418 if (TMath::Abs(vtx[2]) < 20)
420 fPartEta[0]->Fill(eta);
423 fPartEta[1]->Fill(eta);
425 fPartEta[2]->Fill(eta);
429 // for event count per vertex
430 fdNdEtaAnalysisESD->FillEvent(vtx[2], inputCount);
433 fEvents->Fill(vtx[2]);
437 // from tracks is only done for triggered and vertex reconstructed events
438 for (Int_t i=0; i<inputCount; ++i)
440 Int_t label = labelArr[i];
445 //Printf("Getting particle of track %d", label);
446 TParticle* particle = stack->Particle(label);
450 AliDebug(AliLog::kError, Form("ERROR: Could not retrieve particle %d.", label));
454 fdNdEtaAnalysisTracks->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
455 } // end of track loop
457 // for event count per vertex
458 fdNdEtaAnalysisTracks->FillEvent(vtxMC[2], inputCount);
468 if (fReadMC) // Processing of MC information (optional)
470 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
472 Printf("ERROR: Could not retrieve MC event handler");
476 AliMCEvent* mcEvent = eventHandler->MCEvent();
478 Printf("ERROR: Could not retrieve MC event");
482 AliStack* stack = mcEvent->Stack();
485 AliDebug(AliLog::kError, "Stack not available");
489 AliHeader* header = mcEvent->Header();
492 AliDebug(AliLog::kError, "Header not available");
497 AliGenEventHeader* genHeader = header->GenEventHeader();
500 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
505 genHeader->PrimaryVertex(vtxMC);
508 Int_t processType = AliPWG0Helper::GetEventProcessType(header);
509 AliDebug(AliLog::kDebug+1, Form("Found process type %d", processType));
511 if (processType==AliPWG0Helper::kInvalidProcess)
512 AliDebug(AliLog::kError, Form("Unknown process type %d.", processType));
514 // loop over mc particles
515 Int_t nPrim = stack->GetNprimary();
517 Int_t nAcceptedParticles = 0;
519 for (Int_t iMc = 0; iMc < nPrim; ++iMc)
521 if (!stack->IsPhysicalPrimary(iMc))
524 //Printf("Getting particle %d", iMc);
525 TParticle* particle = stack->Particle(iMc);
530 if (particle->GetPDG()->Charge() == 0)
533 //AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID()));
534 Float_t eta = particle->Eta();
535 Float_t pt = particle->Pt();
537 // make a rough eta cut (so that nAcceptedParticles is not too far off the true measured value (NB: this histograms are only gathered for comparison))
538 if (TMath::Abs(eta) < 1.5) // && pt > 0.3)
539 nAcceptedParticles++;
541 fdNdEtaAnalysis->FillTrack(vtxMC[2], eta, pt);
542 fVertex->Fill(particle->Vx(), particle->Vy(), particle->Vz());
544 if (processType != AliPWG0Helper::kSD )
545 fdNdEtaAnalysisNSD->FillTrack(vtxMC[2], eta, pt);
549 fdNdEtaAnalysisTr->FillTrack(vtxMC[2], eta, pt);
551 fdNdEtaAnalysisTrVtx->FillTrack(vtxMC[2], eta, pt);
554 if (TMath::Abs(eta) < 0.8)
558 fdNdEtaAnalysis->FillEvent(vtxMC[2], nAcceptedParticles);
559 if (processType != AliPWG0Helper::kSD)
560 fdNdEtaAnalysisNSD->FillEvent(vtxMC[2], nAcceptedParticles);
564 fdNdEtaAnalysisTr->FillEvent(vtxMC[2], nAcceptedParticles);
566 fdNdEtaAnalysisTrVtx->FillEvent(vtxMC[2], nAcceptedParticles);
571 void AlidNdEtaTask::Terminate(Option_t *)
573 // The Terminate() function is the last function to be called during
574 // a query. It always runs on the client, it can be used to present
575 // the results graphically or save the results to file.
577 fOutput = dynamic_cast<TList*> (GetOutputData(0));
579 Printf("ERROR: fOutput not available");
583 fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("fdNdEtaAnalysisESD"));
584 fMult = dynamic_cast<TH1F*> (fOutput->FindObject("fMult"));
585 fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx"));
586 for (Int_t i=0; i<3; ++i)
587 fPartEta[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("dndeta_check_%d", i)));
588 fEvents = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_vertex"));
589 fVertexResolution = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_vertex_resolution_z"));
591 fPhi = dynamic_cast<TH1F*> (fOutput->FindObject("fPhi"));
592 fEtaPhi = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaPhi"));
593 fDeltaPhi = dynamic_cast<TH1F*> (fOutput->FindObject("fDeltaPhi"));
595 fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCuts"));
598 if (!fdNdEtaAnalysisESD)
600 AliDebug(AliLog::kError, "ERROR: fdNdEtaAnalysisESD not available");
604 if (fMult && fMultVtx)
608 fMultVtx->SetLineColor(2);
609 fMultVtx->Draw("SAME");
614 Int_t events1 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(-19.9), fEvents->GetXaxis()->FindBin(-0.001));
615 Int_t events2 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(0.001), fEvents->GetXaxis()->FindBin(19.9));
617 Printf("%d events with vertex used", events1 + events2);
619 if (events1 > 0 && events2 > 0)
621 fPartEta[0]->Scale(1.0 / (events1 + events2));
622 fPartEta[1]->Scale(1.0 / events1);
623 fPartEta[2]->Scale(1.0 / events2);
625 for (Int_t i=0; i<3; ++i)
626 fPartEta[i]->Scale(1.0 / fPartEta[i]->GetBinWidth(1));
628 new TCanvas("control", "control", 500, 500);
629 for (Int_t i=0; i<3; ++i)
631 fPartEta[i]->SetLineColor(i+1);
632 fPartEta[i]->Draw((i==0) ? "" : "SAME");
639 new TCanvas("control3", "control3", 500, 500);
643 TFile* fout = new TFile("analysis_esd_raw.root", "RECREATE");
645 if (fdNdEtaAnalysisESD)
646 fdNdEtaAnalysisESD->SaveHistograms();
649 fEsdTrackCuts->SaveHistograms("esd_tracks_cuts");
657 for (Int_t i=0; i<3; ++i)
659 fPartEta[i]->Write();
664 if (fVertexResolution)
665 fVertexResolution->Write();
679 Printf("Writting result to analysis_esd_raw.root");
685 fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
686 fdNdEtaAnalysisNSD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaNSD"));
687 fdNdEtaAnalysisTr = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTr"));
688 fdNdEtaAnalysisTrVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTrVtx"));
689 fdNdEtaAnalysisTracks = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTracks"));
690 fPartPt = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_pt"));
691 fVertex = dynamic_cast<TH3F*> (fOutput->FindObject("vertex_check"));
694 if (!fdNdEtaAnalysis || !fdNdEtaAnalysisTr || !fdNdEtaAnalysisTrVtx || !fPartPt)
696 AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p %p", (void*) fdNdEtaAnalysis, (void*) fPartPt));
700 fdNdEtaAnalysis->Finish(0, -1, AlidNdEtaCorrection::kNone);
701 fdNdEtaAnalysisNSD->Finish(0, -1, AlidNdEtaCorrection::kNone);
702 fdNdEtaAnalysisTr->Finish(0, -1, AlidNdEtaCorrection::kNone);
703 fdNdEtaAnalysisTrVtx->Finish(0, -1, AlidNdEtaCorrection::kNone);
704 fdNdEtaAnalysisTracks->Finish(0, -1, AlidNdEtaCorrection::kNone);
706 Int_t events = (Int_t) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Integral();
707 fPartPt->Scale(1.0/events);
708 fPartPt->Scale(1.0/fPartPt->GetBinWidth(1));
710 fout = new TFile("analysis_mc.root","RECREATE");
712 fdNdEtaAnalysis->SaveHistograms();
713 fdNdEtaAnalysisNSD->SaveHistograms();
714 fdNdEtaAnalysisTr->SaveHistograms();
715 fdNdEtaAnalysisTrVtx->SaveHistograms();
716 fdNdEtaAnalysisTracks->SaveHistograms();
727 Printf("Writting result to analysis_mc.root");
731 new TCanvas("control2", "control2", 500, 500);