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