AliPhysics  64f4410 (64f4410)
dNdyAnalysis.C
Go to the documentation of this file.
1 #ifndef DNDYANALYSIS_C
2 #define DNDYANALYSIS_C
3 #ifndef __CINT__
4 # include <TH1.h>
5 # include <TH2.h>
6 # include <TProfile.h>
7 # include <TMath.h>
8 # include <TParticle.h>
9 # include <TObjArray.h>
10 # include <TString.h>
11 # include <TGraphErrors.h>
12 # include <TMultiGraph.h>
13 #else
14 class TH1;
15 class TParticle;
16 class TMap;
17 #endif
18 #include "FastAnalysis.C"
19 #include "FastCentHelper.C"
20 namespace dNdy {
21  //====================================================================
25  struct Base : public FastAnalysis
26  {
28 
35  Base(Bool_t verbose=false, Int_t monitor=0)
36  : FastAnalysis(verbose,monitor), fdNdy(0)
37  {}
43  static TH1D* CreatedNdy()
44  {
45  Double_t maxY = 9; // 10;
46  Double_t dY = .25;
47  TH1D* h = new TH1D("dNdy", "Charged particle-rapidity density",
48  Int_t(2*maxY/dY+.5), -maxY, +maxY);
49  h->Sumw2();
50  h->SetXTitle("\\mathit{y}");
51  h->SetYTitle("\\hbox{d}N_{\\hbox{ch}}/\\hbox{d}y");
52  h->SetMarkerColor(kRed+2);
53  h->SetMarkerStyle(20);
54  h->SetDirectory(0);
55 
56  return h;
57  }
58  static TObject* CreateOutput() { return CreatedNdy(); }
62  virtual void SlaveBegin(TTree*)
63  {
64  Info("SlaveBegin", "Making dN/dy histogram");
65  fdNdy = CreatedNdy();
66  fdNdy->SetMarkerColor(kBlack);
67  fdNdy->SetMarkerStyle(21);
68  fdNdy->SetTitle(GetName());
69  fOutput->Add(fdNdy);
70  }
76  virtual Bool_t ProcessParticle(const TParticle* p)
77  {
78  Fill(p->Y());
79  return true;
80  }
81  virtual void Fill(Double_t y) { fdNdy->Fill(y); }
86  virtual void Terminate()
87  {
88  fOK = GetEventCount();
89  if (fOK <= 0) {
90  SetStatus(-1);
91  Warning("Terminate", "No events selected");
92  return;
93  }
94  fdNdy = static_cast<TH1*>(GetOutputObject("dNdy", TH1::Class()));
95  if (!fdNdy) {
96  SetStatus(-1);
97  Warning("Terminate", "No dN/deta histogram found");
98  return;
99  }
100  Printf("A total of %ld events", fOK);
101  TH1* copy = static_cast<TH1*>(fdNdy->Clone("before"));
102  fOutput->Add(copy);
103 
104  fdNdy->Scale(1. / fOK, "width");
105  fdNdy->Draw();
106  SetStatus(0);
107  }
114  {
115  TObject* m1 = new TNamed("dNdy", "");
116  m1->SetUniqueID(0x8); // Scale
117 
118  TList* ret = new TList;
119  ret->Add(m1);
120 
121  return ret;
122  }
123  ClassDef(Base,1);
124  };
125 
126  //====================================================================
131  struct NSD : public Base
132  {
139  NSD(Bool_t verbose=false, Int_t monitor=0)
140  : Base(verbose, monitor)
141  {}
149  {
150  if (fHeader->fType & 0x1) return false;
151  return true;
152  }
153  ClassDef(NSD,1);
154  };
155 
156  //====================================================================
161  struct INEL : public Base
162  {
169  INEL(Bool_t verbose=false, Int_t monitor=0)
170  : Base(verbose,monitor)
171  {}
177  virtual Bool_t ProcessHeader() { return true; }
178  ClassDef(INEL,1);
179  };
180 
181  //====================================================================
186  struct INELGt0 : public Base
187  {
194  INELGt0(Bool_t verbose=false, Int_t monitor=0)
195  : Base(verbose,monitor)
196  {}
203  {
204  if (fHeader->fType & 0x4) return true;
205  return false;
206  }
207  ClassDef(INELGt0,1);
208  };
209  //====================================================================
214  struct V0AND : public Base
215  {
222  V0AND(Bool_t verbose=false, Int_t monitor=0)
223  : Base(verbose,monitor)
224  {
225  SetTrigger("V0AND");
226  }
233  {
234  return CheckTrigger();
235  }
236  ClassDef(V0AND,1);
237  };
238 
239  //====================================================================
240  struct Cent : public Base
241  {
243  const Long_t fMinEvents;
252  Cent(const char* method="V0M", Bool_t verbose=true, Int_t monitor=0)
253  : Base(verbose,monitor),
254  fHelper(method),
255  fMinEvents(100),
256  fCentBin(0)
257  {
258  fCentMethod = fHelper.fCentMeth;
259  }
265  virtual TH1* CreateSingle() { return Base::CreatedNdy(); }
271  virtual void SlaveBegin(TTree* t)
272  {
273  Base::SlaveBegin(t);
274  // We use the parent fdNdy histogram as a cache
275  fdNdy->SetName("cache");
276  fdNdy->SetTitle("Cache");
277 
278  fHelper.CreateHistos(fOutput, Base::CreateOutput);
279  fOutput->ls();
280  }
282  {
283  if (!Base::SetupEstimator()) return false;
284 
285  fHelper.CreateDiagnostics(fOutput, fCentHist);
286  return true;
287  }
288 
292  virtual void Clear(Option_t* option="")
293  {
294  Base::Clear(option);
295  fdNdy->Reset(); // Clear the cache
296  }
303  {
304  Double_t cent = GetCentrality();
305  fCentBin = fHelper.CheckCentrality(cent,
306  fEventMult,
307  fHeader->fB,
309  fHeader->fNbin);
310  return fCentBin >= 0;
311  }
320  virtual void ProcessParticles()
321  {
322  // Check we got a bin
323  if (fCentBin < 0) return;
324 
325  // Find the histogram to update
326  TH1* out = static_cast<TH1*>(fHelper.fCentList->At(fCentBin-1));
327  // If we still have no histogram, return immediately
328  if (!out) return;
329 
330  // Use parent function to fill cache
332  if (fVerbose) {
333  Double_t n0 = fdNdy->GetBinContent(fdNdy->GetNbinsX()/2);
334  Double_t d0 = fdNdy->GetXaxis()->GetBinWidth(fdNdy->GetNbinsX()/2);
335  Double_t y0 = (n0 / d0);
336  Printf("Centrality %6.2f-%6.2f (bin %4d) "
337  "Nch=%8.1f/%4.2f=%8.1f out=%s (%d)",
338  fHelper.fCentAxis->GetBinLowEdge(fCentBin),
339  fHelper.fCentAxis->GetBinUpEdge(fCentBin),
340  fCentBin,
341  n0, d0, y0,
342  out->GetName(),
343  Int_t(out->GetBinContent(0)));
344  }
345 
346  // Make sure we have nothing in the underflow bin.
347  fdNdy->SetBinContent(0,0);
348 
349  // Add our cache to the appropriate bin
350  out->Add(fdNdy);
351  // Increment underflow bin for the event count
352  out->SetBinContent(0, out->GetBinContent(0)+1);
353  }
360  virtual void Terminate()
361  {
362  fOK = GetEventCount();
363  if (fOK <= 0) {
364  SetStatus(-1);
365  Warning("Terminate", "No events selected");
366  return;
367  }
368 
369  if (!fHelper.Finalize(fOutput, fMinEvents, Cent::Normalize))
370  SetStatus(-1);
371 
372  TH1* sigma = static_cast<TH1*>(fHelper.fCentAll->Clone("sigma"));
373  sigma->Reset();
374  sigma->SetTitle("\\sigma");
375  sigma->SetDirectory(0);
376  sigma->SetMarkerColor(kRed+2);
377  sigma->SetMarkerStyle(20);
378  sigma->SetLineColor(kRed+2);
379  fOutput->Add(sigma);
380 
381  TGraphErrors* gsigma = new TGraphErrors;
382  gsigma->SetName("gsigma");
383  gsigma->SetTitle("\\sigma");
384  gsigma->SetMarkerColor(kRed+2);
385  gsigma->SetMarkerStyle(20);
386  gsigma->SetLineColor(kRed+2);
387  fOutput->Add(gsigma);
388 
389  TGraphErrors* asigma = new TGraphErrors;
390  asigma->SetName("asigma");
391  asigma->SetTitle("\\sigma");
392  asigma->SetMarkerColor(kGreen+2);
393  asigma->SetMarkerStyle(22);
394  asigma->SetLineColor(kGreen+2);
395  fOutput->Add(asigma);
396 
397  TGraphErrors* grms = new TGraphErrors;
398  grms->SetName("grms");
399  grms->SetTitle("RMS");
400  grms->SetMarkerColor(kBlue+2);
401  grms->SetMarkerStyle(21);
402  grms->SetLineColor(kBlue+2);
403  fOutput->Add(grms);
404 
405  TProfile* nPart =
406  static_cast<TProfile*>(fOutput->FindObject("centNPartMean"));
407  if (!nPart) {
408  Warning("Terminate", "Output \"centNpartMean\" not found");
409  return;
410  }
411  Int_t j = 0;
412  for (Int_t i = 1; i <= sigma->GetNbinsX(); i++) {
413  TH1* dndy = fHelper.CentHist(i);
414  if (!dndy) {
415  // Warning("Terminate", "dNdy histogram not found for bin %d", i);
416  continue;
417  }
418  TF1* func = dndy->GetFunction("f");
419  if (!func) {
420  Warning("Terminate", "Function not found for bin %d", i);
421  continue;
422  }
423  TF1* afun = dndy->GetFunction("g");
424  if (!func) {
425  Warning("Terminate", "Function not found for bin %d", i);
426  continue;
427  }
428  Double_t sig = func->GetParameter(2);
429  Double_t esig = func->GetParError (2);
430  Double_t asig = afun->GetParameter(2);
431  Double_t easig = afun->GetParError (2);
432  Double_t nP = nPart->GetBinContent(i);
433  Double_t eP = nPart->GetBinError (i);
434  Double_t rms = dndy->GetRMS();
435  Double_t erms = dndy->GetRMSError();
436  Printf("%2d: %s %6.2f +/- %6.2f -> "
437  "%4.2f +/- %4.2f [%4.2f +/- %4.2f] (%4.2f +/- %4.2f)",
438  i, dndy->GetName(), nP, eP, sig, esig, rms, erms);
439  sigma->SetBinContent (i, sig);
440  sigma->SetBinError (i, esig);
441  gsigma->SetPoint (j, nP, sig);
442  gsigma->SetPointError(j, eP, esig);
443  asigma->SetPoint (j, nP, asig);
444  asigma->SetPointError(j, eP, easig);
445  grms ->SetPoint (j, nP, rms);
446  grms ->SetPointError(j, eP, erms);
447  j++;
448  }
449  TMultiGraph* widths = new TMultiGraph("widths", "Widths");
450  widths->Add(gsigma);
451  widths->Add(asigma);
452  widths->Add(grms);
453  fOutput->Add(widths);
454  gsigma->SaveAs("sigma.C");
455  asigma->SaveAs("asigma.C");
456  grms ->SaveAs("rms.C");
457  }
458  static TH1* Normalize(TObject* o, Int_t n)
459  {
460  if (!o->IsA()->InheritsFrom(TH1::Class())) return 0;
461  TH1* h = static_cast<TH1*>(o);
462  Int_t m = h->GetBinContent(0);
463  if (m != n)
464  ::Warning("Normalize", "Different event count here:%d helper:%d",m,n);
465  h->SetBinContent(0,0);
466  h->Scale(1. / n, "width");
467  TF1* f = new TF1("f", "gausn", -5, 5);
468  f->SetParameters(1,0,4);
469  f->FixParameter(1,0);
470  h->Fit(f,"0QBR+");
471  TF1* g = new TF1("g","gausn(0)*pol1(3)", -5, 5);
472  g->SetParameters(1,0,4,1,1);
473  g->FixParameter(0,1);
474  g->FixParameter(1,0);
475  g->SetParLimits(1,1,10);
476  h->Fit(g, "0QBR+");
477  return h;
478  }
479 
486  {
487  TObject* m1 = new TNamed("cent", "hist text30");
488  TObject* m2 = new TNamed("centAcc", "hist text30");
489  TObject* m3 = new TNamed("byCent", "e");
490 
491  m3->SetUniqueID(0x8); // Scale
492  TList* ret = new TList;
493  ret->Add(m1);
494  ret->Add(m2);
495  ret->Add(m3);
496 
497  return ret;
498  }
499  ClassDef(Cent,1);
500  };
501 
502  //====================================================================
508  struct Mult : public Cent
509  {
517  Mult(const char* method="RefMult00d80",
518  Bool_t verbose=false,
519  Int_t monitor=0)
520  : Cent(method, verbose,monitor)
521  {
522  // +1 +2 +3 +3 +5, +5, +5, +5,+10,+10,+10,+10,+10,+10,+10
523  Double_t bins[]={ 0, 3, 6, 9, 14, 19, 24, 29, 39, 49, 59, 69, 79, 89, 99};
524  fHelper.SetCentAxis(14, bins);
525  }
526  ClassDef(Mult,1);
527  };
528 } // End of namespace
529 //====================================================================
530 /*
531  * The function to make our analyser
532  */
534 {
538  dNdyMaker() : FastAnalysis::Maker("dNdy") {}
539 
550  FastAnalysis* Make(const TString& subtype,
551  Int_t monitor,
552  Bool_t verbose,
553  TMap& uopt)
554  {
555  TString t(subtype);
556  if (t.EqualTo("INEL")) return new dNdy::INEL(verbose,monitor);
557  else if (t.EqualTo("NSD")) return new dNdy::NSD(verbose,monitor);
558  else if (t.EqualTo("INELGt0")) return new dNdy::INELGt0(verbose,monitor);
559  else if (t.EqualTo("V0AND")) return new dNdy::V0AND(verbose,monitor);
560  else if (t.BeginsWith("MULT") || t.BeginsWith("CENT")) {
561  TString w(t(4, t.Length()-4));
562  if (!(w.BeginsWith("RefMult") ||
563  w.BeginsWith("ZNA") ||
564  w.BeginsWith("ZNC") ||
565  w.BeginsWith("ZPA") ||
566  w.BeginsWith("ZPC") ||
567  w.BeginsWith("V0M") ||
568  w.BeginsWith("V0A") ||
569  w.BeginsWith("V0C") ||
570  w.BeginsWith("B") ||
571  w.IsNull())) {
572  Printf("Warning: dNdyMaker::Make: Unknown estimator: %s",
573  w.Data());
574  return 0;
575  }
576  if (t.BeginsWith("MULT"))
577  return new dNdy::Mult(w, verbose, monitor);
578  else
579  return new dNdy::Cent(w, verbose, monitor);
580  }
581  Printf("Error: dNdyMaker::Run: Invalid spec: %s", t.Data());
582  return 0;
583  }
588  void List() const
589  {
590  Printf(" INEL - inelastic");
591  Printf(" INELGt0 - inelastic with at least 1 particle in |eta|<1");
592  Printf(" NSD - Non-single diffractive");
593  Printf(" V0AND - Visible X-section");
594  Printf(" CENT<est> - Centrality classes. <est> is one of ");
595  Printf(" ZNA - ZNA signal");
596  Printf(" ZNC - ZNC signal");
597  Printf(" ZPA - ZPA signal");
598  Printf(" ZPC - ZPC signal");
599  Printf(" V0M - V0-A + -C");
600  Printf(" V0A - V0-A");
601  Printf(" V0C - V0-C");
602  }
606  const char* Script() const { return __FILE__; }
607 };
608 
609 // ------------------------------------------------------------------
610 // Create instance of maker
612 
613 #ifdef __MAKECINT__
614 #pragma link C++ nestedclasses;
615 #pragma link C++ namespace dNdy;
616 #endif
617 #endif
618 //
619 // EOF
620 //
TObject * GetOutputObject(const char *name, TClass *cls)
Definition: FastAnalysis.C:490
Bool_t CheckTrigger() const
Definition: FastAnalysis.C:173
double Double_t
Definition: External.C:58
virtual void Clear(Option_t *option="")
Definition: dNdyAnalysis.C:292
FastShortHeader * fHeader
Definition: FastAnalysis.C:70
virtual Bool_t ProcessHeader()
Definition: dNdyAnalysis.C:148
INELGt0(Bool_t verbose=false, Int_t monitor=0)
Definition: dNdyAnalysis.C:194
virtual Bool_t SetupEstimator()
Definition: FastAnalysis.C:250
Bool_t Finalize(TCollection *output, Long_t minEvents, TH1 *(*callback)(TObject *, Int_t))
TH1 * CentHist(Int_t bin) const
virtual Bool_t ProcessHeader()
Definition: dNdyAnalysis.C:177
TH1 * fCentHist
Definition: FastAnalysis.C:84
virtual TList * GetMonitorObjects()
Definition: dNdyAnalysis.C:113
This script defines classes for looping over the data produced by FastSim.C.
virtual Bool_t SetupEstimator()
Definition: dNdyAnalysis.C:281
dNdyMaker * _dNdyMaker
Definition: dNdyAnalysis.C:611
virtual void ProcessParticles()
Definition: dNdyAnalysis.C:320
Double_t * sigma
FastCentHelper fHelper
Definition: dNdyAnalysis.C:242
int Int_t
Definition: External.C:63
Base(Bool_t verbose=false, Int_t monitor=0)
Definition: dNdyAnalysis.C:35
Int_t fCentBin
Definition: dNdyAnalysis.C:244
Int_t method
Definition: External.C:212
virtual Bool_t ProcessHeader()
Definition: dNdyAnalysis.C:202
virtual void Terminate()
Definition: dNdyAnalysis.C:86
virtual void SlaveBegin(TTree *)
Definition: dNdyAnalysis.C:62
virtual void SlaveBegin(TTree *t)
Definition: dNdyAnalysis.C:271
static TH1D * CreatedNdy()
Definition: dNdyAnalysis.C:43
void CreateHistos(TCollection *output, TObject *(*callback)())
Long_t GetEventCount()
Definition: FastAnalysis.C:472
Double_t GetCentrality() const
Definition: FastAnalysis.C:442
virtual void ProcessParticles()
Definition: FastAnalysis.C:561
FastAnalysis * Make(const TString &subtype, Int_t monitor, Bool_t verbose, TMap &uopt)
Definition: dNdyAnalysis.C:550
const char * Script() const
Definition: dNdyAnalysis.C:606
Bool_t fVerbose
Definition: FastAnalysis.C:74
Cent(const char *method="V0M", Bool_t verbose=true, Int_t monitor=0)
Definition: dNdyAnalysis.C:252
virtual TList * GetMonitorObjects()
Definition: dNdyAnalysis.C:485
const Long_t fMinEvents
Definition: dNdyAnalysis.C:243
ULong64_t fEventMult
Definition: FastAnalysis.C:80
ULong_t fOK
Definition: FastAnalysis.C:86
void CreateDiagnostics(TCollection *output, TH1 *centHist)
static TObject * CreateOutput()
Definition: dNdyAnalysis.C:58
Mult(const char *method="RefMult00d80", Bool_t verbose=false, Int_t monitor=0)
Definition: dNdyAnalysis.C:517
static TH1 * Normalize(TObject *o, Int_t n)
Definition: dNdyAnalysis.C:458
TH1 * fdNdy
Definition: dNdyAnalysis.C:27
virtual Bool_t ProcessHeader()
Definition: dNdyAnalysis.C:232
Int_t CheckCentrality(Double_t cent, Double_t mult, Double_t b, Double_t nPart, Double_t nBin)
virtual void Terminate()
Definition: dNdyAnalysis.C:360
INEL(Bool_t verbose=false, Int_t monitor=0)
Definition: dNdyAnalysis.C:169
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
void SetTrigger(UInt_t mask)
Definition: FastAnalysis.C:125
virtual Bool_t ProcessParticle(const TParticle *p)
Definition: dNdyAnalysis.C:76
void List() const
Definition: dNdyAnalysis.C:588
V0AND(Bool_t verbose=false, Int_t monitor=0)
Definition: dNdyAnalysis.C:222
TString fCentMethod
Definition: FastAnalysis.C:82
virtual void Fill(Double_t y)
Definition: dNdyAnalysis.C:81
NSD(Bool_t verbose=false, Int_t monitor=0)
Definition: dNdyAnalysis.C:139
Definition: External.C:196
virtual TH1 * CreateSingle()
Definition: dNdyAnalysis.C:265
virtual void Clear(Option_t *option="")
Definition: FastAnalysis.C:546
virtual Bool_t ProcessHeader()
Definition: dNdyAnalysis.C:302