dNdyAnalysis.C

Go to the documentation of this file.

#ifndef __CINT__
# include <TH1.h>
# include <TH2.h>
# include <TMath.h>
# include <TParticle.h>
# include <TObjArray.h>
# include <TString.h>
#else
class TH1;
class TParticle;
#endif
#include "FastAnalysis.C"
#include "FastCentHelper.C"
namespace dNdy {
  //====================================================================
  struct Base : public FastAnalysis
  {
    TH1* fdNdy; 
  
    Base(Bool_t verbose=false, Int_t monitor=0)
      : FastAnalysis(verbose,monitor), fdNdy(0)
    {}
    static TH1D* CreatedNdy()
    {
      Double_t maxY   = 9; // 10;
      Double_t dY     = .25; 
      TH1D*    h      = new TH1D("dNdy", "Charged particle-rapidity density",
         Int_t(2*maxY/dY+.5), -maxY, +maxY);
      h->Sumw2();
      h->SetXTitle("\\mathit{y}");
      h->SetYTitle("\\hbox{d}N_{\\hbox{ch}}/\\hbox{d}y");
      h->SetMarkerColor(kRed+2);
      h->SetMarkerStyle(20);
      h->SetDirectory(0);
    
      return h;
    }
    virtual void SlaveBegin(TTree*)
    {
      Info("SlaveBegin", "Making dN/dy histogram");
      fdNdy = CreatedNdy();
      fdNdy->SetMarkerColor(kBlack);
      fdNdy->SetMarkerStyle(21);
      fdNdy->SetTitle(GetName());
      fOutput->Add(fdNdy);
    }
    virtual Bool_t ProcessParticle(const TParticle* p)
    {
      Fill(p->Y());
      return true;
    }
    virtual void Fill(Double_t y) { fdNdy->Fill(y); }
    virtual void Terminate()
    {
      fOK = GetEventCount();
      if (fOK <= 0) {
  SetStatus(-1);
  Warning("Terminate", "No events selected");
  return;
      }
      fdNdy = static_cast<TH1*>(GetOutputObject("dNdy", TH1::Class()));
      if (!fdNdy) {
  SetStatus(-1);
  Warning("Terminate", "No dN/deta histogram found");
  return;
      }
      Printf("A total of %ld events", fOK);
      TH1* copy = static_cast<TH1*>(fdNdy->Clone("before"));
      fOutput->Add(copy);
            
      fdNdy->Scale(1. / fOK, "width");
      fdNdy->Draw();
      SetStatus(0);
    }
    virtual TList* GetMonitorObjects()
    {
      TObject* m1 = new TNamed("dNdy", "");
      m1->SetUniqueID(0x8); // Scale 

      TList* ret = new TList;
      ret->Add(m1);
    
      return ret;
    }
    ClassDef(Base,1);
  };

  //====================================================================
  struct NSD : public Base
  {
    NSD(Bool_t verbose=false, Int_t monitor=0)
      : Base(verbose, monitor)
    {}
    virtual Bool_t ProcessHeader()
    {
      if (fHeader->fType & 0x1) return false;
      return true;
    }
    ClassDef(NSD,1);
  };

  //====================================================================
  struct INEL : public Base
  {
    INEL(Bool_t verbose=false, Int_t monitor=0)
      : Base(verbose,monitor)
    {}
    virtual Bool_t ProcessHeader() { return true; }
    ClassDef(INEL,1);
  };

  //====================================================================
  struct INELGt0 : public Base
  {
    INELGt0(Bool_t verbose=false, Int_t monitor=0)
      : Base(verbose,monitor)
    {}
    virtual Bool_t ProcessHeader()
    {
      if (fHeader->fType & 0x4) return true;
      return false;
    }
    ClassDef(INELGt0,1);
  };

  //====================================================================
  struct Cent : public Base 
  {
    FastCentHelper fHelper;
    const Long_t   fMinEvents;
    Int_t          fCentBin;
    Cent(const char* method="V0M", Bool_t verbose=true, Int_t monitor=0)
      : Base(verbose,monitor),
  fHelper(method),
  fMinEvents(100),
  fCentBin(0) 
    {
      fCentMethod = fHelper.fCentMeth;
    }
    virtual TH1* CreateSingle() { return Base::CreatedNdy(); }
    virtual void SlaveBegin(TTree* t)
    {
      Base::SlaveBegin(t);
      // We use the parent fdNdy histogram as a cache 
      fdNdy->SetName("cache");
      fdNdy->SetTitle("Cache");

      fHelper.CreateHistos(fOutput, Base::CreatedNdy);
      fOutput->ls();
    }
    virtual Bool_t SetupEstimator()
    {
      if (!Base::SetupEstimator()) return false;
    
      fHelper.CreateDiagnostics(fOutput, fCentHist);
      return true;
    }    
         
    virtual void Clear(Option_t* option="") 
    {
      Base::Clear(option);
      fdNdy->Reset(); // Clear the cache
    }
    virtual Bool_t ProcessHeader()
    {
      Double_t cent = GetCentrality();
      fCentBin      = fHelper.CheckCentrality(cent,
                fEventMult,
                fHeader->fB,
                fHeader->fNtgt+fHeader->fNproj,
                fHeader->fNbin);
      return fCentBin >= 0;
    }
    virtual void ProcessParticles()
    {
      // Check we got a bin 
      if (fCentBin < 0) return;

      // Find the histogram to update 
      TH1*     out  = static_cast<TH1*>(fHelper.fCentList->At(fCentBin-1));
      // If we still have no histogram, return immediately 
      if (!out) return;

      // Use parent function to fill cache 
      Base::ProcessParticles();
      if (fVerbose) {
  Double_t n0   = fdNdy->GetBinContent(fdNdy->GetNbinsX()/2);
  Double_t d0   = fdNdy->GetXaxis()->GetBinWidth(fdNdy->GetNbinsX()/2);
  Double_t y0   = (n0 / d0);
  Printf("Centrality %6.2f-%6.2f (bin %4d) "
         "Nch=%8.1f/%4.2f=%8.1f out=%s (%d)",
         fHelper.fCentAxis->GetBinLowEdge(fCentBin),
         fHelper.fCentAxis->GetBinUpEdge(fCentBin),
         fCentBin,
         n0, d0, y0,
         out->GetName(),
         Int_t(out->GetBinContent(0)));
      }

      // Make sure we have nothing in the underflow bin. 
      fdNdy->SetBinContent(0,0);

      // Add our cache to the appropriate bin 
      out->Add(fdNdy);
      // Increment underflow bin for the event count 
      out->SetBinContent(0, out->GetBinContent(0)+1);
    }
    virtual void Terminate()
    {
      fOK = GetEventCount();
      if (fOK <= 0) {
  SetStatus(-1);
  Warning("Terminate", "No events selected");
  return;
      }

      if (!fHelper.Finalize(fOutput, fMinEvents))
  SetStatus(-1);
    }
    virtual TList* GetMonitorObjects()
    {
      TObject* m1 = new TNamed("cent",     "hist text30");
      TObject* m2 = new TNamed("centAcc",  "hist text30");
      TObject* m3 = new TNamed("byCent",   "e");
    
      m3->SetUniqueID(0x8); // Scale 
      TList* ret = new TList;
      ret->Add(m1);
      ret->Add(m2);
      ret->Add(m3);
    
      return ret;
    }
    ClassDef(Cent,1);
  };
  
  //====================================================================
  struct Mult : public Cent
  {
    Mult(const char* method="RefMult00d80",
   Bool_t verbose=false,
   Int_t monitor=0)
      : Cent(method, verbose,monitor)
    {
      //              +1 +2 +3 +3  +5, +5, +5, +5,+10,+10,+10,+10,+10,+10,+10
      Double_t bins[]={ 0, 3, 6, 9, 14, 19, 24, 29, 39, 49, 59, 69, 79, 89, 99};
      fHelper.SetCentAxis(14, bins);
    }
    ClassDef(Mult,1);
  };
} // End of namespace 
//====================================================================
/*
 * The function to make our analyser 
 */
struct dNdyMaker : public FastAnalysis::Maker
{
  dNdyMaker() : FastAnalysis::Maker("dNdy") {}
  
  FastAnalysis*  Make(const TString& subtype,
          Int_t          monitor,
          Bool_t         verbose,
          TString&       uopt)
  {
    TString t(subtype);
    if      (t.EqualTo("INEL"))    return new dNdy::INEL(verbose,monitor);
    else if (t.EqualTo("NSD"))     return new dNdy::NSD(verbose,monitor);
    else if (t.EqualTo("INELGt0")) return new dNdy::INELGt0(verbose,monitor);
    else if (t.BeginsWith("MULT") || t.BeginsWith("CENT")) {
      TString w(t(4, t.Length()-4));
      if (!(w.BeginsWith("RefMult") ||
      w.BeginsWith("ZNA") || 
      w.BeginsWith("ZNC") || 
      w.BeginsWith("ZPA") || 
      w.BeginsWith("ZPC") || 
      w.BeginsWith("V0M") ||
      w.BeginsWith("V0A") ||
      w.BeginsWith("V0C") ||
      w.BeginsWith("B")   ||
      w.IsNull())) {
  Printf("Warning: dNdyMaker::Make: Unknown estimator: %s",
         w.Data());
  return 0;
      }
      if (t.BeginsWith("MULT"))
  return new dNdy::Mult(w, verbose, monitor);
      else
  return new dNdy::Cent(w, verbose, monitor);
    }
    Printf("Error: dNdyMaker::Run: Invalid spec: %s", t.Data());
    return 0;
  }
  void List() const
  {
    Printf(" INEL            - inelastic");
    Printf(" INELGt0         - inelastic with at least 1 particle in |eta|<1");
    Printf(" NSD             - Non-single diffractive");
    Printf(" CENT<est>       - Centrality classes. <est> is one of ");
    Printf("  ZNA            - ZNA signal");
    Printf("  ZNC            - ZNC signal");
    Printf("  ZPA            - ZPA signal");
    Printf("  ZPC            - ZPC signal");
    Printf("  V0M            - V0-A + -C");
    Printf("  V0A            - V0-A");
    Printf("  V0C            - V0-C");
  }
  const char* Script() const {   return __FILE__;  }
};

// ------------------------------------------------------------------
// Create instance of maker
dNdyMaker* _dNdyMaker = new dNdyMaker;

#ifdef __MAKECINT__
#pragma link C++ nestedclasses;
#pragma link C++ namespace dNdy;
#endif

//
// EOF
//