AliPhysics  8d00e07 (8d00e07)
SPDComparison.C
Go to the documentation of this file.
1 
39 #ifdef BUILD
40 #include <AliAnalysisTaskSE.h>
41 #include <AliESDEvent.h>
42 #include <AliLog.h>
43 #include <AliMultiplicity.h>
44 #include "AliFMDEventInspector.h"
45 #include "AliAODForwardMult.h"
46 #include <TH1D.h>
47 #include <TH2D.h>
48 #include <TGraphErrors.h>
49 #include <TList.h>
50 #include <TObjArray.h>
51 #include <TMath.h>
52 
53 
66 class SPDComparisonTask : public AliAnalysisTaskSE
67 {
68 public:
74  struct VtxBin : public TNamed
75  {
84  static const char* BinName(Double_t low, Double_t high)
85  {
86  static TString n;
87  n = (Form("vtx%+05.1f_%+05.1f", low, high));
88  n.ReplaceAll("+", "p");
89  n.ReplaceAll("-", "m");
90  n.ReplaceAll(".", "d");
91  return n.Data();
92  }
99  VtxBin(Double_t low, Double_t high)
100  : TNamed(BinName(low,high), ""),
101  fClusters(0),
102  fTracklets(0)
103  {
104  }
108  VtxBin() : TNamed(), fClusters(0), fTracklets(0) {}
114  VtxBin(const VtxBin& o) :
115  TNamed(o), fClusters(o.fClusters), fTracklets(o.fTracklets)
116  {}
124  VtxBin& operator=(const VtxBin& o)
125  {
126  TNamed::operator=(o);
127  fClusters = o.fClusters;
128  fTracklets = o.fTracklets;
129  return *this;
130  }
137  void DefineOutput(TList* l, const TAxis& eta)
138  {
139  TList* ll = new TList;
140  ll->SetName(GetName());
141  ll->SetOwner();
142  l->Add(ll);
143 
144  fClusters = new TH1D("clusters", "Clusters",
145  eta.GetNbins(),
146  eta.GetXmin(),
147  eta.GetXmax());
148  fClusters->SetXTitle("#eta");
149  fClusters->SetYTitle("dN/d#eta");
150  fClusters->SetDirectory(0);
151  fClusters->Sumw2();
152  ll->Add(fClusters);
153 
154  fTracklets = static_cast<TH1D*>(fClusters->Clone("tracklets"));
155  fTracklets->SetTitle("Tracklets");
156  fTracklets->SetDirectory(0);
157  ll->Add(fTracklets);
158  }
164  void Process(const AliMultiplicity* spdmult)
165  {
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);
171  }
172 
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);
177  }
178  }
185  void Finish(TList* l, Int_t nEvents)
186  {
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);
196  ll->Add(ratio);
197  }
198  TH1D* fClusters; // Cluster histogram
199  TH1D* fTracklets; // Tracklet histogram
200  };
201 
202  //__________________________________________________________________
206  SPDComparisonTask()
207  : AliAnalysisTaskSE(),
208  fInspector(),
209  fBins(0),
210  fEvents(0),
211  fList(0),
212  fFirstEvent(true),
213  fVertexAxis(20, -20, 20)
214  {}
218  SPDComparisonTask(const char*)
219  : AliAnalysisTaskSE("spdComparision"),
220  fInspector("eventInspector"),
221  fBins(0),
222  fEvents(0),
223  fList(0),
224  fFirstEvent(true),
225  fVertexAxis(20, -20, 20)
226  {
227  // Declare our output container
228  DefineOutput(1,TList::Class());
229  }
235  SPDComparisonTask(const SPDComparisonTask& o)
236  : AliAnalysisTaskSE(o),
237  fInspector(o.fInspector),
238  fBins(o.fBins),
239  fEvents(o.fEvents),
240  fList(o.fList),
241  fFirstEvent(o.fFirstEvent),
242  fVertexAxis(20, -20, 20)
243  {
244  SetVertexAxis(o.fVertexAxis);
245  }
253  SPDComparisonTask& operator=(const SPDComparisonTask& o)
254  {
255  AliAnalysisTaskSE::operator=(o);
256  fInspector = o.fInspector;
257  fEvents = o.fEvents;
258  fList = o.fList;
259  fFirstEvent = o.fFirstEvent;
260  SetVertexAxis(o.fVertexAxis);
261  return *this;
262  }
266  ~SPDComparisonTask() {}
272  void SetVertexAxis(const TAxis& a)
273  {
274  SetVertexAxis(a.GetNbins(),
275  a.GetXmin(),
276  a.GetXmax());
277  }
285  void SetVertexAxis(Int_t n, Double_t xmin, Double_t xmax)
286  {
287  fVertexAxis.Set(n, xmin, xmax);
288  }
293  void UserCreateOutputObjects()
294  {
295  fList = new TList;
296  fList->SetName(GetName());
297  fList->SetOwner();
298 
299 
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]");
306  fList->Add(fEvents);
307 
308  TAxis& vtxAxis = fVertexAxis;
309  fBins = new TObjArray(vtxAxis.GetNbins(),1);
310 
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);
316  fBins->AddAt(bin,i);
317  }
318 
319 
320  fInspector.DefineOutput(fList);
321  fInspector.Init(*(fEvents->GetXaxis()));
322 
323  PostData(1, fList);
324  }
330  void UserExec(Option_t* option="")
331  {
332  // Get the input data - ESD event
333  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
334  if (!esd) {
335  AliWarning("No ESD event found for input event");
336  return;
337  }
338 
339  //--- Read run information -----------------------------------------
340  if (fFirstEvent && esd->GetESDRun()) {
341  fInspector.ReadRunDetails(esd);
342 
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(),
350  esd->GetBeamType(),
351  esd->GetCurrentL3(),
352  esd->GetMagneticField(),
353  esd->GetRunNumber()));
354 
355  Print();
356  fFirstEvent = false;
357  }
358 
359  // Some variables
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
366  // Process the data
367  UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ,
368  cent, nClusters);
369  Bool_t isInel = triggers & AliAODForwardMult::kB;
370  Bool_t hasVtx = retESD == AliFMDEventInspector::kOk;
371 
372  if (!isInel || !hasVtx) return;
373  fEvents->Fill(vZ);
374 
375  const AliMultiplicity* spdmult = esd->GetMultiplicity();
376 
377  VtxBin* bin = static_cast<VtxBin*>(fBins->At(iVz));
378  if (!bin) {
379  AliError(Form("No bin @ %d (%fcm)", iVz, vZ));
380  return;
381  }
382  bin->Process(spdmult);
383 
384  // Post data to output container
385  PostData(1,fList);
386  }
392  void Terminate(Option_t* option="")
393  {
394  fList = dynamic_cast<TList*>(GetOutputData(1));
395  if (!fList) {
396  AliError("No output list defined");
397  return;
398  }
399 
400  fEvents = static_cast<TH1D*>(fList->FindObject("events"));
401 
402  Int_t nEvents = fEvents->GetEntries();
403  AliInfo(Form("Got a total of %d events", nEvents));
404 
405  TH1* clusters = 0;
406  TH1* tracklets = 0;
407  TIter next(fBins);
408  VtxBin* bin = 0;
409  Int_t i = 1;
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);
416  }
417  clusters->SetDirectory(0);
418  tracklets->SetDirectory(0);
419  clusters->Scale(1. / i);
420  tracklets->Scale(1. / i);
421 
422 
423  TH1D* ratio = static_cast<TH1D*>(clusters->Clone("ratio"));
424  ratio->SetDirectory(0);
425  ratio->SetTitle("Clusters/Tracklets");
426  ratio->Divide(tracklets);
427 
428  fList->Add(clusters);
429  fList->Add(tracklets);
430  fList->Add(ratio);
431  }
432 protected:
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
438  TAxis fVertexAxis;
439 
440  ClassDef(SPDComparisonTask,1); // BF analysis
441 };
442 #else
443 
444 //====================================================================
451 void
452 SPDComparison(const char* esddir, Int_t nEvents=-1)
453 {
454  // --- Libraries to load -------------------------------------------
455  gROOT->Macro("$ALICE_PHYSICS/PWGLF/FORWARD/analysis2/scripts/LoadLibs.C");
456 
457  // --- Our data chain ----------------------------------------------
458  gROOT->LoadMacro("$ALICE_PHYSICS/PWGLF/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();
462 
463  // --- Manager -----------------------------------------------------
464  AliAnalysisManager* mgr = new AliAnalysisManager("FB", "FB train");
465  AliAnalysisManager::SetCommonFileName("spd_comps.root");
466 
467  // --- AOD input handler -------------------------------------------
468  AliESDInputHandler *inputHandler = new AliESDInputHandler();
469  mgr->SetInputEventHandler(inputHandler);
470 
471  // --- compile our code --------------------------------------------
472  gSystem->AddIncludePath("-I${ALICE_PHYSICS}/PWGLF/FORWARD/analysis2 "
473  "-I${ALICE_PHYSICS}/include "
474  "-I${ALICE_ROOT}/ANALYSIS "
475  "-I${ALICE_ROOT}/include -DBUILD=1");
476  gROOT->LoadMacro("./SPDComparison.C++g");
477 
478  // --- Make our object ---------------------------------------------
479  SPDComparisonTask* task = new SPDComparisonTask("SPD_COMP");
480  task->SetVertexAxis(10, -10, 10);
481  mgr->AddTask(task);
482 
483  // --- create containers for input/output --------------------------
484  AliAnalysisDataContainer *sums =
485  mgr->CreateContainer("spdComp", TList::Class(),
486  AliAnalysisManager::kOutputContainer,
487  AliAnalysisManager::GetCommonFileName());
488  // --- connect input/output ----------------------------------------
489  mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
490  mgr->ConnectOutput(task, 1, sums);
491 
492  // --- Run the analysis --------------------------------------------
493  TStopwatch t;
494  if (!mgr->InitAnalysis()) {
495  Error("SPDComparison", "Failed to initialize analysis train!");
496  return;
497  }
498  // Some informative output
499  mgr->PrintStatus();
500  // mgr->SetDebugLevel(3);
501  if (mgr->GetDebugLevel() < 1) mgr->SetUseProgressBar(kTRUE,100);
502 
503  // Run the train
504  t.Start();
505  Printf("=== RUNNING ANALYSIS on %9d events ==================", nEvents);
506  mgr->StartAnalysis("local", chain, nEvents);
507  t.Stop();
508  t.Print();
509 
511 }
512 
513 //====================================================================
519 void
520 DrawSPDComparison(const char* filename="spd_comps.root")
521 {
522  // --- Open the file -----------------------------------------------
523  TFile* file = TFile::Open(filename, "READ");
524  if (!file) {
525  Error("DrawSPDComparison", "Failed to open file %s", filename);
526  return;
527  }
528 
529  // --- Get our list ------------------------------------------------
530  TList* spd = static_cast<TList*>(file->Get("spdComp"));
531  if (!spd) {
532  Error("DrawSPDComparison", "Failed to get list SPD_COMP from file %s",
533  filename);
534  return;
535  }
536 
537  // --- Loop over list and get correlation plots --------------------
538  TH1* clusters = static_cast<TH1*>(spd->FindObject("clusters"));
539  TH1* tracklets = static_cast<TH1*>(spd->FindObject("tracklets"));
540  TH1* ratio = static_cast<TH1*>(spd->FindObject("ratio"));
541  TH1* events = static_cast<TH1*>(spd->FindObject("events"));
542 
543  // --- Set style parameters ----------------------------------------
544  gStyle->SetPalette(1);
545  gStyle->SetTitleFillColor(0);
546  gStyle->SetTitleStyle(0);
547  gStyle->SetTitleBorderSize(0);
548  gStyle->SetOptStat(0);
549  gStyle->SetTitleX(0.69);
550  gStyle->SetTitleY(0.99);
551  gStyle->SetTitleW(0.30);
552  gStyle->SetTitleH(0.10);
553 
554  // --- Make canvas for correlation plots and plot them -------------
555  TCanvas* c1 = new TCanvas("c1", "c1", 900, 600);
556  c1->SetTopMargin(0.05);
557  c1->SetRightMargin(0.05);
558  c1->SetFillColor(0);
559  c1->SetFillStyle(0);
560  c1->SetBorderSize(0);
561  c1->SetBorderMode(0);
562 
563  c1->cd();
564  TPad* p1 = new TPad("p1","p1", 0.0, 0.3, 1.0, 1.0, 0, 0, 0);
565  p1->SetFillColor(0);
566  p1->SetFillStyle(0);
567  p1->SetBorderSize(0);
568  p1->SetBorderMode(0);
569  p1->SetTopMargin(0.05);
570  p1->SetBottomMargin(0.001);
571  p1->SetRightMargin(0.05);
572  p1->Draw();
573  p1->cd();
574  clusters->SetMarkerStyle(20);
575  clusters->SetMarkerColor(kRed+1);
576  clusters->Draw();
577 
578  tracklets->SetMarkerStyle(20);
579  tracklets->SetMarkerColor(kBlue+1);
580  tracklets->Draw("same");
581 
582  TString v(Form("%+5.1f<|v_{z}|<%+5.1f",
583  events->GetXaxis()->GetXmin(),
584  events->GetXaxis()->GetXmax()));
585  TLatex* ltx = new TLatex(.2,.80, v.Data());
586  ltx->Draw();
587 
588  TLegend* l = p1->BuildLegend(.6,.6,.94,.94);
589  l->SetFillColor(0);
590  l->SetFillStyle(0);
591  l->SetBorderSize(0);
592 
593  c1->cd();
594  TPad* p2 = new TPad("p2","p2", 0.0, 0.0, 1.0, 0.3, 0, 0, 0);
595  p2->SetFillColor(0);
596  p2->SetFillStyle(0);
597  p2->SetBorderSize(0);
598  p2->SetBorderMode(0);
599  p2->SetTopMargin(0.001);
600  p2->SetBottomMargin(0.2);
601  p2->SetRightMargin(0.05);
602  p2->Draw();
603  p2->cd();
604  ratio->SetMarkerStyle(20);
605  ratio->SetMarkerColor(kGray+1);
606  ratio->Draw();
607  ratio->GetYaxis()->SetRangeUser(0,2);
608 
609  c1->SaveAs("SPDComparison.png");
610  c1->SaveAs("SPDComparison.root");
611 }
612 
613 
614 
615 
616 
617 #endif
618 //
619 // EOF
620 //
621 
const char * filename
Definition: TestFCM.C:1
void Print(std::ostream &o, const char *name, Double_t dT, Double_t dVM, Double_t alldT, Double_t alldVM)
Definition: PlotSysInfo.C:121
TChain * MakeChain(const char *what, const char *datadir, bool recursive=false)
Definition: MakeChain.C:151
double Double_t
Definition: External.C:58
void SPDComparison(const char *esddir, Int_t nEvents=-1)
TSystem * gSystem
Per-event per bin.
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
Definition: External.C:212
Double_t nEvents
plot quality messages
void DrawSPDComparison(const char *filename="spd_comps.root")
TFile * file
TList with histograms for a given trigger.
unsigned short UShort_t
Definition: External.C:28
const char Option_t
Definition: External.C:48
void Process(Int_t *pflag[23040][7], TH1 *inhisto, Double_t Nsigma=4., Int_t dnbins=200, Double_t dmaxval=-1., Int_t compteur=1)
bool Bool_t
Definition: External.C:53
Definition: External.C:196