2 * A task to do a comparison between tracklets and clusers in the SPD
4 * Since the class SPDComparisonTask derives from a compiled class
5 * (AliAnalysisTaskSE) we need to compile that code. The script will,
6 * when executed in the AliROOT prompt load it self again and byte
7 * compile it with the preprocessor flag BUILD defined. \
9 * The preprocessor define BUILD is not defined at first, when the
10 * script is loaded using
13 * Root> .x SPDComparison.C
16 * which means that CINT will only see the function SPDComparison.
17 * In that function, we define the BUILD preprocessor symbol
20 * gSystem->AddIncludePath("-DBUILD=1 ...");
23 * and then ACLic compile ourselves
26 * gROOT->LoadMacro("./SPDComparison.C++");
29 * But since BUILD is now defined, it means that ACLic will only see
30 * the class and not the function SPDComparison.
32 * This trick hinges on that when you initially load the script and
33 * when it is done inside the script it is done using two distinct
34 * paths - otherwise ROOT will try to unload the script first, and
40 #include <AliAnalysisTaskSE.h>
41 #include <AliESDEvent.h>
43 #include <AliMultiplicity.h>
44 #include "AliFMDEventInspector.h"
45 #include "AliAODForwardMult.h"
48 #include <TGraphErrors.h>
50 #include <TObjArray.h>
55 * @class SPDComparisonTask
57 * Task to do forward/backward correlations using FMD and SPD
60 * The class contains the sub-structure Bin that represents an
61 * @f$\eta@f$ range. One can add (possibly overlapping) @f$\eta@f$
62 * ranges by calling the member function AddBin
64 * @ingroup pwg2_forward_scripts
66 class SPDComparisonTask : public AliAnalysisTaskSE
70 * A single vertex bin - contains summed histogram of clusters and
74 struct VtxBin : public TNamed
77 * Static function to get a bin name
79 * @param low Low limit
80 * @param high High limit
82 * @return Pointer to static array
84 static const char* BinName(Double_t low, Double_t high)
87 n = (Form("vtx%+05.1f_%+05.1f", low, high));
88 n.ReplaceAll("+", "p");
89 n.ReplaceAll("-", "m");
90 n.ReplaceAll(".", "d");
96 * @param low Low limit
97 * @param high High limit
99 VtxBin(Double_t low, Double_t high)
100 : TNamed(BinName(low,high), ""),
108 VtxBin() : TNamed(), fClusters(0), fTracklets(0) {}
112 * @param o Object to copy from
114 VtxBin(const VtxBin& o) :
115 TNamed(o), fClusters(o.fClusters), fTracklets(o.fTracklets)
118 * Assignment operator
120 * @param o Object to assign from
122 * @return Reference to this.
124 VtxBin& operator=(const VtxBin& o)
126 TNamed::operator=(o);
127 fClusters = o.fClusters;
128 fTracklets = o.fTracklets;
134 * @param l Output list
135 * @param eta Eta axis to use
137 void DefineOutput(TList* l, const TAxis& eta)
139 TList* ll = new TList;
140 ll->SetName(GetName());
144 fClusters = new TH1D("clusters", "Clusters",
148 fClusters->SetXTitle("#eta");
149 fClusters->SetYTitle("dN/d#eta");
150 fClusters->SetDirectory(0);
154 fTracklets = static_cast<TH1D*>(fClusters->Clone("tracklets"));
155 fTracklets->SetTitle("Tracklets");
156 fTracklets->SetDirectory(0);
162 * @param spdmult Input
164 void Process(const AliMultiplicity* spdmult)
166 //Filling clusters in layer 1 used for tracklets...
167 for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++) {
168 Double_t eta = spdmult->GetEta(j);
169 fClusters->Fill(eta);
170 fTracklets->Fill(eta);
173 //...and then the unused ones in layer 1
174 for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++) {
175 Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.));
176 fClusters->Fill(eta);
180 * End of job processing
182 * @param l Output list
183 * @param nEvents Number of events to normalise to
185 void Finish(TList* l, Int_t nEvents)
187 TList* ll = static_cast<TList*>(l->FindObject(GetName()));
188 fClusters = static_cast<TH1D*>(ll->FindObject("clusters"));
189 fTracklets = static_cast<TH1D*>(ll->FindObject("tracklets"));
190 fClusters->Scale(1. / nEvents, "width");
191 fTracklets->Scale(1. / nEvents, "width");
192 TH1* ratio = static_cast<TH1D*>(fClusters->Clone("ratio"));
193 ratio->SetDirectory(0);
194 ratio->SetTitle("Clusters/Tracklets");
195 ratio->Divide(fTracklets);
198 TH1D* fClusters; // Cluster histogram
199 TH1D* fTracklets; // Tracklet histogram
202 //__________________________________________________________________
204 * Default constructor - do not use
207 : AliAnalysisTaskSE(),
213 fVertexAxis(20, -20, 20)
216 * Constructor with a single argument.
218 SPDComparisonTask(const char*)
219 : AliAnalysisTaskSE("spdComparision"),
220 fInspector("eventInspector"),
225 fVertexAxis(20, -20, 20)
227 // Declare our output container
228 DefineOutput(1,TList::Class());
233 * @param o Object to copy from
235 SPDComparisonTask(const SPDComparisonTask& o)
236 : AliAnalysisTaskSE(o),
237 fInspector(o.fInspector),
241 fFirstEvent(o.fFirstEvent),
242 fVertexAxis(20, -20, 20)
244 SetVertexAxis(o.fVertexAxis);
249 * @param o Object to assign from
251 * @return Reference to this
253 SPDComparisonTask& operator=(const SPDComparisonTask& o)
255 AliAnalysisTaskSE::operator=(o);
256 fInspector = o.fInspector;
259 fFirstEvent = o.fFirstEvent;
260 SetVertexAxis(o.fVertexAxis);
266 ~SPDComparisonTask() {}
268 * Set the vertex axis to use
270 * @param a Axis to set from
272 void SetVertexAxis(const TAxis& a)
274 SetVertexAxis(a.GetNbins(),
279 * Set the vertex axis to use
281 * @param n Number of bins
282 * @param xmin Least @f$v_z@f$
283 * @param xmax Most @f$v_z@f$
285 void SetVertexAxis(Int_t n, Double_t xmin, Double_t xmax)
287 fVertexAxis.Set(n, xmin, xmax);
290 * Create output objects
293 void UserCreateOutputObjects()
296 fList->SetName(GetName());
300 fEvents = new TH1D("events", "Events",
301 fVertexAxis.GetNbins(),
302 fVertexAxis.GetXmin(),
303 fVertexAxis.GetXmax());
304 fEvents->SetDirectory(0);
305 fEvents->SetXTitle("v_{z} [cm]");
308 TAxis& vtxAxis = fVertexAxis;
309 fBins = new TObjArray(vtxAxis.GetNbins(),1);
311 TAxis etaAxis(120, -3, 3);
312 for (Int_t i = 1; i <= vtxAxis.GetNbins(); i++) {
313 VtxBin* bin = new VtxBin(vtxAxis.GetBinLowEdge(i),
314 vtxAxis.GetBinUpEdge(i));
315 bin->DefineOutput(fList, etaAxis);
320 fInspector.DefineOutput(fList);
321 fInspector.Init(*(fEvents->GetXaxis()));
328 * @param option Not used
330 void UserExec(Option_t* option="")
332 // Get the input data - ESD event
333 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
335 AliWarning("No ESD event found for input event");
339 //--- Read run information -----------------------------------------
340 if (fFirstEvent && esd->GetESDRun()) {
341 fInspector.ReadRunDetails(esd);
343 AliInfo(Form("Initializing with parameters from the ESD:\n"
344 " AliESDEvent::GetBeamEnergy() ->%f\n"
345 " AliESDEvent::GetBeamType() ->%s\n"
346 " AliESDEvent::GetCurrentL3() ->%f\n"
347 " AliESDEvent::GetMagneticField()->%f\n"
348 " AliESDEvent::GetRunNumber() ->%d\n",
349 esd->GetBeamEnergy(),
352 esd->GetMagneticField(),
353 esd->GetRunNumber()));
360 UInt_t triggers; // Trigger bits
361 Bool_t lowFlux; // Low flux flag
362 UShort_t iVz; // Vertex bin from ESD
363 Double_t vZ; // Z coordinate from ESD
364 Double_t cent; // Centrality
365 UShort_t nClusters;// Number of SPD clusters
367 UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ,
369 Bool_t isInel = triggers & AliAODForwardMult::kB;
370 Bool_t hasVtx = retESD == AliFMDEventInspector::kOk;
372 if (!isInel || !hasVtx) return;
375 const AliMultiplicity* spdmult = esd->GetMultiplicity();
377 VtxBin* bin = static_cast<VtxBin*>(fBins->At(iVz));
379 AliError(Form("No bin @ %d (%fcm)", iVz, vZ));
382 bin->Process(spdmult);
384 // Post data to output container
388 * Called at the end of the processing
392 void Terminate(Option_t* option="")
394 fList = dynamic_cast<TList*>(GetOutputData(1));
396 AliError("No output list defined");
400 fEvents = static_cast<TH1D*>(fList->FindObject("events"));
402 Int_t nEvents = fEvents->GetEntries();
403 AliInfo(Form("Got a total of %d events", nEvents));
410 while ((bin = static_cast<VtxBin*>(next()))) {
411 bin->Finish(fList, fEvents->GetBinContent(i++));
412 if (!clusters) clusters = static_cast<TH1D*>(bin->fClusters->Clone());
413 else clusters->Add(bin->fClusters);
414 if (!tracklets) tracklets = static_cast<TH1D*>(bin->fTracklets->Clone());
415 else tracklets->Add(bin->fTracklets);
417 clusters->SetDirectory(0);
418 tracklets->SetDirectory(0);
419 clusters->Scale(1. / i);
420 tracklets->Scale(1. / i);
423 TH1D* ratio = static_cast<TH1D*>(clusters->Clone("ratio"));
424 ratio->SetDirectory(0);
425 ratio->SetTitle("Clusters/Tracklets");
426 ratio->Divide(tracklets);
428 fList->Add(clusters);
429 fList->Add(tracklets);
433 AliFMDEventInspector fInspector; // Inspector
434 TObjArray* fBins; // Vertex bins
435 TH1D* fEvents; // Events
436 TList* fList; // Output list
437 Bool_t fFirstEvent;// First event flag
440 ClassDef(SPDComparisonTask,1); // BF analysis
444 //====================================================================
448 * @param esddir Input directory
449 * @param nEvents Number of events, negative means all
452 SPDComparison(const char* esddir, Int_t nEvents=-1)
454 // --- Libraries to load -------------------------------------------
455 gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
457 // --- Our data chain ----------------------------------------------
458 gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/MakeChain.C");
459 TChain* chain = MakeChain("ESD", esddir, true);
460 // If 0 or less events is select, choose all
461 if (nEvents <= 0) nEvents = chain->GetEntries();
463 // --- Manager -----------------------------------------------------
464 AliAnalysisManager* mgr = new AliAnalysisManager("FB", "FB train");
465 AliAnalysisManager::SetCommonFileName("spd_comps.root");
467 // --- AOD input handler -------------------------------------------
468 AliESDInputHandler *inputHandler = new AliESDInputHandler();
469 mgr->SetInputEventHandler(inputHandler);
471 // --- compile our code --------------------------------------------
472 gSystem->AddIncludePath("-I${ALICE_ROOT}/PWG2/FORWARD/analysis2 "
473 "-I${ALICE_ROOT}/ANALYSIS "
474 "-I${ALICE_ROOT}/include -DBUILD=1");
475 gROOT->LoadMacro("./SPDComparison.C++g");
477 // --- Make our object ---------------------------------------------
478 SPDComparisonTask* task = new SPDComparisonTask("SPD_COMP");
479 task->SetVertexAxis(10, -10, 10);
482 // --- create containers for input/output --------------------------
483 AliAnalysisDataContainer *sums =
484 mgr->CreateContainer("spdComp", TList::Class(),
485 AliAnalysisManager::kOutputContainer,
486 AliAnalysisManager::GetCommonFileName());
487 // --- connect input/output ----------------------------------------
488 mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
489 mgr->ConnectOutput(task, 1, sums);
491 // --- Run the analysis --------------------------------------------
493 if (!mgr->InitAnalysis()) {
494 Error("SPDComparison", "Failed to initialize analysis train!");
497 // Some informative output
499 // mgr->SetDebugLevel(3);
500 if (mgr->GetDebugLevel() < 1) mgr->SetUseProgressBar(kTRUE,100);
504 Printf("=== RUNNING ANALYSIS on %9d events ==================", nEvents);
505 mgr->StartAnalysis("local", chain, nEvents);
512 //====================================================================
519 DrawSPDComparison(const char* filename="spd_comps.root")
521 // --- Open the file -----------------------------------------------
522 TFile* file = TFile::Open(filename, "READ");
524 Error("DrawSPDComparison", "Failed to open file %s", filename);
528 // --- Get our list ------------------------------------------------
529 TList* spd = static_cast<TList*>(file->Get("spdComp"));
531 Error("DrawSPDComparison", "Failed to get list SPD_COMP from file %s",
536 // --- Loop over list and get correlation plots --------------------
537 TH1* clusters = static_cast<TH1*>(spd->FindObject("clusters"));
538 TH1* tracklets = static_cast<TH1*>(spd->FindObject("tracklets"));
539 TH1* ratio = static_cast<TH1*>(spd->FindObject("ratio"));
540 TH1* events = static_cast<TH1*>(spd->FindObject("events"));
542 // --- Set style parameters ----------------------------------------
543 gStyle->SetPalette(1);
544 gStyle->SetTitleFillColor(0);
545 gStyle->SetTitleStyle(0);
546 gStyle->SetTitleBorderSize(0);
547 gStyle->SetOptStat(0);
548 gStyle->SetTitleX(0.69);
549 gStyle->SetTitleY(0.99);
550 gStyle->SetTitleW(0.30);
551 gStyle->SetTitleH(0.10);
553 // --- Make canvas for correlation plots and plot them -------------
554 TCanvas* c1 = new TCanvas("c1", "c1", 900, 600);
555 c1->SetTopMargin(0.05);
556 c1->SetRightMargin(0.05);
559 c1->SetBorderSize(0);
560 c1->SetBorderMode(0);
563 TPad* p1 = new TPad("p1","p1", 0.0, 0.3, 1.0, 1.0, 0, 0, 0);
566 p1->SetBorderSize(0);
567 p1->SetBorderMode(0);
568 p1->SetTopMargin(0.05);
569 p1->SetBottomMargin(0.001);
570 p1->SetRightMargin(0.05);
573 clusters->SetMarkerStyle(20);
574 clusters->SetMarkerColor(kRed+1);
577 tracklets->SetMarkerStyle(20);
578 tracklets->SetMarkerColor(kBlue+1);
579 tracklets->Draw("same");
581 TString v(Form("%+5.1f<|v_{z}|<%+5.1f",
582 events->GetXaxis()->GetXmin(),
583 events->GetXaxis()->GetXmax()));
584 TLatex* ltx = new TLatex(.2,.80, v.Data());
587 TLegend* l = p1->BuildLegend(.6,.6,.94,.94);
593 TPad* p2 = new TPad("p2","p2", 0.0, 0.0, 1.0, 0.3, 0, 0, 0);
596 p2->SetBorderSize(0);
597 p2->SetBorderMode(0);
598 p2->SetTopMargin(0.001);
599 p2->SetBottomMargin(0.2);
600 p2->SetRightMargin(0.05);
603 ratio->SetMarkerStyle(20);
604 ratio->SetMarkerColor(kGray+1);
606 ratio->GetYaxis()->SetRangeUser(0,2);
608 c1->SaveAs("SPDComparison.png");
609 c1->SaveAs("SPDComparison.root");