PWG2/FORWARD/analysis2/AliBasedNdetaTask.cxx

   1 //====================================================================
   2 #include "AliBasedNdetaTask.h"
   3 #include <TMath.h>
   4 #include <TH2D.h>
   5 #include <TH1D.h>
   6 #include <THStack.h>
   7 #include <TList.h>
   8 #include <TParameter.h>
   9 #include <AliAnalysisManager.h>
  10 #include <AliAODEvent.h>
  11 #include <AliAODHandler.h>
  12 #include <AliAODInputHandler.h>
  13 #include "AliForwardUtil.h"
  14 #include "AliAODForwardMult.h"
  15 #include <TFile.h>
  16
  17 //____________________________________________________________________
  18 AliBasedNdetaTask::AliBasedNdetaTask()
  19   : AliAnalysisTaskSE(),
  20     fSum(0),            //  Sum of histograms
  21     fSumMC(0),          //  Sum of MC histograms (if any)
  22     fSums(0),           // Container of sums
  23     fOutput(0),         // Container of outputs
  24     fTriggers(0),       // Histogram of triggers
  25     fVtxMin(0),         // Minimum v_z
  26     fVtxMax(0),         // Maximum v_z
  27     fTriggerMask(0),    // Trigger mask
  28     fRebin(0),          // Rebinning factor
  29     fCutEdges(false),
  30     fSymmetrice(true),
  31     fCorrEmpty(true),
  32     fTriggerEff(1),
  33     fShapeCorr(0),
  34     fCentLow(0),
  35     fCentHigh(100)
  36 {
  37   //
  38   // Constructor
  39   //
  40 }
  41
  42 //____________________________________________________________________
  43 AliBasedNdetaTask::AliBasedNdetaTask(const char* name)
  44   : AliAnalysisTaskSE(name),
  45     fSum(0),            // Sum of histograms
  46     fSumMC(0),          // Sum of MC histograms (if any)
  47     fSums(0),           // Container of sums
  48     fOutput(0),         // Container of outputs
  49     fTriggers(0),       // Histogram of triggers
  50     fVtxMin(-10),       // Minimum v_z
  51     fVtxMax(10),        // Maximum v_z
  52     fTriggerMask(AliAODForwardMult::kInel),
  53     fRebin(5),          // Rebinning factor
  54     fCutEdges(false),
  55     fSymmetrice(true),
  56     fCorrEmpty(true),
  57     fTriggerEff(1),
  58     fShapeCorr(0),
  59     fCentLow(0),
  60     fCentHigh(100)
  61 {
  62   //
  63   // Constructor
  64   //
  65
  66   // Output slot #1 writes into a TH1 container
  67   DefineOutput(1, TList::Class());
  68   DefineOutput(2, TList::Class());
  69 }
  70
  71 //____________________________________________________________________
  72 AliBasedNdetaTask::AliBasedNdetaTask(const AliBasedNdetaTask& o)
  73   : AliAnalysisTaskSE(o),
  74     fSum(o.fSum),       // TH2D* -  Sum of histograms
  75     fSumMC(o.fSumMC),// TH2D* -  Sum of MC histograms (if any)
  76     fSums(o.fSums),             // TList* - Container of sums
  77     fOutput(o.fOutput),         // TList* - Container of outputs
  78     fTriggers(o.fTriggers),     // TH1D* - Histogram of triggers
  79     fVtxMin(o.fVtxMin),         // Double_t - Minimum v_z
  80     fVtxMax(o.fVtxMax),         // Double_t - Maximum v_z
  81     fTriggerMask(o.fTriggerMask),// Int_t - Trigger mask
  82     fRebin(o.fRebin),           // Int_t - Rebinning factor
  83     fCutEdges(o.fCutEdges),     // Bool_t - Whether to cut edges when rebinning
  84     fSymmetrice(true),
  85     fCorrEmpty(true),
  86     fTriggerEff(o.fTriggerEff),
  87     fShapeCorr(o.fShapeCorr),
  88     fCentLow(o.fCentLow),
  89     fCentHigh(o.fCentHigh)
  90 {}
  91
  92 //____________________________________________________________________
  93 AliBasedNdetaTask::~AliBasedNdetaTask()
  94 {
  95   //
  96   // Destructor
  97   //
  98   if (fSums) {
  99     fSums->Delete();
 100     delete fSums;
 101     fSums = 0;
 102   }
 103   if (fOutputs) {
 104     fOutputs->Delete();
 105     delete fOutputs;
 106     fOutputs = 0;
 107   }
 108 }
 109
 110 //________________________________________________________________________
 111 void
 112 AliBasedNdetaTask::SetTriggerMask(const char* mask)
 113 {
 114   //
 115   // Set the trigger maskl
 116   //
 117   // Parameters:
 118   //    mask Trigger mask
 119   //
 120   UShort_t    trgMask = 0;
 121   TString     trgs(mask);
 122   trgs.ToUpper();
 123   TObjString* trg;
 124   TIter       next(trgs.Tokenize(" ,|"));
 125   while ((trg = static_cast<TObjString*>(next()))) {
 126     TString s(trg->GetString());
 127     if      (s.IsNull()) continue;
 128     if      (s.CompareTo("INEL")  == 0) trgMask = AliAODForwardMult::kInel;
 129     else if (s.CompareTo("INEL>0")== 0) trgMask = AliAODForwardMult::kInelGt0;
 130     else if (s.CompareTo("NSD")   == 0) trgMask = AliAODForwardMult::kNSD;
 131     else
 132       Warning("SetTriggerMask", "Unknown trigger %s", s.Data());
 133   }
 134   if (trgMask == 0) trgMask = 1;
 135   SetTriggerMask(trgMask);
 136 }
 137
 138 //________________________________________________________________________
 139 void
 140 AliBasedNdetaTask::SetShapeCorrection(const TH1* c)
 141 {
 142   //
 143   // Set the shape correction (a.k.a., track correction) for selected
 144   // trigger(s)
 145   //
 146   // Parameters:
 147   //    h Correction
 148   //
 149   if (!c) return;
 150   fShapeCorr = static_cast<TH1*>(c->Clone());
 151   fShapeCorr->SetDirectory(0);
 152 }
 153
 154 //________________________________________________________________________
 155 void
 156 AliBasedNdetaTask::UserCreateOutputObjects()
 157 {
 158   //
 159   // Create output objects.
 160   //
 161   // This is called once per slave process
 162   //
 163   fOutput = new TList;
 164   fOutput->SetName(Form("%s_result", GetName()));
 165   fOutput->SetOwner();
 166
 167   fSums = new TList;
 168   fSums->SetName(Form("%s_sums", GetName()));
 169   fSums->SetOwner();
 170
 171
 172   fTriggers = new TH1D("triggers", "Number of triggers",
 173                        kAccepted, 1, kAccepted);
 174   fTriggers->SetYTitle("# of events");
 175   fTriggers->GetXaxis()->SetBinLabel(kAll,         "All events");
 176   fTriggers->GetXaxis()->SetBinLabel(kB,           "w/B trigger");
 177   fTriggers->GetXaxis()->SetBinLabel(kA,           "w/A trigger");
 178   fTriggers->GetXaxis()->SetBinLabel(kC,           "w/C trigger");
 179   fTriggers->GetXaxis()->SetBinLabel(kE,           "w/E trigger");
 180   fTriggers->GetXaxis()->SetBinLabel(kMB,          "w/Collision trigger");
 181   fTriggers->GetXaxis()->SetBinLabel(kWithVertex,  "w/Vertex");
 182   fTriggers->GetXaxis()->SetBinLabel(kWithTrigger, "w/Selected trigger");
 183   fTriggers->GetXaxis()->SetBinLabel(kPileUp,      "with pileup");
 184   fTriggers->GetXaxis()->SetBinLabel(kAccepted,    "Accepted by cut");
 185   fTriggers->GetXaxis()->SetNdivisions(kAccepted, false);
 186   fTriggers->SetFillColor(kRed+1);
 187   fTriggers->SetFillStyle(3001);
 188   fTriggers->SetStats(0);
 189   fSums->Add(fTriggers);
 190
 191   // Check that we have an AOD input handler
 192   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 193   AliAODInputHandler* ah =
 194     dynamic_cast<AliAODInputHandler*>(am->GetInputEventHandler());
 195   if (!ah) AliFatal("No AOD input handler set in analysis manager");
 196
 197   // Post data for ALL output slots >0 here, to get at least an empty histogram
 198   PostData(1, fSums);
 199 }
 200
 201 //____________________________________________________________________
 202 TH2D*
 203 AliBasedNdetaTask::CloneHist(const TH2D* in, const char* name)
 204 {
 205   //
 206   // Clone a 2D histogram
 207   //
 208   // Parameters:
 209   //    in    Histogram to clone.
 210   //    name  New name of clone.
 211   //
 212   // Return:
 213   //    The clone
 214   //
 215   if (!in) return 0;
 216   TH2D* ret = static_cast<TH2D*>(in->Clone(name));
 217   ret->SetDirectory(0);
 218   ret->Sumw2();
 219   fSums->Add(ret);
 220
 221   return ret;
 222 }
 223
 224 //____________________________________________________________________
 225 Bool_t
 226 AliBasedNdetaTask::CheckEvent(const AliAODEvent* aod)
 227 {
 228   //
 229   // Check the trigger and vertex
 230   //
 231   // Parameters:
 232   //    aod
 233   //
 234   // Return:
 235   //
 236   //
 237   TObject*           oForward   = aod->FindListObject("Forward");
 238   if (!oForward) {
 239     AliWarning("No forward object found");
 240     return false;
 241   }
 242   AliAODForwardMult* forward   = static_cast<AliAODForwardMult*>(oForward);
 243
 244   // Count event
 245   fTriggers->AddBinContent(kAll);
 246   if (forward->IsTriggerBits(AliAODForwardMult::kB))
 247     fTriggers->AddBinContent(kB);
 248   if (forward->IsTriggerBits(AliAODForwardMult::kA))
 249     fTriggers->AddBinContent(kA);
 250   if (forward->IsTriggerBits(AliAODForwardMult::kC))
 251     fTriggers->AddBinContent(kC);
 252   if (forward->IsTriggerBits(AliAODForwardMult::kE))
 253     fTriggers->AddBinContent(kE);
 254   if (forward->IsTriggerBits(AliAODForwardMult::kInel))
 255     fTriggers->AddBinContent(kMB);
 256   if (forward->IsTriggerBits(AliAODForwardMult::kMCNSD))
 257     fTriggers->AddBinContent(kMCNSD);
 258   if (forward->IsTriggerBits(AliAODForwardMult::kPileUp))
 259     fTriggers->AddBinContent(kPileUp);
 260
 261
 262
 263   // Check if we have an event of interest.
 264   if (!forward->IsTriggerBits(fTriggerMask)) return false;
 265
 266   //Check for pileup
 267   if (forward->IsTriggerBits(AliAODForwardMult::kPileUp)) return false;
 268
 269   fTriggers->AddBinContent(kWithTrigger);
 270
 271   // Check that we have a valid vertex
 272   if (!forward->HasIpZ()) return false;
 273   fTriggers->AddBinContent(kWithVertex);
 274
 275   // Check that vertex is within cuts
 276   if (!forward->InRange(fVtxMin, fVtxMax)) return false;
 277   fTriggers->AddBinContent(kAccepted);
 278
 279   return true;
 280 }
 281 //____________________________________________________________________
 282 void
 283 AliBasedNdetaTask::UserExec(Option_t *)
 284 {
 285   //
 286   // Process a single event
 287   //
 288   // Parameters:
 289   //    option Not used
 290   //
 291   // Main loop
 292   AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
 293   if (!aod) {
 294     AliError("Cannot get the AOD event");
 295     return;
 296   }
 297
 298   TObject* obj = 0;
 299   obj = aod->FindListObject("Forward");
 300
 301   AliAODForwardMult* forward = static_cast<AliAODForwardMult*>(obj);
 302   Float_t cent = forward->GetCentrality();
 303
 304   if(cent > 0) {
 305     if( cent < fCentLow || cent > fCentHigh ) return;
 306     else std::cout<<"selecting event with cent "<<cent<<std::endl;
 307   }
 308   //std::cout<<fCentLow<<"  "<<fCentHigh<<std::endl;
 309   // Get the histogram(s)
 310   TH2D* data   = GetHistogram(aod, false);
 311   TH2D* dataMC = GetHistogram(aod, true);
 312
 313   // We should have a base object at least
 314   if (!data) {
 315     AliWarning("No data object found in AOD");
 316     return;
 317   }
 318
 319
 320
 321   // Check the event
 322   if (!CheckEvent(aod)) return;
 323
 324   // Create our sum histograms
 325   if (!fSum)  {
 326     fSum   = CloneHist(data,  GetName());
 327     if(forward->GetSystem()  < 2)
 328       LoadNormalizationData(forward->GetSystem(),forward->GetSNN());
 329   }
 330   else fSum->Add((data));
 331
 332   //for(Int_t i = 1; i<=data->GetNbinsX(); i++) {
 333   //  std::cout<<fSum->GetBinContent(i,0) <<"   "<<data->GetBinContent(i,0)<<std::endl;
 334   // }
 335
 336   if (!fSumMC && dataMC) fSumMC = CloneHist(dataMC,Form("%sMC", GetName()));
 337   else if (fSumMC) fSumMC->Add((dataMC));
 338   // Add contribution
 339
 340   PostData(1, fSums);
 341 }
 342
 343 //________________________________________________________________________
 344 void
 345 AliBasedNdetaTask::SetHistogramAttributes(TH1D* h, Int_t colour, Int_t marker,
 346                                           const char* title, const char* ytitle)
 347 {
 348   //
 349   // Set histogram graphical options, etc.
 350   //
 351   // Parameters:
 352   //    h       Histogram to modify
 353   //    colour  Marker color
 354   //    marker  Marker style
 355   //    title   Title of histogram
 356   //    ytitle  Title on y-axis.
 357   //
 358   h->SetTitle(title);
 359   h->SetMarkerColor(colour);
 360   h->SetMarkerStyle(marker);
 361   h->SetMarkerSize(1);
 362   h->SetFillStyle(0);
 363   h->SetYTitle(ytitle);
 364   h->GetXaxis()->SetTitleFont(132);
 365   h->GetXaxis()->SetLabelFont(132);
 366   h->GetXaxis()->SetNdivisions(10);
 367   h->GetYaxis()->SetTitleFont(132);
 368   h->GetYaxis()->SetLabelFont(132);
 369   h->GetYaxis()->SetNdivisions(10);
 370   h->GetYaxis()->SetDecimals();
 371   h->SetStats(0);
 372 }
 373
 374 //________________________________________________________________________
 375 TH1D*
 376 AliBasedNdetaTask::ProjectX(const TH2D* h,
 377                             const char* name,
 378                             Int_t firstbin,
 379                             Int_t lastbin,
 380                             bool  corr,
 381                             bool  error) const
 382 {
 383   //
 384   // Project onto the X axis
 385   //
 386   // Parameters:
 387   //    h         2D histogram
 388   //    name      New name
 389   //    firstbin  First bin to use
 390   //    lastbin   Last bin to use
 391   //    error     Whether to calculate errors
 392   //
 393   // Return:
 394   //    Newly created histogram or null
 395   //
 396   if (!h) return 0;
 397   //#if USE_ROOT_PROJECT
 398   return h->ProjectionX(name, firstbin, lastbin, (error ? "e" : ""));
 399   //#endif
 400
 401   TAxis* xaxis = h->GetXaxis();
 402   TAxis* yaxis = h->GetYaxis();
 403   TH1D*  ret   = new TH1D(name, h->GetTitle(), xaxis->GetNbins(),
 404                           xaxis->GetXmin(), xaxis->GetXmax());
 405   static_cast<const TAttLine*>(h)->Copy(*ret);
 406   static_cast<const TAttFill*>(h)->Copy(*ret);
 407   static_cast<const TAttMarker*>(h)->Copy(*ret);
 408   ret->GetXaxis()->ImportAttributes(xaxis);
 409
 410   Int_t first = firstbin;
 411   Int_t last  = lastbin;
 412   if (first < 0)                         first = 0;
 413   else if (first >= yaxis->GetNbins()+1) first = yaxis->GetNbins();
 414   if      (last  < 0)                    last  = yaxis->GetNbins();
 415   else if (last  >  yaxis->GetNbins()+1) last  = yaxis->GetNbins();
 416   if (last-first < 0) {
 417     AliWarning(Form("Nothing to project [%d,%d]", first, last));
 418     return 0;
 419
 420   }
 421
 422   // Loop over X bins
 423   // AliInfo(Form("Projecting bins [%d,%d] of %s", first, last, h->GetName()));
 424   Int_t ybins = (last-first+1);
 425   for (Int_t xbin = 0; xbin <= xaxis->GetNbins()+1; xbin++) {
 426     Double_t content = 0;
 427     Double_t error2  = 0;
 428     Int_t    nbins   = 0;
 429
 430
 431     for (Int_t ybin = first; ybin <= last; ybin++) {
 432       Double_t c1 = h->GetCellContent(xbin, ybin);
 433       Double_t e1 = h->GetCellError(xbin, ybin);
 434
 435       // Ignore empty bins
 436       if (c1 < 1e-12) continue;
 437       if (e1 < 1e-12) {
 438         if (error) continue;
 439         e1 = 1;
 440       }
 441
 442       content    += c1;
 443       error2     += e1*e1;
 444       nbins++;
 445     } // for (ybin)
 446     if(content > 0 && nbins > 0) {
 447       Double_t factor = (corr ? Double_t(ybins) / nbins : 1);
 448       if (error) {
 449         // Calculate weighted average
 450         ret->SetBinContent(xbin, content * factor);
 451         ret->SetBinError(xbin, factor * TMath::Sqrt(error2));
 452       }
 453       else
 454         ret->SetBinContent(xbin, factor * content);
 455     }
 456   } // for (xbin)
 457
 458   return ret;
 459 }
 460
 461 //________________________________________________________________________
 462 void
 463 AliBasedNdetaTask::Terminate(Option_t *)
 464 {
 465   //
 466   // Called at end of event processing..
 467   //
 468   // This is called once in the master
 469   //
 470   // Parameters:
 471   //    option Not used
 472   //
 473   // Draw result to screen, or perform fitting, normalizations Called
 474   // once at the end of the query
 475
 476   fSums = dynamic_cast<TList*> (GetOutputData(1));
 477   if(!fSums) {
 478     AliError("Could not retrieve TList fSums");
 479     return;
 480   }
 481
 482   if (!fOutput) {
 483     fOutput = new TList;
 484     fOutput->SetName(Form("%s_result", GetName()));
 485     fOutput->SetOwner();
 486   }
 487
 488   fSum     = static_cast<TH2D*>(fSums->FindObject(GetName()));
 489   fSumMC   = static_cast<TH2D*>(fSums->FindObject(Form("%sMC", GetName())));
 490   fTriggers       = static_cast<TH1D*>(fSums->FindObject("triggers"));
 491
 492   if (!fTriggers) {
 493     AliError("Couldn't find histogram 'triggers' in list");
 494     return;
 495   }
 496   if (!fSum) {
 497     AliError("Couldn't find histogram 'base' in list");
 498     return;
 499   }
 500
 501   Int_t nAll        = Int_t(fTriggers->GetBinContent(kAll));
 502   Int_t nB          = Int_t(fTriggers->GetBinContent(kB));
 503   Int_t nA          = Int_t(fTriggers->GetBinContent(kA));
 504   Int_t nC          = Int_t(fTriggers->GetBinContent(kC));
 505   Int_t nE          = Int_t(fTriggers->GetBinContent(kE));
 506   Int_t nMB         = Int_t(fTriggers->GetBinContent(kMB));
 507   Int_t nTriggered  = Int_t(fTriggers->GetBinContent(kWithTrigger));
 508   Int_t nWithVertex = Int_t(fTriggers->GetBinContent(kWithVertex));
 509   Int_t nAccepted   = Int_t(fTriggers->GetBinContent(kAccepted));
 510   Int_t nGood       = nB - nA - nC + 2 * nE;
 511
 512   //  Int_t nBg         = nA + nC -nE;
 513   //std::cout<<"nBackground "<<nBg<<"   , with vertex in range "<<fNbgValid<<" ,  total "<<fNbgVertex<<"  "<<std::endl<<std::endl;
 514   Float_t alpha            = ((Float_t)nAccepted) / ((Float_t)nWithVertex);
 515   //Float_t alphaBG  = 1;
 516   //if(fNbgVertex > 0) alphaBG= (Float_t)fNbgValid / (Float_t)fNbgVertex;
 517   //Float_t nNormalizationJF = nAccepted + alpha*(nMB - nAccepted - nBg);
 518   Float_t nNormalizationJF = nAccepted + alpha*(nTriggered - nWithVertex) ; // -alphaBG*(nBg-fNbgVertex);
 519
 520   //Float_t nNormalizationBg = fNbgValid + alphaBG*nBg;
 521
 522   //std::cout<<nTriggered -nWithVertex<<std::endl;
 523
 524   //Double_t vtxEff   = Double_t(nMB) / nTriggered * Double_t(nAccepted) / nGood;
 525   Double_t vtxEff   =  (Float_t)nAccepted / nNormalizationJF;
 526   vtxEff = vtxEff*fTriggerEff;
 527
 528   //if ( fTriggerMask == AliAODForwardMult::kNSD ) vtxEff = 0.97; //From paper...
 529   //std::cout<<
 530
 531   //Double_t vtxEffBg   =  (Float_t)fNbgValid / nNormalizationBg;
 532   //Double_t vtxEff   =  1.;//(Float_t)nAccepted / nNormalizationJF;
 533
 534
 535
 536
 537   if (nTriggered <= 0) {
 538     AliError("Number of triggered events <= 0");
 539     return;
 540   }
 541   if (nGood <= 0) {
 542     AliWarning(Form("Number of good events=%d=%d-%d-%d+2*%d<=0",
 543                     nGood, nB, nA, nC, nE));
 544     nGood = nMB;
 545   }
 546   //Double_t vtxEff   =  Double_t(nMB) / nTriggered * Double_t(nAccepted) / nGood;
 547
 548   Double_t vNorm    = Double_t(nAccepted) / nGood;
 549   AliInfo(Form("Total of %9d events\n"
 550                "                   of these %9d are minimum bias\n"
 551                "                   of these %9d has a %s trigger\n"
 552                "                   of these %9d has a vertex\n"
 553                "                   of these %9d were in [%+4.1f,%+4.1f]cm\n"
 554                "                   Triggers by type:\n"
 555                "                     B   = %9d\n"
 556                "                     A|C = %9d (%9d+%-9d)\n"
 557                "                     E   = %9d\n"
 558                "                   Implies %9d good triggers\n"
 559                "                   Vertex efficiency: %f (%f)",
 560                nAll, nMB, nTriggered,
 561                AliAODForwardMult::GetTriggerString(fTriggerMask),
 562                nWithVertex, nAccepted,
 563                fVtxMin, fVtxMax,
 564                nB, nA+nC, nA, nC, nE, nGood, vtxEff, vNorm));
 565
 566   const char* name = GetName();
 567
 568   // Get acceptance normalisation from underflow bins
 569   TH1D* norm   = ProjectX(fSum, Form("norm%s",name), 0, 0, fCorrEmpty, false);
 570
 571
 572
 573
 574   //
 575   //std::cout<<norm->GetMaximumBin()<<"    "<< (Float_t)nAccepted / norm->GetBinContent((Float_t)norm->GetMaximumBin()) <<std::endl;
 576   Float_t scaleForNorm =  (Float_t)nAccepted / (Float_t)norm->GetBinContent(norm->GetMaximumBin()) ;
 577   //Float_t scaleForNorm =  norm->Integral() ;
 578   std::cout<<norm->Integral()<<std::endl;
 579
 580   norm->Scale(scaleForNorm);
 581   //norm->Scale(1., "width");
 582   //norm->Scale(norm->GetNbinsX() / (norm->GetXaxis()->GetXmax() - norm->GetXaxis()->GetXmin() ));
 583
 584   //norm->Scale(1.,"width");
 585   //
 586   // Project onto eta axis - _ignoring_underflow_bins_!
 587
 588   TH2D* tmpNorm = (TH2D*)fSum->Clone("tmpNorm");
 589
 590   if(fShapeCorr)
 591     fSum->Divide(fShapeCorr);
 592
 593   TH1D* dndeta = ProjectX(fSum, Form("dndeta%s",name), 1, fSum->GetNbinsY(),
 594                           fCorrEmpty);
 595
 596   // Normalize to the acceptance
 597   //dndeta->Divide(norm);
 598
 599   for(Int_t i = 1; i<=tmpNorm->GetNbinsX(); i++) {
 600
 601     if( tmpNorm->GetBinContent(i,0) > 0 ) {
 602       //std::cout<<tmpNorm->GetBinContent(i,0) <<std::endl;
 603       dndeta->SetBinContent(i,dndeta->GetBinContent(i) / tmpNorm->GetBinContent(i,0));
 604       dndeta->SetBinError(i,dndeta->GetBinError(i) / tmpNorm->GetBinContent(i,0));
 605     }
 606   }
 607
 608   // Scale by the vertex efficiency
 609   //dndeta->Scale(vNorm, "width");
 610
 611   dndeta->Scale(vtxEff, "width");
 612
 613   //dndeta->Scale(1./(Float_t)nAccepted);
 614   //dndeta->Scale(dndeta->GetNbinsX() / (dndeta->GetXaxis()->GetXmax() - dndeta->GetXaxis()->GetXmin() ));
 615
 616   //norm->Scale(1. / nAccepted);
 617   //norm->Scale(1. / nNormalizationJF);
 618   SetHistogramAttributes(dndeta, kRed+1, 20, Form("ALICE %s", name));
 619   SetHistogramAttributes(norm, kRed+1,20,Form("ALICE %s normalisation", name));
 620
 621   fOutput->Add(fTriggers->Clone());
 622   if (fSymmetrice)   fOutput->Add(Symmetrice(dndeta));
 623   fOutput->Add(dndeta);
 624   fOutput->Add(norm);
 625   fOutput->Add(Rebin(dndeta));
 626   if (fSymmetrice)   fOutput->Add(Symmetrice(Rebin(dndeta)));
 627
 628   if (fSumMC) {
 629     // Get acceptance normalisation from underflow bins
 630     TH1D* normMC   = ProjectX(fSumMC,Form("norm%sMC", name), 0, 1,
 631                               fCorrEmpty, false);
 632     // Project onto eta axis - _ignoring_underflow_bins_!
 633     TH1D* dndetaMC = ProjectX(fSumMC,Form("dndeta%sMC", name),1,
 634                               fSum->GetNbinsY(), fCorrEmpty);
 635     // Normalize to the acceptance
 636     dndetaMC->Divide(normMC);
 637     // Scale by the vertex efficiency
 638     dndetaMC->Scale(vNorm, "width");
 639     normMC->Scale(1. / nAccepted);
 640
 641     SetHistogramAttributes(dndetaMC, kRed+3, 21, Form("ALICE %s (MC)",name));
 642     SetHistogramAttributes(normMC, kRed+3, 21, Form("ALICE %s (MC) "
 643                                                     "normalisation",name));
 644
 645     fOutput->Add(dndetaMC);
 646     if (fSymmetrice)   fOutput->Add(Symmetrice(dndetaMC));
 647     fOutput->Add(normMC);
 648     fOutput->Add(Rebin(dndetaMC));
 649     if (fSymmetrice)   fOutput->Add(Symmetrice(Rebin(dndetaMC)));
 650   }
 651
 652   TNamed* trigString =
 653     new TNamed("trigString", AliAODForwardMult::GetTriggerString(fTriggerMask));
 654   trigString->SetUniqueID(fTriggerMask);
 655   fOutput->Add(trigString);
 656
 657   TAxis* vtxAxis = new TAxis(1,fVtxMin,fVtxMax);
 658   vtxAxis->SetName("vtxAxis");
 659   vtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", fVtxMin,fVtxMax));
 660   fOutput->Add(vtxAxis);
 661   fOutput->Add(new TParameter<Double_t>("triggerEff", fTriggerEff));
 662   if (fShapeCorr) fOutput->Add(fShapeCorr);
 663
 664   PostData(2, fOutput);
 665 }
 666 //________________________________________________________________________
 667 void
 668 AliBasedNdetaTask::LoadNormalizationData(UShort_t sys, UShort_t energy)
 669 {
 670   sys = 1;
 671   energy = 900;
 672   TString type("pp");
 673   if(sys == 2) type.Form("PbPb");
 674   TString snn("900");
 675   if(energy == 7000) snn.Form("7000");
 676   if(energy == 2750) snn.Form("2750");
 677
 678   if(fShapeCorr && (fTriggerEff != 1)) {AliInfo("Objects already set for normalization - no action taken"); return; }
 679
 680   TFile* fin = TFile::Open(Form("$ALICE_ROOT/PWG2/FORWARD/corrections/Normalization/normalizationHists_%s_%s.root",type.Data(),snn.Data()));
 681   if(!fin) AliWarning("no file for normalization");
 682
 683   if ( fTriggerMask == AliAODForwardMult::kInel ) {
 684     TH2F* shapeCor = dynamic_cast<TH2F*>(fin->Get("hInelNormalization"));
 685     if(shapeCor) SetShapeCorrection(shapeCor);
 686     TParameter<float>* ineleff = (TParameter<float>*)fin->Get("inelTriggerEff");
 687   if ( ineleff) SetTriggerEff(ineleff->GetVal());
 688   if (shapeCor && ineleff) AliInfo("Loaded objects for normalization.");
 689
 690 }
 691
 692 if ( fTriggerMask == AliAODForwardMult::kNSD ) {
 693   TH2F* shapeCor = dynamic_cast<TH2F*>(fin->Get("hNSDNormalization"));
 694   if(shapeCor) SetShapeCorrection(shapeCor);
 695   TParameter<float>* nsdeff = (TParameter<float>*)fin->Get("nsdTriggerEff");
 696 if ( nsdeff) SetTriggerEff(nsdeff->GetVal());
 697 if (shapeCor && nsdeff) AliInfo("Loaded objects for normalization.");
 698 }
 699
 700 if ( fTriggerMask == AliAODForwardMult::kInelGt0 ) {
 701   TH2F* shapeCor = dynamic_cast<TH2F*>(fin->Get("hInelGt0Normalization"));
 702     if(shapeCor) SetShapeCorrection(shapeCor);
 703     TParameter<float>* inelgt0eff = (TParameter<float>*)fin->Get("inelgt0TriggerEff");
 704 if ( inelgt0eff) SetTriggerEff(inelgt0eff->GetVal());
 705 if (shapeCor && inelgt0eff) AliInfo("Loaded objects for normalization.");
 706   }
 707
 708
 709 }
 710 //________________________________________________________________________
 711 TH1D*
 712 AliBasedNdetaTask::Rebin(const TH1D* h) const
 713 {
 714   //
 715   // Make a copy of the input histogram and rebin that histogram
 716   //
 717   // Parameters:
 718   //    h  Histogram to rebin
 719   //
 720   // Return:
 721   //    New (rebinned) histogram
 722   //
 723   if (fRebin <= 1) return 0;
 724
 725   Int_t nBins = h->GetNbinsX();
 726   if(nBins % fRebin != 0) {
 727     Warning("Rebin", "Rebin factor %d is not a devisor of current number "
 728             "of bins %d in the histogram %s", fRebin, nBins, h->GetName());
 729     return 0;
 730   }
 731
 732   // Make a copy
 733   TH1D* tmp = static_cast<TH1D*>(h->Clone(Form("%s_rebin%02d",
 734                                                h->GetName(), fRebin)));
 735   tmp->Rebin(fRebin);
 736   tmp->SetDirectory(0);
 737
 738   // The new number of bins
 739   Int_t nBinsNew = nBins / fRebin;
 740   for(Int_t i = 1;i<= nBinsNew; i++) {
 741     Double_t content = 0;
 742     Double_t sumw    = 0;
 743     Double_t wsum    = 0;
 744     Int_t    nbins   = 0;
 745     for(Int_t j = 1; j<=fRebin;j++) {
 746       Int_t    bin = (i-1)*fRebin + j;
 747       Double_t c   =  h->GetBinContent(bin);
 748       if (c <= 0) continue;
 749
 750       if (fCutEdges) {
 751         if (h->GetBinContent(bin+1)<=0 ||
 752             h->GetBinContent(bin-1)<=0) {
 753           Warning("Rebin", "removing bin %d=%f of %s (%d=%f,%d=%f)",
 754                   bin, c, h->GetName(),
 755                   bin+1, h->GetBinContent(bin+1),
 756                   bin-1, h->GetBinContent(bin-1));
 757           continue;
 758         }
 759       }
 760       Double_t e =  h->GetBinError(bin);
 761       Double_t w =  1 / (e*e); // 1/c/c
 762       content    += c;
 763       sumw       += w;
 764       wsum       += w * c;
 765       nbins++;
 766     }
 767
 768     if(content > 0 && nbins > 0) {
 769       tmp->SetBinContent(i, wsum / sumw);
 770       tmp->SetBinError(i,1./TMath::Sqrt(sumw));
 771     }
 772   }
 773
 774   return tmp;
 775 }
 776
 777 //__________________________________________________________________
 778 TH1*
 779 AliBasedNdetaTask::Symmetrice(const TH1* h) const
 780 {
 781   //
 782   // Make an extension of @a h to make it symmetric about 0
 783   //
 784   // Parameters:
 785   //    h Histogram to symmertrice
 786   //
 787   // Return:
 788   //    Symmetric extension of @a h
 789   //
 790   Int_t nBins = h->GetNbinsX();
 791   TH1* s = static_cast<TH1*>(h->Clone(Form("%s_mirror", h->GetName())));
 792   s->SetTitle(Form("%s (mirrored)", h->GetTitle()));
 793   s->Reset();
 794   s->SetBins(nBins, -h->GetXaxis()->GetXmax(), -h->GetXaxis()->GetXmin());
 795   s->SetMarkerColor(h->GetMarkerColor());
 796   s->SetMarkerSize(h->GetMarkerSize());
 797   s->SetMarkerStyle(h->GetMarkerStyle()+4);
 798   s->SetFillColor(h->GetFillColor());
 799   s->SetFillStyle(h->GetFillStyle());
 800   s->SetDirectory(0);
 801
 802   // Find the first and last bin with data
 803   Int_t first = nBins+1;
 804   Int_t last  = 0;
 805   for (Int_t i = 1; i <= nBins; i++) {
 806     if (h->GetBinContent(i) <= 0) continue;
 807     first = TMath::Min(first, i);
 808     last  = TMath::Max(last,  i);
 809   }
 810
 811   Double_t xfirst = h->GetBinCenter(first-1);
 812   Int_t    f1     = h->GetXaxis()->FindBin(-xfirst);
 813   Int_t    l2     = s->GetXaxis()->FindBin(xfirst);
 814   for (Int_t i = f1, j=l2; i <= last; i++,j--) {
 815     s->SetBinContent(j, h->GetBinContent(i));
 816     s->SetBinError(j, h->GetBinError(i));
 817   }
 818   // Fill in overlap bin
 819   s->SetBinContent(l2+1, h->GetBinContent(first));
 820   s->SetBinError(l2+1, h->GetBinError(first));
 821   return s;
 822 }
 823 //
 824 // EOF
 825 //