PWGLF/FORWARD/analysis2/AliForwardMultiplicityDistribution.cxx

   1
   2 #include <TH1D.h>
   3 #include "AliForwardMultiplicityDistribution.h"
   4 #include "AliAODForwardMult.h"
   5 #include "AliAODCentralMult.h"
   6 #include "AliAODEvent.h"
   7 #include "TString.h"
   8 using namespace std;
   9
  10 ClassImp(AliForwardMultiplicityDistribution)
  11 #if 0
  12 ; // This is for Emacs - do not delete
  13 #endif
  14 //______________________________________________________________________
  15 AliForwardMultiplicityDistribution::AliForwardMultiplicityDistribution()
  16   : AliBasedNdetaTask(),
  17     fTrigger(0),   // trigger histogram
  18     fBins(),       // eta bin list
  19     fOutput(0),    // output list
  20     fLowCent(-1),  // lower centrality limit
  21     fHighCent(-1), // upper centrality limit
  22     fNBins(-1),    // multiplicity axis' runs from 0 to fNbins
  23     fCent(0)       // centrality
  24 {
  25   //
  26   // Default Constructor
  27   //
  28 }
  29
  30 //_____________________________________________________________________
  31 AliForwardMultiplicityDistribution::AliForwardMultiplicityDistribution(const char* name)
  32   : AliBasedNdetaTask(name),
  33     fTrigger(0),   // trigger histogram
  34     fBins(),       // eta bin list
  35     fOutput(0),    // output list
  36     fLowCent(-1),  // lower centrality limit
  37     fHighCent(-1), // upper centrality limit
  38     fNBins(-1),    // multiplicity axis' runs from 0 to fNbins
  39     fCent(0)       // centrality
  40 {
  41   //
  42   // Constructor
  43   //
  44   DefineOutput(1, TList::Class());
  45 }
  46
  47 //_____________________________________________________________________
  48 void AliForwardMultiplicityDistribution::UserCreateOutputObjects()
  49 {
  50   //
  51   // Create Output Objects
  52   //
  53   fOutput = new TList;
  54   fOutput->SetOwner();
  55   fOutput->SetName(GetName());
  56
  57   TH1D* dndetaSumForward = new TH1D("dndetaSumForward","dndetaSumForward", 200,-4,6);
  58   TH1D* dndetaSumCentral = new TH1D("dndetaSumCentral","dndetaSumCentral", 200,-4,6);
  59   fOutput->Add(dndetaSumForward);
  60   fOutput->Add(dndetaSumCentral);
  61
  62   fCent = new TH1D("cent", "Centrality", 100, 0, 100);
  63   fCent->SetDirectory(0);
  64   fCent->SetXTitle(0);
  65   fOutput->Add(fCent);
  66
  67   TH1D* dndetaEventForward = new TH1D();
  68   TH1D* dndetaEventCentral = new TH1D();
  69   fOutput->Add(dndetaEventForward);
  70   fOutput->Add(dndetaEventCentral);
  71
  72   //make trigger diagnostics histogram
  73   fTrigger= new TH1I();
  74   fTrigger= AliAODForwardMult::MakeTriggerHistogram("triggerHist");
  75   fOutput->Add(fTrigger);
  76
  77   //Loop over all individual eta bins, and define their hisograms
  78   TIter next(&fBins);
  79   Bin * bin = 0;
  80   while ((bin = static_cast<Bin*>(next()))) {
  81     bin->DefineOutputs(fOutput, fNBins);
  82   }
  83   //fOutput->ls();
  84   PostData(1, fOutput);
  85 }
  86
  87
  88 //_____________________________________________________________________
  89 void AliForwardMultiplicityDistribution::UserExec(Option_t */*option*/)
  90 {
  91   //
  92   // User Exec
  93   //
  94   AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
  95   if (!aod) {
  96     AliError("Cannot get the AOD event");
  97     return;
  98   }
  99   // Check AOD event
 100   AliAODForwardMult* AODforward = static_cast<AliAODForwardMult*>(aod->FindListObject("Forward"));
 101   if (!AODforward->CheckEvent(fTriggerMask, fVtxMin, fVtxMax, 0,0,fTrigger)){
 102     AliError("Invalid Event");
 103     return;
 104   }
 105
 106   // Fill centrality histogram, only useful for PbPb
 107   Float_t cent = AODforward->GetCentrality();
 108   if(fLowCent<fHighCent){
 109     if(cent<fLowCent||cent>fHighCent) return;
 110   }
 111   fCent->Fill(cent);
 112
 113   TH2D forward;
 114   TH2D central;
 115
 116   // Get 2D eta-phi histograms for each event
 117   GetHistograms(aod, forward, central);
 118
 119   TH1D* dndetaSumForward    = (TH1D*)fOutput->FindObject("dndetaSumForward");
 120   TH1D* dndetaSumCentral    = (TH1D*)fOutput->FindObject("dndetaSumCentral");
 121   TH1D* dndetaEventForward  = (TH1D*)fOutput->FindObject("dndetaEventForward");
 122   TH1D* dndetaEventCentral  = (TH1D*)fOutput->FindObject("dndetaEventCentral");
 123   TH1D* normEventForward    = 0;
 124   TH1D* normEventCentral    = 0;
 125
 126   dndetaEventForward = forward.ProjectionX("dndetaForward",1,forward.GetNbinsY(),"");
 127   dndetaEventCentral = central.ProjectionX("dndetaCentral",1,central.GetNbinsY(),"");
 128   dndetaSumForward->Add(dndetaEventForward);
 129   dndetaSumCentral->Add(dndetaEventCentral);
 130
 131   // underflow eta bin of forward/central carry information on whether there is acceptance. 1= acceptance, 0= no acceptance
 132   normEventForward   = forward.ProjectionX("dndetaEventForwardNorm",0,0,"");
 133   normEventCentral   = central.ProjectionX("dndetaEventCentralNorm",0,0,"");
 134
 135   Double_t VtxZ = AODforward->GetIpZ();
 136
 137
 138   // loop over all eta bins, and fill multiplicity histos
 139   TIter next(&fBins);
 140   Bin * bin = 0;
 141   while ((bin = static_cast<Bin*>(next()))) {
 142     bin->Process(dndetaEventForward, dndetaEventCentral,
 143                  normEventForward,   normEventCentral, VtxZ);
 144   }
 145
 146   PostData(1, fOutput);
 147
 148 }
 149
 150 //_____________________________________________________________________
 151 void AliForwardMultiplicityDistribution::Terminate(Option_t */*option*/)
 152 {
 153   //
 154   // Terminate
 155   //
 156 }
 157
 158 //_________________________________________________________________-
 159 TH2D* AliForwardMultiplicityDistribution::GetHistogram(const AliAODEvent*, Bool_t)
 160 {
 161   //
 162   // implementation of pure virtual function, always returning 0
 163   //
 164   return 0;
 165 }
 166
 167 void AliForwardMultiplicityDistribution::GetHistograms(const AliAODEvent* aod, TH2D& forward, TH2D& central, Bool_t /*mc*/)
 168 {
 169   //
 170   // Get single event forward and central dN²/dEta dPhi histograms
 171   //
 172   TObject* forwardObj = 0;
 173   TObject* centralObj = 0;
 174
 175
 176   forwardObj = aod->FindListObject("Forward");
 177   centralObj = aod->FindListObject("CentralClusters");
 178
 179   AliAODForwardMult* AODforward = static_cast<AliAODForwardMult*>(forwardObj);
 180   AliAODCentralMult* AODcentral = static_cast<AliAODCentralMult*>(centralObj);
 181
 182   forward= AODforward->GetHistogram();
 183   central= AODcentral->GetHistogram();
 184
 185 }
 186
 187 //______________________________________________________________________
 188 #if 1
 189 const Char_t*
 190 #else
 191 TString
 192 #endif
 193 AliForwardMultiplicityDistribution::FormBinName(Double_t etaLow, Double_t etaHigh)
 194 {
 195   //
 196   //  Form name of eta bin
 197   //
 198   TString sLow(TString::Format("%+03d",int(10*etaLow)));
 199   TString sHigh(TString::Format("%+03d",int(10*etaHigh)));
 200   sLow.ReplaceAll("+", "plus");
 201   sLow.ReplaceAll("-", "minus");
 202   sHigh.ReplaceAll("+", "plus");
 203   sHigh.ReplaceAll("-", "minus");
 204  #if 0
 205   TString* combined= new TString();
 206   combined->Append(TString::Format("%s_%s", sLow.Data(), sHigh.Data()));
 207   return combined->Data();
 208 #else
 209   static TString tmp;
 210   tmp = TString::Format("%s_%s", sLow.Data(), sHigh.Data());
 211   return tmp.Data();
 212 #endif
 213 }
 214
 215 //=====================================================================
 216 AliForwardMultiplicityDistribution::Bin::Bin()
 217   : TNamed(),
 218     fEtaLow(-1111),     // low eta limit
 219     fEtaHigh(-1111),    // high eta limit
 220     fHist(0),           // multiplicity distribution hist
 221     fHistPlus05(0),     // multiplicity distribution hist scaled up with 5%
 222     fHistPlus075(0),    // multiplicity distribution hist scaled up with 7.5%
 223     fHistPlus10(0),     // multiplicity distribution hist scaled up with 10%
 224     fHistMinus05(0),    // multiplicity distribution hist scaled down with 5%
 225     fHistMinus075(0),   // multiplicity distribution hist scaled down with 7.5%
 226     fHistMinus10(0),    // multiplicity distribution hist scaled down with 10%
 227     fHistPlusSys(0),    // multiplicity distribution hist scaled up with the event uncertainty
 228     fHistMinusSys(0),   // multiplicity distribution hist scaled down with the event uncertainty
 229     fAcceptance(0),     // histogram showing the 'holes' in acceptance.
 230     fVtxZvsNdataBins(0) // VtxZ vs. number of data acceptance bins (normalised to the eta range)
 231 {
 232   //
 233   // Default Constructor
 234   //
 235 }
 236 //_____________________________________________________________________
 237 AliForwardMultiplicityDistribution::Bin::Bin(Double_t etaLow, Double_t etaHigh)
 238   : TNamed("", ""),
 239     fEtaLow(etaLow),    // low eta limit
 240     fEtaHigh(etaHigh),  // high eta limit
 241     fHist(0),           // multiplicity distribution hist
 242     fHistPlus05(0),     // multiplicity distribution hist scaled up with 5%
 243     fHistPlus075(0),    // multiplicity distribution hist scaled up with 7.5%
 244     fHistPlus10(0),     // multiplicity distribution hist scaled up with 10%
 245     fHistMinus05(0),    // multiplicity distribution hist scaled down with 5%
 246     fHistMinus075(0),   // multiplicity distribution hist scaled down with 7.5%
 247     fHistMinus10(0),    // multiplicity distribution hist scaled down with 10%
 248     fHistPlusSys(0),    // multiplicity distribution hist scaled up with the event uncertainty
 249     fHistMinusSys(0),   // multiplicity distribution hist scaled down with the event uncertainty
 250     fAcceptance(0),     // histogram showing the 'holes' in acceptance.
 251     fVtxZvsNdataBins(0) // VtxZ vs. number of data acceptance bins (normalised to the eta range)
 252 {
 253   //
 254   // Constructor
 255   //
 256   const char* name = AliForwardMultiplicityDistribution::FormBinName(fEtaLow,fEtaHigh);
 257
 258   SetName(name);
 259   SetTitle(TString::Format("%+4.1f < #eta < %+4.1f", fEtaLow, fEtaHigh));
 260
 261 }
 262 //_____________________________________________________________________
 263 void AliForwardMultiplicityDistribution::Bin::DefineOutputs(TList* cont,  Int_t max)
 264 {
 265   //
 266   // Define eta bin output histograms
 267   //
 268   TList* out = new TList;
 269   out->SetName(GetName());
 270   cont->Add(out);
 271
 272   fHist            = new TH1D("mult", GetTitle(), max, -0.5, max-.5);
 273   fHistPlus05      = new TH1D("multPlus05", GetTitle(), max, -0.5, max-.5);
 274   fHistPlus075     = new TH1D("multPlus075", GetTitle(), max, -0.5, max-.5);
 275   fHistPlus10      = new TH1D("multPlus10", GetTitle(), max, -0.5, max-.5);
 276   fHistMinus05     = new TH1D("multMinus05", GetTitle(), max, -0.5, max-.5);
 277   fHistMinus075    = new TH1D("multMinus075", GetTitle(), max, -0.5, max-.5);
 278   fHistMinus10     = new TH1D("multMinus10", GetTitle(), max, -0.5, max-.5);
 279   fHistPlusSys     = new TH1D("multPlusSys", GetTitle(), max, -0.5, max-.5);
 280   fHistMinusSys    = new TH1D("multMinusSys", GetTitle(), max, -0.5, max-.5);
 281   fVtxZvsNdataBins = new TH2D("VtxZvsNdataBins", "VtxZ vs dataAcceptance/etaRange;z-vtz;dataAcceptance/etaRange", 20, -10,10, 130,0,1.3);
 282   fAcceptance      = new TH2D("Acceptance","Acceptance;#eta;z-vtx",200,-4, 6 , 20,-10,10);
 283
 284   out->Add(fHist);
 285   out->Add(fHistPlus05);
 286   out->Add(fHistPlus075);
 287   out->Add(fHistPlus10);
 288   out->Add(fHistMinus05);
 289   out->Add(fHistMinus075);
 290   out->Add(fHistMinus10);
 291   out->Add(fHistPlusSys);
 292   out->Add(fHistMinusSys);
 293   out->Add(fVtxZvsNdataBins);
 294   out->Add(fAcceptance);
 295
 296 }
 297
 298 //_____________________________________________________________________
 299 void AliForwardMultiplicityDistribution::Bin::Process(TH1D* dndetaForward, TH1D* dndetaCentral,
 300                                                       TH1D* normForward,   TH1D* normCentral, Double_t VtxZ)
 301 {
 302   //
 303   // Process single eta bin
 304   //
 305
 306   // Diagnostics on event acceptance
 307   Int_t    first = dndetaForward->GetXaxis()->FindBin(fEtaLow);
 308   Int_t    last  = dndetaForward->GetXaxis()->FindBin(fEtaHigh-.0001);
 309   Double_t acceptanceBins=0;
 310   for(Int_t n= first;n<=last;n++){
 311     if(normForward->GetBinContent(n)>0||normCentral->GetBinContent(n)>0){
 312       acceptanceBins++;
 313     }
 314     fAcceptance->SetBinContent(n, fAcceptance->GetYaxis()->FindBin(VtxZ), 1);
 315     if(normForward->GetBinContent(n)>0||normCentral->GetBinContent(n)>0)
 316       fAcceptance->SetBinContent(n, fAcceptance->GetYaxis()->FindBin(VtxZ),10);
 317   }
 318
 319   //std::cout << "% of bins covered in eta-range: " << acceptanceBins/(last-first+1) << std::endl;
 320   fVtxZvsNdataBins->Fill(VtxZ, (Double_t)acceptanceBins/(last-first+1));
 321
 322
 323   Double_t c     = 0;
 324   Double_t e2    = 0;
 325   for (Int_t i = first; i <= last; i++) {
 326     Double_t cForward  = 0;
 327     Double_t cCentral  = 0;
 328     Double_t e2Forward = 0;
 329     Double_t e2Central = 0;
 330     if (normForward->GetBinContent(i) > 0) {
 331       cForward  = dndetaForward->GetBinContent(i);
 332       e2Forward = dndetaForward->GetBinError(i) * dndetaForward->GetBinError(i);
 333     }
 334     if (normCentral->GetBinContent(i) > 0) {
 335       cCentral  = dndetaCentral->GetBinContent(i);
 336       e2Central = dndetaCentral->GetBinError(i) * dndetaCentral->GetBinError(i);
 337     }
 338     Double_t cc  = 0;
 339     Double_t ee2 = 0;
 340     // if overlap between central/forward, use average of the two
 341     if (cCentral > 0 && cForward > 0) {
 342       cc  = 0.5 * (cForward  + cCentral);
 343       ee2 = 0.5 * (e2Forward + e2Central);
 344     }
 345     else if (cCentral > 0) {
 346       cc  = cCentral;
 347       ee2 = e2Central;
 348     }
 349     else if (cForward > 0) {
 350       cc  = cForward;
 351       ee2 = e2Forward;
 352     }
 353     c  += cc;
 354     e2 += ee2;
 355   }
 356
 357   //Systematic errors from here
 358
 359   Double_t fmd=0;
 360   Double_t spd=0;
 361   Double_t overlap=0;
 362
 363   // number of eta bins in fmd, spd and overlap respectively
 364   for(Int_t i = first;i<=last;i++){
 365     if(normForward->GetBinContent(i)>0&&normCentral->GetBinContent(i)<1)
 366       fmd++;
 367     if(normForward->GetBinContent(i)>0&&normCentral->GetBinContent(i)>0)
 368       overlap++;
 369     if(normCentral->GetBinContent(i)>0&&normForward->GetBinContent(i)<1)
 370       spd++;
 371   }
 372
 373   Double_t sysErrorSquared = 0;
 374   // estimate of systematic uncertainty on the event multiplicity - hardcoded :(. estimates taken from Hans Hjersing Dalsgaard or Casper Nygaard phd theses.
 375   Double_t fmdSysError= 0.08;
 376   Double_t spdSysError= 0.067;
 377
 378   Double_t total = 0;
 379   total= fmd+ spd + overlap;
 380   // Combined systematc event uncertainty, by weighting with the number of eta-bins of fmd-only, spd-only and the overlap.
 381   if(total){
 382     sysErrorSquared= (fmd*TMath::Power(fmdSysError,2)+ spd*TMath::Power(spdSysError,2)+
 383                       0.5*overlap*TMath::Power(fmdSysError,2)+ 0.5*overlap*TMath::Power(spdSysError,2))/total;
 384   }
 385
 386   //Strangeness correction, assumed to be 2.0 +- 1% (syst),  taken from HHD and CN thesis
 387   c=0.98*c;
 388
 389
 390   // correction for missing material description. Correction is based on separate PbPb analysis of difference between nominel and displaced vertices
 391   if(fEtaHigh<1.55&&fEtaHigh>1.45)
 392     c=0.98*c;
 393   if(fEtaHigh<2.05&&fEtaHigh>1.95)
 394     c=0.963*c;
 395   if(fEtaHigh<2.45&&fEtaHigh>2.35)
 396     c=0.956*c;
 397   if(fEtaHigh>2.95)
 398     c=0.955*c;
 399
 400   fHist->Fill(c);
 401   fHistPlusSys->Fill((1+TMath::Sqrt(sysErrorSquared))*c);
 402   fHistMinusSys->Fill((1-TMath::Sqrt(sysErrorSquared))*c);
 403   fHistPlus05->Fill(1.05*c);
 404   fHistPlus075->Fill(1.075*c);
 405   fHistPlus10->Fill(1.1*c);
 406   fHistMinus05->Fill(0.95*c);
 407   fHistMinus075->Fill(0.925*c);
 408   fHistMinus10->Fill(0.9*c);
 409
 410
 411 }
 412
 413
 414
 415
 416 //_____________________________________________________________________
 417 //
 418 //
 419 // EOF