PWGLF/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 <AliAnalysisManager.h>
   9 #include <AliAODEvent.h>
  10 #include <AliAODHandler.h>
  11 #include <AliAODInputHandler.h>
  12 #include "AliForwardUtil.h"
  13 #include "AliAODForwardMult.h"
  14 #include <TFile.h>
  15 #include <TStyle.h>
  16
  17 //____________________________________________________________________
  18 AliBasedNdetaTask::AliBasedNdetaTask()
  19   : AliAnalysisTaskSE(),
  20     fSums(0),           // Container of sums
  21     fOutput(0),         // Container of output
  22     fVtxMin(0),         // Minimum v_z
  23     fVtxMax(0),         // Maximum v_z
  24     fTriggerMask(0),    // Trigger mask
  25     fRebin(0),          // Rebinning factor
  26     fCutEdges(false),
  27     fSymmetrice(true),
  28     fCorrEmpty(true),
  29     fUseROOTProj(false),
  30     fTriggerEff(1),
  31     fTriggerEff0(1),
  32     fShapeCorr(0),
  33     fListOfCentralities(0),
  34     fSNNString(0),
  35     fSysString(0),
  36     fCent(0),
  37     fCentAxis(0),
  38     fVtx(0),
  39     fNormalizationScheme(kFull),
  40     fSchemeString(0),
  41     fTriggerString(0),
  42     fFinalMCCorrFile(""),
  43     fSatelliteVertices(0),
  44     fglobalempiricalcorrection(0),
  45     fmeabsignalvscentr(0)
  46 {
  47   //
  48   // Constructor
  49   //
  50   DGUARD(fDebug,3,"Default CTOR of AliBasedNdetaTask");
  51 }
  52
  53 //____________________________________________________________________
  54 AliBasedNdetaTask::AliBasedNdetaTask(const char* name)
  55   : AliAnalysisTaskSE(name),
  56     fSums(0),           // Container of sums
  57     fOutput(0),         // Container of output
  58     fVtxMin(-10),       // Minimum v_z
  59     fVtxMax(10),        // Maximum v_z
  60     fTriggerMask(AliAODForwardMult::kInel),
  61     fRebin(5),          // Rebinning factor
  62     fCutEdges(false),
  63     fSymmetrice(true),
  64     fCorrEmpty(true),
  65     fUseROOTProj(false),
  66     fTriggerEff(1),
  67     fTriggerEff0(1),
  68     fShapeCorr(0),
  69     fListOfCentralities(0),
  70     fSNNString(0),
  71     fSysString(0),
  72     fCent(0),
  73     fCentAxis(0),
  74     fVtx(0),
  75     fNormalizationScheme(kFull),
  76     fSchemeString(0),
  77     fTriggerString(0),
  78     fFinalMCCorrFile(""),
  79     fSatelliteVertices(0),
  80     fglobalempiricalcorrection(0),
  81     fmeabsignalvscentr(0)
  82 {
  83   //
  84   // Constructor
  85   //
  86   DGUARD(fDebug, 3,"Named CTOR of AliBasedNdetaTask: %s", name);
  87   fListOfCentralities = new TObjArray(1);
  88
  89   // Set the normalisation scheme
  90   SetNormalizationScheme(kFull);
  91
  92   // Set the trigger mask
  93   SetTriggerMask(AliAODForwardMult::kInel);
  94
  95   // Output slot #1 writes into a TH1 container
  96   DefineOutput(1, TList::Class());
  97   DefineOutput(2, TList::Class());
  98 }
  99
 100 //____________________________________________________________________
 101 AliBasedNdetaTask::AliBasedNdetaTask(const AliBasedNdetaTask& o)
 102   : AliAnalysisTaskSE(o),
 103     fSums(o.fSums),             // TList* - Container of sums
 104     fOutput(o.fOutput),         // Container of output
 105     fVtxMin(o.fVtxMin),         // Double_t - Minimum v_z
 106     fVtxMax(o.fVtxMax),         // Double_t - Maximum v_z
 107     fTriggerMask(o.fTriggerMask),// Int_t - Trigger mask
 108     fRebin(o.fRebin),           // Int_t - Rebinning factor
 109     fCutEdges(o.fCutEdges),     // Bool_t - Whether to cut edges when rebinning
 110     fSymmetrice(o.fSymmetrice),
 111     fCorrEmpty(o.fCorrEmpty),
 112     fUseROOTProj(o.fUseROOTProj),
 113     fTriggerEff(o.fTriggerEff),
 114     fTriggerEff0(o.fTriggerEff0),
 115     fShapeCorr(o.fShapeCorr),
 116     fListOfCentralities(o.fListOfCentralities),
 117     fSNNString(o.fSNNString),
 118     fSysString(o.fSysString),
 119     fCent(o.fCent),
 120     fCentAxis(o.fCentAxis),
 121     fVtx(o.fVtx),
 122     fNormalizationScheme(o.fNormalizationScheme),
 123     fSchemeString(o.fSchemeString),
 124     fTriggerString(o.fTriggerString),
 125     fFinalMCCorrFile(o.fFinalMCCorrFile),
 126     fSatelliteVertices(o.fSatelliteVertices),
 127     fglobalempiricalcorrection(o.fglobalempiricalcorrection),
 128   fmeabsignalvscentr(o.fmeabsignalvscentr)
 129 {
 130   DGUARD(fDebug, 3,"Copy CTOR of AliBasedNdetaTask");
 131 }
 132
 133 //____________________________________________________________________
 134 AliBasedNdetaTask::~AliBasedNdetaTask()
 135 {
 136   //
 137   // Destructor
 138   //
 139   DGUARD(fDebug,3,"Destruction of AliBasedNdetaTask");
 140   if (fSums) {
 141     fSums->Delete();
 142     delete fSums;
 143     fSums = 0;
 144   }
 145   if (fOutput) {
 146     fOutput->Delete();
 147     delete fOutput;
 148     fOutput = 0;
 149   }
 150 }
 151
 152 //________________________________________________________________________
 153 void
 154 AliBasedNdetaTask::SetDebugLevel(Int_t lvl)
 155 {
 156   AliAnalysisTaskSE::SetDebugLevel(lvl);
 157   for (Int_t i = 0; i < fListOfCentralities->GetEntries(); i++) {
 158     CentralityBin* bin =
 159       static_cast<CentralityBin*>(fListOfCentralities->At(i));
 160     bin->SetDebugLevel(lvl);
 161   }
 162 }
 163
 164 //________________________________________________________________________
 165 void
 166 AliBasedNdetaTask::SetCentralityAxis(UShort_t n, Short_t* bins)
 167 {
 168   DGUARD(fDebug,3,"Set centrality axis, %d bins", n);
 169   if (!fCentAxis) {
 170     fCentAxis = new TAxis();
 171     fCentAxis->SetName("centAxis");
 172     fCentAxis->SetTitle("Centrality [%]");
 173   }
 174   TArrayD dbins(n+1);
 175   for (UShort_t i = 0; i <= n; i++)
 176     dbins[i] = (bins[i] == 100 ? 100.1 : bins[i]);
 177   fCentAxis->Set(n, dbins.GetArray());
 178 }
 179
 180 //________________________________________________________________________
 181 void
 182 AliBasedNdetaTask::AddCentralityBin(UShort_t at, Short_t low, Short_t high)
 183 {
 184   //
 185   // Add a centrality bin
 186   //
 187   // Parameters:
 188   //    low  Low cut
 189   //    high High cut
 190   //
 191   DGUARD(fDebug,3,"Add a centrality bin [%d,%d] @ %d", low, high, at);
 192   CentralityBin* bin = MakeCentralityBin(GetName(), low, high);
 193   if (!bin) {
 194     Error("AddCentralityBin",
 195           "Failed to create centrality bin for %s [%d,%d] @ %d",
 196           GetName(), low, high, at);
 197     return;
 198   }
 199   bin->SetSatelliteVertices(fSatelliteVertices);
 200   bin->SetDebugLevel(fDebug);
 201   fListOfCentralities->AddAtAndExpand(bin, at);
 202 }
 203
 204 //________________________________________________________________________
 205 AliBasedNdetaTask::CentralityBin*
 206 AliBasedNdetaTask::MakeCentralityBin(const char* name,
 207                                      Short_t low, Short_t high) const
 208 {
 209   //
 210   // Make a centrality bin
 211   //
 212   // Parameters:
 213   //    name  Name used for histograms
 214   //    low   Low cut in percent
 215   //    high  High cut in percent
 216   //
 217   // Return:
 218   //    A newly created centrality bin
 219   //
 220   DGUARD(fDebug,3,"Make a centrality bin %s [%d,%d]", name, low, high);
 221   return new CentralityBin(name, low, high);
 222 }
 223
 224 #define TESTAPPEND(SCHEME,BIT,STRING) \
 225   do { if (!(SCHEME & BIT)) break;                                      \
 226     if (!s.IsNull()) s.Append(","); s.Append(STRING); } while(false)
 227
 228 //________________________________________________________________________
 229 const Char_t*
 230 AliBasedNdetaTask::NormalizationSchemeString(UShort_t scheme)
 231 {
 232   // Create a string from normalization scheme bits
 233   static TString s;
 234   s = "";
 235
 236   if (scheme == kNone)
 237     return s.Data();
 238   if (scheme == kFull) {
 239     s = "FULL";
 240     return s.Data();
 241   }
 242   TESTAPPEND(scheme, kEventLevel,        "EVENT");
 243   TESTAPPEND(scheme, kShape,             "SHAPE");
 244   TESTAPPEND(scheme, kBackground,        "BACKGROUND");
 245   TESTAPPEND(scheme, kTriggerEfficiency, "TRIGGER");
 246   TESTAPPEND(scheme, kZeroBin,           "ZEROBIN");
 247
 248   return s.Data();
 249 }
 250 //________________________________________________________________________
 251 UShort_t
 252 AliBasedNdetaTask::ParseNormalizationScheme(const char* what)
 253 {
 254   UShort_t    scheme = 0;
 255   TString     twhat(what);
 256   twhat.ToUpper();
 257   TObjString* opt;
 258   TObjArray* token = twhat.Tokenize(" ,|");
 259   TIter       next(token);
 260   while ((opt = static_cast<TObjString*>(next()))) {
 261     TString s(opt->GetString());
 262     if      (s.IsNull()) continue;
 263     Bool_t add = true;
 264     switch (s[0]) {
 265     case '-': add = false; // Fall through
 266     case '+': s.Remove(0,1);  // Remove character
 267     }
 268     UShort_t bit = 0;
 269     if      (s.CompareTo("EVENT")     == 0) bit = kEventLevel;
 270     else if (s.CompareTo("SHAPE")     == 0) bit = kShape;
 271     else if (s.CompareTo("BACKGROUND")== 0) bit = kBackground;
 272     else if (s.CompareTo("TRIGGER")   == 0) bit = kTriggerEfficiency;
 273     else if (s.CompareTo("FULL")      == 0) bit = kFull;
 274     else if (s.CompareTo("NONE")      == 0) bit = kNone;
 275     else if (s.CompareTo("ZEROBIN")   == 0) bit = kZeroBin;
 276     else
 277       ::Warning("SetNormalizationScheme", "Unknown option %s", s.Data());
 278     if (add) scheme |= bit;
 279     else     scheme ^= bit;
 280   }
 281   delete token;
 282   return scheme;
 283 }
 284 //________________________________________________________________________
 285 void
 286 AliBasedNdetaTask::SetNormalizationScheme(const char* what)
 287 {
 288   //
 289   // Set normalisation scheme
 290   //
 291   DGUARD(fDebug,3,"Set the normalization scheme: %s", what);
 292   SetNormalizationScheme(ParseNormalizationScheme(what));
 293 }
 294 //________________________________________________________________________
 295 void
 296 AliBasedNdetaTask::SetNormalizationScheme(UShort_t scheme)
 297 {
 298   DGUARD(fDebug,3,"Set the normalization scheme: 0x%x", scheme);
 299   fNormalizationScheme = scheme;
 300   if (fSchemeString) delete fSchemeString;
 301   fSchemeString = AliForwardUtil::MakeParameter("scheme", scheme);
 302 }
 303 //________________________________________________________________________
 304 void
 305 AliBasedNdetaTask::SetTriggerMask(const char* mask)
 306 {
 307   //
 308   // Set the trigger maskl
 309   //
 310   // Parameters:
 311   //    mask Trigger mask
 312   //
 313   DGUARD(fDebug,3,"Set the trigger mask: %s", mask);
 314   SetTriggerMask(AliAODForwardMult::MakeTriggerMask(mask));
 315 }
 316 //________________________________________________________________________
 317 void
 318 AliBasedNdetaTask::SetTriggerMask(UShort_t mask)
 319 {
 320   DGUARD(fDebug,3,"Set the trigger mask: 0x%0x", mask);
 321   fTriggerMask = mask;
 322   if (fTriggerString) delete fTriggerString;
 323   fTriggerString = AliForwardUtil::MakeParameter("trigger", fTriggerMask);
 324 }
 325
 326 //________________________________________________________________________
 327 void
 328 AliBasedNdetaTask::SetShapeCorrection(const TH2F* c)
 329 {
 330   //
 331   // Set the shape correction (a.k.a., track correction) for selected
 332   // trigger(s)
 333   //
 334   // Parameters:
 335   //    h Correction
 336   //
 337   DGUARD(fDebug,3,"Set the shape correction: %p", c);
 338   if (!c) return;
 339   fShapeCorr = static_cast<TH2F*>(c->Clone());
 340   fShapeCorr->SetDirectory(0);
 341 }
 342 //________________________________________________________________________
 343 void
 344 AliBasedNdetaTask::InitializeCentBins()
 345 {
 346   if (fListOfCentralities->GetEntries() > 0) return;
 347
 348   // Automatically add 'all' centrality bin if nothing has been defined.
 349   AddCentralityBin(0, 0, 0);
 350   if (fCentAxis && fCentAxis->GetNbins() > 0 && fCentAxis->GetXbins()) {
 351     const TArrayD* bins = fCentAxis->GetXbins();
 352     Int_t          nbin = fCentAxis->GetNbins();
 353     for (Int_t i = 0; i < nbin; i++)
 354       AddCentralityBin(i+1,  Short_t((*bins)[i]), Short_t((*bins)[i+1]));
 355   }
 356 }
 357
 358 //________________________________________________________________________
 359 void
 360 AliBasedNdetaTask::UserCreateOutputObjects()
 361 {
 362   //
 363   // Create output objects.
 364   //
 365   // This is called once per slave process
 366   //
 367   DGUARD(fDebug,1,"Create user ouput object");
 368   fSums = new TList;
 369   fSums->SetName(Form("%s_sums", GetName()));
 370   fSums->SetOwner();
 371
 372   InitializeCentBins();
 373   if (fCentAxis) fSums->Add(fCentAxis);
 374
 375   fSums->Add(AliForwardUtil::MakeParameter("alirootRev",
 376                                            AliForwardUtil::AliROOTRevision()));
 377   fSums->Add(AliForwardUtil::MakeParameter("alirootBranch",
 378                                            AliForwardUtil::AliROOTBranch()));
 379   fSums->Add(AliForwardUtil::MakeParameter("empirical",
 380                                            fglobalempiricalcorrection != 0));
 381   // Centrality histogram
 382   fCent = new TH1D("cent", "Centrality", 100, 0, 100);
 383   fCent->SetDirectory(0);
 384   fCent->SetXTitle(0);
 385   fSums->Add(fCent);
 386
 387   // Custom vertex axis that will include satellite vertices
 388   // Satellite vertices are at k*37.5 where k=-10,-9,...,9,10
 389   // Nominal vertices are usually in -10 to 10 and we should have
 390   // 10 bins in that range.  That gives us a total of
 391   //
 392   //   10+10+10=30 bins
 393   //
 394   // or 31 bin boundaries
 395   const Int_t   nCenter = 20; // Max possible is 56
 396   const Int_t   nSat    = 10;
 397   const Int_t   nBins   = 2*nSat + nCenter;
 398   Float_t       bins[nBins+1];
 399   for (Int_t i = 0; i < nSat;i++) {
 400     bins[i] = (i-nSat-.5) * 37.5;
 401     // printf("bins[%2d]=%+6.2f\n", i, bins[i]);
 402   }
 403   for (Int_t i = nSat; i < nSat+nCenter+1; i++) {
 404     bins[i] = -nCenter + (i-nSat) * 2;
 405     // printf("bins[%2d]=%+6.2f\n", i, bins[i]);
 406   }
 407   for (Int_t i = nSat+nCenter+1; i < 2*nSat+nCenter; i++) {
 408     bins[i] = (i-nSat-nCenter +.5) * 37.5;
 409     // printf("bins[%2d]=%+6.2f\n", i, bins[i]);
 410   }
 411   bins[nBins] = (nSat + .5) * 37.5;
 412
 413   fVtx           = new TH1D("vtx", "Vertex dist", nBins, bins);
 414   fVtx->SetDirectory(0);
 415   fVtx->SetFillColor(kRed+1);
 416   fVtx->SetFillStyle(3001);
 417   fVtx->SetXTitle("IP_{z} [cm]");
 418   fVtx->SetYTitle("Events");
 419   fSums->Add(fVtx);
 420 #if 0
 421   TAxis* a = fVtx->GetXaxis();
 422   for (Int_t i = 1; i <= nBins; i++)
 423     printf("%2d/%2d: %+6.2f - %+6.2f: %+6.2f\n",
 424            i,nBins,a->GetBinLowEdge(i),a->GetBinUpEdge(i),a->GetBinCenter(i));
 425 #endif
 426
 427   // Loop over centrality bins
 428   TIter next(fListOfCentralities);
 429   CentralityBin* bin = 0;
 430   while ((bin = static_cast<CentralityBin*>(next())))
 431     bin->CreateOutputObjects(fSums, fTriggerMask);
 432
 433   fmeabsignalvscentr=new TH2D("meanAbsSignalVsCentr",
 434                               "Mean absolute signal versus centrality",
 435                               400, 0, 20, 100, 0, 100);
 436   fSums->Add(fmeabsignalvscentr);
 437
 438   Print();
 439   // Post data for ALL output slots >0 here, to get at least an empty
 440   // histogram
 441   PostData(1, fSums);
 442 }
 443
 444 //____________________________________________________________________
 445 void
 446 AliBasedNdetaTask::UserExec(Option_t *)
 447 {
 448   //
 449   // Process a single event
 450   //
 451   // Parameters:
 452   //    option Not used
 453   //
 454   // Main loop
 455   DGUARD(fDebug,1,"Analyse the AOD event");
 456   AliAODEvent* aod = AliForwardUtil::GetAODEvent(this);
 457   if (!aod) return;
 458
 459
 460   TObject* obj = aod->FindListObject("Forward");
 461   if (!obj) {
 462     AliWarning("No forward object found");
 463     return;
 464   }
 465   AliAODForwardMult* forward = static_cast<AliAODForwardMult*>(obj);
 466
 467   // Fill centrality histogram
 468
 469   Double_t vtx = forward->GetIpZ();
 470   Float_t cent = forward->GetCentrality();
 471
 472   // Get the histogram(s)
 473   TH2D* data   = GetHistogram(aod, false);
 474   TH2D* dataMC = GetHistogram(aod, true);
 475   if (!data) return;
 476
 477   CheckEventData(vtx, data, dataMC);
 478
 479   if (!ApplyEmpiricalCorrection(forward,data))
 480     return;
 481
 482 #if 0
 483   // Code disabled - breaks execution
 484   Int_t notemptybins=0;
 485   Double_t sum=0.0;
 486   for (Int_t ix=1;ix<=data->GetXaxis()->GetNbins();ix++)
 487     {
 488       Double_t sumy=0.0;
 489       for(Int_t iy=1;iy<=data->GetYaxis()->GetNbins();iy++)
 490         {
 491           if(data->GetBinContent(ix,iy)>0.0)
 492             {
 493               sumy+=data->GetBinContent(ix,iy);
 494               notemptybins++;
 495             }
 496
 497         }
 498       sum+=sumy;
 499     }
 500
 501   if(notemptybins>0)
 502     {
 503       sum=sum/((Double_t)notemptybins);
 504     }
 505   else
 506     sum=-1.0;
 507   fmeabsignalvscentr->Fill(sum,cent);
 508 #endif
 509
 510   Bool_t isZero = ((fNormalizationScheme & kZeroBin) &&
 511                    !forward->IsTriggerBits(AliAODForwardMult::kNClusterGt0));
 512   Bool_t taken = false;
 513
 514   // Loop over centrality bins
 515   CentralityBin* allBin =
 516     static_cast<CentralityBin*>(fListOfCentralities->At(0));
 517   if (allBin->ProcessEvent(forward, fTriggerMask, isZero,
 518                            fVtxMin, fVtxMax, data, dataMC)) taken = true;
 519
 520   // Find this centrality bin
 521   if (fCentAxis && fCentAxis->GetNbins() > 0) {
 522     Int_t          icent   = fCentAxis->FindBin(cent);
 523     CentralityBin* thisBin = 0;
 524     if (icent >= 1 && icent <= fCentAxis->GetNbins())
 525       thisBin = static_cast<CentralityBin*>(fListOfCentralities->At(icent));
 526     if (thisBin)
 527       if (thisBin->ProcessEvent(forward, fTriggerMask, isZero, fVtxMin,
 528                                 fVtxMax, data, dataMC)) taken = true;
 529   }
 530
 531   // Fill diagnostics only if we took this event
 532   if (taken) {
 533     fCent->Fill(cent);
 534     fVtx->Fill(vtx);
 535   }
 536
 537   // Here, we get the update
 538   if (!fSNNString) {
 539     fSNNString = AliForwardUtil::MakeParameter("sNN", forward->GetSNN());
 540     fSysString = AliForwardUtil::MakeParameter("sys", forward->GetSystem());
 541
 542     fSums->Add(fSNNString);
 543     fSums->Add(fSysString);
 544     fSums->Add(fSchemeString);
 545     fSums->Add(fTriggerString);
 546
 547     // Show stuff on first event
 548     Print();
 549   }
 550   PostData(1, fSums);
 551 }
 552
 553 //________________________________________________________________________
 554 void AliBasedNdetaTask::CheckEventData(Double_t,
 555                                        TH2*,
 556                                        TH2*)
 557 {
 558 }
 559
 560 //________________________________________________________________________
 561 void
 562 AliBasedNdetaTask::SetHistogramAttributes(TH1D* h, Int_t colour, Int_t marker,
 563                                           const char* title, const char* ytitle)
 564 {
 565   //
 566   // Set histogram graphical options, etc.
 567   //
 568   // Parameters:
 569   //    h       Histogram to modify
 570   //    colour  Marker color
 571   //    marker  Marker style
 572   //    title   Title of histogram
 573   //    ytitle  Title on y-axis.
 574   //
 575   h->SetTitle(title);
 576   h->SetMarkerColor(colour);
 577   h->SetMarkerStyle(marker);
 578   h->SetMarkerSize(marker == 29 || marker == 30 ? 1.2 : 1);
 579   h->SetFillStyle(0);
 580   h->SetYTitle(ytitle);
 581   h->GetXaxis()->SetTitleFont(132);
 582   h->GetXaxis()->SetLabelFont(132);
 583   h->GetXaxis()->SetNdivisions(10);
 584   h->GetYaxis()->SetTitleFont(132);
 585   h->GetYaxis()->SetLabelFont(132);
 586   h->GetYaxis()->SetNdivisions(10);
 587   h->GetYaxis()->SetDecimals();
 588   h->SetStats(0);
 589 }
 590
 591 //________________________________________________________________________
 592 void
 593 AliBasedNdetaTask::ScaleToCoverage(TH2D* copy, const TH1D* norm)
 594 {
 595   // Normalize to the acceptance -
 596   // dndeta->Divide(accNorm);
 597   for (Int_t i = 1; i <= copy->GetNbinsX(); i++) {
 598     Double_t a = norm->GetBinContent(i);
 599     for (Int_t j = 1; j <= copy->GetNbinsY(); j++) {
 600       if (a <= 0) {
 601         copy->SetBinContent(i,j,0);
 602         copy->SetBinError(i,j,0);
 603         continue;
 604       }
 605       Double_t c = copy->GetBinContent(i, j);
 606       Double_t e = copy->GetBinError(i, j);
 607       copy->SetBinContent(i, j, c / a);
 608       copy->SetBinError(i, j, e / a);
 609     }
 610   }
 611 }
 612 //________________________________________________________________________
 613 void
 614 AliBasedNdetaTask::ScaleToCoverage(TH1D* copy, const TH1D* norm)
 615 {
 616   // Normalize to the acceptance -
 617   // dndeta->Divide(accNorm);
 618   for (Int_t i = 1; i <= copy->GetNbinsX(); i++) {
 619     Double_t a = norm->GetBinContent(i);
 620     if (a <= 0) {
 621       copy->SetBinContent(i,0);
 622       copy->SetBinError(i,0);
 623       continue;
 624     }
 625     Double_t c = copy->GetBinContent(i);
 626     Double_t e = copy->GetBinError(i);
 627     copy->SetBinContent(i, c / a);
 628     copy->SetBinError(i, e / a);
 629   }
 630 }
 631
 632 //________________________________________________________________________
 633 TH1D*
 634 AliBasedNdetaTask::ProjectX(const TH2D* h,
 635                             const char* name,
 636                             Int_t firstbin,
 637                             Int_t lastbin,
 638                             bool  useRoot,
 639                             bool  corr,
 640                             bool  error)
 641 {
 642   //
 643   // Project onto the X axis
 644   //
 645   // Parameters:
 646   //    h         2D histogram
 647   //    name      New name
 648   //    firstbin  First bin to use
 649   //    lastbin   Last bin to use
 650   //    error     Whether to calculate errors
 651   //
 652   // Return:
 653   //    Newly created histogram or null
 654   //
 655   if (!h) return 0;
 656   if (useRoot)
 657     return h->ProjectionX(name, firstbin, lastbin, (error ? "e" : ""));
 658
 659   TAxis* xaxis = h->GetXaxis();
 660   TAxis* yaxis = h->GetYaxis();
 661   TH1D*  ret   = new TH1D(name, h->GetTitle(), xaxis->GetNbins(),
 662                           xaxis->GetXmin(), xaxis->GetXmax());
 663   static_cast<const TAttLine*>(h)->Copy(*ret);
 664   static_cast<const TAttFill*>(h)->Copy(*ret);
 665   static_cast<const TAttMarker*>(h)->Copy(*ret);
 666   ret->GetXaxis()->ImportAttributes(xaxis);
 667
 668   Int_t first = firstbin;
 669   Int_t last  = lastbin;
 670   if      (first < 0)                    first = 1;
 671   else if (first >= yaxis->GetNbins()+2) first = yaxis->GetNbins()+1;
 672   if      (last  < 0)                    last  = yaxis->GetNbins();
 673   else if (last  >= yaxis->GetNbins()+2) last  = yaxis->GetNbins()+1;
 674   if (last-first < 0) {
 675     AliWarningGeneral("AliBasedNdetaTask",
 676                       Form("Nothing to project [%d,%d]", first, last));
 677     return 0;
 678
 679   }
 680
 681   // Loop over X bins
 682   // AliInfo(Form("Projecting bins [%d,%d] of %s", first, last, h->GetName()));
 683   Int_t ybins = (last-first+1);
 684   for (Int_t xbin = 0; xbin <= xaxis->GetNbins()+1; xbin++) {
 685     Double_t content = 0;
 686     Double_t error2  = 0;
 687     Int_t    nbins   = 0;
 688
 689
 690     for (Int_t ybin = first; ybin <= last; ybin++) {
 691       Double_t c1 = h->GetCellContent(xbin, ybin);
 692       Double_t e1 = h->GetCellError(xbin, ybin);
 693
 694       // Ignore empty bins
 695       if (c1 < 1e-12) continue;
 696       if (e1 < 1e-12) {
 697         if (error) continue;
 698         e1 = 1;
 699       }
 700
 701       content    += c1;
 702       error2     += e1*e1;
 703       nbins++;
 704     } // for (ybin)
 705     if(content > 0 && nbins > 0) {
 706       Double_t factor = (corr ? Double_t(ybins) / nbins : 1);
 707 #if 0
 708       AliWarningGeneral(ret->GetName(),
 709                         Form("factor @ %d is %d/%d -> %f",
 710                              xbin, ybins, nbins, factor));
 711 #endif
 712       if (error) {
 713         // Calculate weighted average
 714         ret->SetBinContent(xbin, content * factor);
 715         ret->SetBinError(xbin, factor * TMath::Sqrt(error2));
 716       }
 717       else
 718         ret->SetBinContent(xbin, factor * content);
 719     }
 720   } // for (xbin)
 721
 722   return ret;
 723 }
 724
 725 //________________________________________________________________________
 726 void
 727 AliBasedNdetaTask::Terminate(Option_t *)
 728 {
 729   //
 730   // Called at end of event processing..
 731   //
 732   // This is called once in the master
 733   //
 734   // Parameters:
 735   //    option Not used
 736   //
 737   // Draw result to screen, or perform fitting, normalizations Called
 738   // once at the end of the query
 739   DGUARD(fDebug,1,"Process final merged results");
 740
 741   fSums = dynamic_cast<TList*> (GetOutputData(1));
 742   if(!fSums) {
 743     AliError("Could not retrieve TList fSums");
 744     return;
 745   }
 746
 747   fOutput = new TList;
 748   fOutput->SetName(Form("%s_result", GetName()));
 749   fOutput->SetOwner();
 750
 751   fSNNString     = fSums->FindObject("sNN");
 752   fSysString     = fSums->FindObject("sys");
 753   fCentAxis      = static_cast<TAxis*>(fSums->FindObject("centAxis"));
 754   fSchemeString  = fSums->FindObject("scheme");
 755   fTriggerString = fSums->FindObject("trigger");
 756
 757   if(fSysString && fSNNString &&
 758      fSysString->GetUniqueID() == AliForwardUtil::kPP)
 759     LoadNormalizationData(fSysString->GetUniqueID(),
 760                           fSNNString->GetUniqueID());
 761
 762   InitializeCentBins();
 763
 764   // Print before we loop
 765   Print();
 766
 767   // Loop over centrality bins
 768   TIter next(fListOfCentralities);
 769   CentralityBin* bin = 0;
 770   gStyle->SetPalette(1);
 771   THStack* dndetaStack        = new THStack("dndeta", "dN/d#eta");
 772   THStack* dndetaStackRebin   = new THStack(Form("dndeta_rebin%02d", fRebin),
 773                                             "dN_{ch}/d#eta");
 774   THStack* dndetaMCStack      = new THStack("dndetaMC", "dN_{ch}/d#eta");
 775   THStack* dndetaMCStackRebin = new THStack(Form("dndetaMC_rebin%02d", fRebin),
 776                                             "dN_{ch}/d#eta");
 777
 778   TList* mclist = 0;
 779   TList* truthlist = 0;
 780
 781   if (fFinalMCCorrFile.Contains(".root")) {
 782     TFile* ftest = TFile::Open(fFinalMCCorrFile.Data());
 783     if(ftest) {
 784       mclist    = dynamic_cast<TList*>(ftest->Get(Form("%sResults",GetName())));
 785       truthlist = dynamic_cast<TList*>(ftest->Get("MCTruthResults"));
 786     }
 787     else
 788       AliWarning("MC analysis file invalid - no final MC correction possible");
 789   }
 790   Int_t style = GetMarker();
 791   Int_t color = GetColor();
 792
 793   AliInfo(Form("Marker style=%d, color=%d", style, color));
 794   while ((bin = static_cast<CentralityBin*>(next()))) {
 795
 796     bin->End(fSums, fOutput, fNormalizationScheme, fShapeCorr,
 797              fTriggerEff, fTriggerEff0,
 798              fSymmetrice, fRebin, fUseROOTProj, fCorrEmpty, fCutEdges,
 799              fTriggerMask, style, color, mclist, truthlist);
 800     if (fCentAxis && bin->IsAllBin()) continue;
 801     TH1* dndeta      = bin->GetResult(0, false, "");
 802     TH1* dndetaSym   = bin->GetResult(0, true,  "");
 803     TH1* dndetaMC    = bin->GetResult(0, false, "MC");
 804     TH1* dndetaMCSym = bin->GetResult(0, true,  "MC");
 805     if (dndeta)      dndetaStack->Add(dndeta);
 806     if (dndetaSym)   dndetaStack->Add(dndetaSym);
 807     if (dndetaMC)    dndetaMCStack->Add(dndetaMC);
 808     if (dndetaMCSym) dndetaMCStack->Add(dndetaMCSym);
 809     if (fRebin > 1) {
 810       dndeta      = bin->GetResult(fRebin, false, "");
 811       dndetaSym   = bin->GetResult(fRebin, true,  "");
 812       dndetaMC    = bin->GetResult(fRebin, false, "MC");
 813       dndetaMCSym = bin->GetResult(fRebin, true,  "MC");
 814       if (dndeta)      dndetaStackRebin->Add(dndeta);
 815       if (dndetaSym)   dndetaStackRebin->Add(dndetaSym);
 816       if (dndetaMC)    dndetaMCStackRebin->Add(dndetaMC);
 817       if (dndetaMCSym) dndetaMCStackRebin->Add(dndetaMCSym);
 818     }
 819   }
 820   // Output the stack
 821   fOutput->Add(dndetaStack);
 822
 823   // If available output rebinned stack
 824   if (!dndetaStackRebin->GetHists() ||
 825       dndetaStackRebin->GetHists()->GetEntries() <= 0) {
 826     AliWarning("No rebinned histograms found");
 827     delete dndetaStackRebin;
 828     dndetaStackRebin = 0;
 829   }
 830   if (dndetaStackRebin) fOutput->Add(dndetaStackRebin);
 831
 832   // If available, output track-ref stack
 833   if (!dndetaMCStack->GetHists() ||
 834       dndetaMCStack->GetHists()->GetEntries() <= 0) {
 835     AliWarning("No MC histograms found");
 836     delete dndetaMCStack;
 837     dndetaMCStack = 0;
 838   }
 839   if (dndetaMCStack) fOutput->Add(dndetaMCStack);
 840
 841   // If available, output rebinned track-ref stack
 842   if (!dndetaMCStackRebin->GetHists() ||
 843       dndetaMCStackRebin->GetHists()->GetEntries() <= 0) {
 844     AliWarning("No rebinned MC histograms found");
 845     delete dndetaMCStackRebin;
 846     dndetaMCStackRebin = 0;
 847   }
 848   if (dndetaMCStackRebin) fOutput->Add(dndetaMCStackRebin);
 849
 850   // Output collision energy string
 851   if (fSNNString) {
 852     UShort_t sNN = fSNNString->GetUniqueID();
 853     TNamed* sNNObj = new TNamed(fSNNString->GetName(),
 854                                 AliForwardUtil::CenterOfMassEnergyString(sNN));
 855     sNNObj->SetUniqueID(sNN);
 856     fOutput->Add(sNNObj); // fSNNString->Clone());
 857   }
 858
 859   // Output collision system string
 860   if (fSysString) {
 861     UShort_t sys = fSysString->GetUniqueID();
 862     TNamed* sysObj = new TNamed(fSysString->GetName(),
 863                                 AliForwardUtil::CollisionSystemString(sys));
 864     sysObj->SetUniqueID(sys);
 865     fOutput->Add(sysObj); // fSysString->Clone());
 866   }
 867
 868   // Output centrality axis
 869   if (fCentAxis) fOutput->Add(fCentAxis);
 870
 871   // Output trigger string
 872   if (fTriggerString) {
 873     UShort_t mask = fTriggerString->GetUniqueID();
 874     TNamed* maskObj = new TNamed(fTriggerString->GetName(),
 875                                  AliAODForwardMult::GetTriggerString(mask));
 876     maskObj->SetUniqueID(mask);
 877     fOutput->Add(maskObj); // fTriggerString->Clone());
 878   }
 879
 880   // Normalization string
 881   if (fSchemeString) {
 882     UShort_t scheme = fSchemeString->GetUniqueID();
 883     TNamed* schemeObj = new TNamed(fSchemeString->GetName(),
 884                                    NormalizationSchemeString(scheme));
 885     schemeObj->SetUniqueID(scheme);
 886     fOutput->Add(schemeObj); // fSchemeString->Clone());
 887   }
 888
 889   // Output vertex axis
 890   TAxis* vtxAxis = new TAxis(1,fVtxMin,fVtxMax);
 891   vtxAxis->SetName("vtxAxis");
 892   vtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", fVtxMin,fVtxMax));
 893   fOutput->Add(vtxAxis);
 894
 895   // Output trigger efficiency and shape correction
 896   fOutput->Add(AliForwardUtil::MakeParameter("triggerEff", fTriggerEff));
 897   fOutput->Add(AliForwardUtil::MakeParameter("triggerEff0", fTriggerEff0));
 898   if (fShapeCorr) fOutput->Add(fShapeCorr);
 899
 900   TNamed* options = new TNamed("options","");
 901   TString str;
 902   str.Append(Form("Edges %scut, ", fCutEdges ? "" : "not "));
 903   str.Append(Form("Empty bins %scorrected for, ", fCorrEmpty ? "" : "not "));
 904   str.Append(Form("TH2::ProjectionX %sused", fUseROOTProj ? "" : "not "));
 905   options->SetTitle(str);
 906   fOutput->Add(options);
 907
 908   PostData(2, fOutput);
 909 }
 910 //________________________________________________________________________
 911 void
 912 AliBasedNdetaTask::LoadNormalizationData(UShort_t sys, UShort_t energy)
 913 {
 914   // Load the normalisation data for dN/deta for pp INEL, INEL>0, and NSD
 915   DGUARD(fDebug,1,"Load the normalization data for sys=%d, energy=%d",
 916          sys, energy);
 917   TString type("pp");
 918   TString snn("900");
 919   if(energy == 7000) snn.Form("7000");
 920   if(energy == 2750) snn.Form("2750");
 921
 922   // Check if shape correction/trigger efficiency was requsted and not
 923   // already set
 924   Bool_t needShape = ((fNormalizationScheme & kShape) && !fShapeCorr);
 925   Bool_t needEff   = ((fNormalizationScheme & kTriggerEfficiency) &&
 926                       ((1 - fTriggerEff) < 1e-6) && fTriggerEff > 0);
 927   if (needShape) AliInfo("Will load shape correction");
 928   if (needEff)   AliInfoF("Will load trigger efficiency, was=%f, %f",
 929                           fTriggerEff, fTriggerEff0);
 930   if(!needShape && !needShape) {
 931     AliInfo("Objects already set for normalization - no action taken");
 932     return;
 933   }
 934
 935   TString fname(Form("$ALICE_ROOT/PWGLF/FORWARD/corrections/"
 936                      "Normalization/normalizationHists_%s_%s.root",
 937                      type.Data(),snn.Data()));
 938   AliWarningF("Using old-style corrections from %s", fname.Data());
 939   TFile* fin = TFile::Open(fname, "READ");
 940   if(!fin) {
 941     AliWarningF("no file for normalization of %d/%d (%s)",
 942                 sys, energy, fname.Data());
 943     return;
 944   }
 945
 946   // Shape correction
 947   if (needShape) {
 948     TString trigName("All");
 949     if (fTriggerMask == AliAODForwardMult::kInel ||
 950         fTriggerMask == AliAODForwardMult::kNClusterGt0)
 951       trigName = "Inel";
 952     else if (fTriggerMask == AliAODForwardMult::kNSD)
 953       trigName = "NSD";
 954     else if (fTriggerMask == AliAODForwardMult::kInelGt0)
 955       trigName = "InelGt0";
 956     else {
 957       AliWarning(Form("Normalization for trigger %s not known, using all",
 958                       AliAODForwardMult::GetTriggerString(fTriggerMask)));
 959     }
 960
 961     TString shapeCorName(Form("h%sNormalization", trigName.Data()));
 962     TH2F*   shapeCor = dynamic_cast<TH2F*>(fin->Get(shapeCorName));
 963     if (shapeCor) SetShapeCorrection(shapeCor);
 964     else {
 965       AliWarning(Form("No shape correction found for %s", trigName.Data()));
 966     }
 967   }
 968
 969   // Trigger efficiency
 970   if (needEff) {
 971     TString effName(Form("%sTriggerEff",
 972                          fTriggerMask == AliAODForwardMult::kInel ? "inel" :
 973                          fTriggerMask == AliAODForwardMult::kNSD ? "nsd" :
 974                          fTriggerMask == AliAODForwardMult::kInelGt0 ?
 975                          "inelgt0" : "all"));
 976     Double_t trigEff = 1;
 977     TObject* eff = fin->Get(effName);
 978     if (eff) AliForwardUtil::GetParameter(eff, trigEff);
 979
 980     if (trigEff <= 0)
 981       AliWarningF("Retrieved trigger efficiency %s is %f<=0, ignoring",
 982                   effName.Data(), trigEff);
 983     else
 984       SetTriggerEff(trigEff);
 985
 986     // Trigger efficiency
 987     TString eff0Name(effName);
 988     eff0Name.Append("0");
 989
 990     Double_t trigEff0 = 1;
 991     TObject* eff0 = fin->Get(eff0Name);
 992     if (eff0) AliForwardUtil::GetParameter(eff, trigEff0);
 993     if (trigEff0 < 0)
 994       AliWarningF("Retrieved trigger efficiency %s is %f<0, ignoring",
 995                   eff0Name.Data(), trigEff0);
 996     else
 997       SetTriggerEff0(trigEff0);
 998   }
 999
1000   // TEMPORARY FIX
1001   // Rescale the shape correction by the trigger efficiency
1002   if (fShapeCorr) {
1003     AliWarning(Form("Rescaling shape correction by trigger efficiency: "
1004                     "1/E_X=1/%f", fTriggerEff));
1005     fShapeCorr->Scale(1. / fTriggerEff);
1006   }
1007   if (fin) fin->Close();
1008
1009   // Print - out
1010   if (fShapeCorr && fTriggerEff) AliInfo("Loaded objects for normalization.");
1011 }
1012
1013
1014 //________________________________________________________________________
1015 void
1016 AliBasedNdetaTask::Print(Option_t*) const
1017 {
1018   //
1019   // Print information
1020   //
1021   TString trigString("none");
1022   TString schemeString("none");
1023   TString sysString("unknown");
1024   TString sNNString("unknown");
1025   if (fTriggerString)
1026     trigString = AliAODForwardMult::GetTriggerString(fTriggerString->
1027                                                      GetUniqueID());
1028   if (fSchemeString)
1029     schemeString = NormalizationSchemeString(fSchemeString->GetUniqueID());
1030   if (fSysString)
1031     sysString = AliForwardUtil::CollisionSystemString(fSysString->
1032                                                       GetUniqueID());
1033   if (fSNNString)
1034     sNNString = AliForwardUtil::CenterOfMassEnergyString(fSNNString->
1035                                                          GetUniqueID());
1036
1037
1038   std::cout << this->ClassName() << ": " << this->GetName() << "\n"
1039             << std::boolalpha
1040             << " Trigger:                    " << trigString << "\n"
1041             << " Vertex range:               [" << fVtxMin << ":"
1042             << fVtxMax << "]\n"
1043             << " Rebin factor:               " << fRebin << "\n"
1044             << " Cut edges:                  " << fCutEdges << "\n"
1045             << " Symmertrice:                " << fSymmetrice << "\n"
1046             << " Use TH2::ProjectionX:       " << fUseROOTProj << "\n"
1047             << " Correct for empty:          " << fCorrEmpty << "\n"
1048             << " Normalization scheme:       " << schemeString <<"\n"
1049             << " Trigger efficiency:         " << fTriggerEff << "\n"
1050             << " Bin-0 Trigger efficiency:   " << fTriggerEff0 << "\n"
1051             << " Shape correction:           " << (fShapeCorr ?
1052                                                    fShapeCorr->GetName() :
1053                                                    "none") << "\n"
1054             << " sqrt(s_NN):                 " << sNNString << "\n"
1055             << " Collision system:           " << sysString << "\n"
1056             << " Centrality bins:            " << (fCentAxis ? "" : "none");
1057   if (fCentAxis) {
1058     Int_t           nBins = fCentAxis->GetNbins();
1059     const Double_t* bins  = fCentAxis->GetXbins()->GetArray();
1060     for (Int_t i = 0; i <= nBins; i++)
1061       std::cout << (i==0 ? "" : "-") << bins[i];
1062   }
1063   std::cout << std::noboolalpha << std::endl;
1064
1065 }
1066
1067 //________________________________________________________________________
1068 TH1D*
1069 AliBasedNdetaTask::Rebin(const TH1D* h, Int_t rebin, Bool_t cutEdges)
1070 {
1071   //
1072   // Make a copy of the input histogram and rebin that histogram
1073   //
1074   // Parameters:
1075   //    h  Histogram to rebin
1076   //
1077   // Return:
1078   //    New (rebinned) histogram
1079   //
1080   if (rebin <= 1) return 0;
1081
1082   Int_t nBins = h->GetNbinsX();
1083   if(nBins % rebin != 0) {
1084     AliWarningGeneral("AliBasedNdetaTask",
1085                       Form("Rebin factor %d is not a devisor of current number "
1086                            "of bins %d in the histogram %s",
1087                            rebin, nBins, h->GetName()));
1088     return 0;
1089   }
1090
1091   // Make a copy
1092   TH1D* tmp = static_cast<TH1D*>(h->Clone(Form("%s_rebin%02d",
1093                                                h->GetName(), rebin)));
1094   tmp->Rebin(rebin);
1095   // Hist should be reset, as it otherwise messes with the cutEdges option
1096   tmp->Reset();
1097   tmp->SetDirectory(0);
1098
1099   // The new number of bins
1100   Int_t nBinsNew = nBins / rebin;
1101   for(Int_t i = 1;i<= nBinsNew; i++) {
1102     Double_t content = 0;
1103     Double_t sumw    = 0;
1104     Double_t wsum    = 0;
1105     Int_t    nbins   = 0;
1106     for(Int_t j = 1; j<=rebin;j++) {
1107       Int_t    bin = (i-1)*rebin + j;
1108       Double_t c   =  h->GetBinContent(bin);
1109       if (c <= 0)  {
1110         continue; // old TODO: check
1111         //content = -1; // new
1112         //break; // also new
1113       }
1114
1115       if (cutEdges) {
1116         if (h->GetBinContent(bin+1)<=0 ||
1117             h->GetBinContent(bin-1)<=0) {
1118 #if 0
1119           AliWarningGeneral("AliBasedNdetaTask",
1120                             Form("removing bin %d=%f of %s (%d=%f,%d=%f)",
1121                                  bin, c, h->GetName(),
1122                                  bin+1, h->GetBinContent(bin+1),
1123                                  bin-1, h->GetBinContent(bin-1)));
1124 #endif
1125           continue;
1126         }
1127       }
1128       Double_t e =  h->GetBinError(bin);
1129       Double_t w =  1 / (e*e); // 1/c/c
1130       content    += c;
1131       sumw       += w;
1132       wsum       += w * c;
1133       nbins++;
1134     }
1135
1136     if(content > 0 && nbins > 0) {
1137       tmp->SetBinContent(i, wsum / sumw);
1138       tmp->SetBinError(i,1./TMath::Sqrt(sumw));
1139     }
1140   }
1141
1142   return tmp;
1143 }
1144
1145 //__________________________________________________________________
1146 TH1*
1147 AliBasedNdetaTask::Symmetrice(const TH1* h)
1148 {
1149   //
1150   // Make an extension of @a h to make it symmetric about 0
1151   //
1152   // Parameters:
1153   //    h Histogram to symmertrice
1154   //
1155   // Return:
1156   //    Symmetric extension of @a h
1157   //
1158   Int_t nBins = h->GetNbinsX();
1159   TH1* s = static_cast<TH1*>(h->Clone(Form("%s_mirror", h->GetName())));
1160   s->SetTitle(Form("%s (mirrored)", h->GetTitle()));
1161   s->Reset();
1162   s->SetBins(nBins, -h->GetXaxis()->GetXmax(), -h->GetXaxis()->GetXmin());
1163   s->SetMarkerColor(h->GetMarkerColor());
1164   s->SetMarkerSize(h->GetMarkerSize());
1165   s->SetMarkerStyle(FlipHollowStyle(h->GetMarkerStyle()));
1166   s->SetFillColor(h->GetFillColor());
1167   s->SetFillStyle(h->GetFillStyle());
1168   s->SetDirectory(0);
1169
1170   // Find the first and last bin with data
1171   Int_t first = nBins+1;
1172   Int_t last  = 0;
1173   for (Int_t i = 1; i <= nBins; i++) {
1174     if (h->GetBinContent(i) <= 0) continue;
1175     first = TMath::Min(first, i);
1176     last  = TMath::Max(last,  i);
1177   }
1178
1179   Double_t xfirst = h->GetBinCenter(first-1);
1180   Int_t    f1     = h->GetXaxis()->FindBin(-xfirst);
1181   Int_t    l2     = s->GetXaxis()->FindBin(xfirst);
1182   for (Int_t i = f1, j=l2; i <= last; i++,j--) {
1183     s->SetBinContent(j, h->GetBinContent(i));
1184     s->SetBinError(j, h->GetBinError(i));
1185   }
1186   // Fill in overlap bin
1187   s->SetBinContent(l2+1, h->GetBinContent(first));
1188   s->SetBinError(l2+1, h->GetBinError(first));
1189   return s;
1190 }
1191
1192 //__________________________________________________________________
1193 Int_t
1194 AliBasedNdetaTask::GetMarkerStyle(UShort_t bits)
1195 {
1196   Int_t  base   = bits & (0xFE);
1197   Bool_t hollow = bits & kHollow;
1198   switch (base) {
1199   case kCircle:       return (hollow ? 24 : 20);
1200   case kSquare:       return (hollow ? 25 : 21);
1201   case kUpTriangle:   return (hollow ? 26 : 22);
1202   case kDownTriangle: return (hollow ? 32 : 23);
1203   case kDiamond:      return (hollow ? 27 : 33);
1204   case kCross:        return (hollow ? 28 : 34);
1205   case kStar:         return (hollow ? 30 : 29);
1206   }
1207   return 1;
1208 }
1209 //__________________________________________________________________
1210 UShort_t
1211 AliBasedNdetaTask::GetMarkerBits(Int_t style)
1212 {
1213   UShort_t bits = 0;
1214   switch (style) {
1215   case 24: case 25: case 26: case 27: case 28: case 30: case 32:
1216     bits |= kHollow; break;
1217   }
1218   switch (style) {
1219   case 20: case 24: bits |= kCircle;       break;
1220   case 21: case 25: bits |= kSquare;       break;
1221   case 22: case 26: bits |= kUpTriangle;   break;
1222   case 23: case 32: bits |= kDownTriangle; break;
1223   case 27: case 33: bits |= kDiamond;      break;
1224   case 28: case 34: bits |= kCross;        break;
1225   case 29: case 30: bits |= kStar;         break;
1226   }
1227   return bits;
1228 }
1229 //__________________________________________________________________
1230 Int_t
1231 AliBasedNdetaTask::FlipHollowStyle(Int_t style)
1232 {
1233   UShort_t bits = GetMarkerBits(style);
1234   Int_t    ret  = GetMarkerStyle(bits ^ kHollow);
1235   return ret;
1236 }
1237
1238 //====================================================================
1239 void
1240 AliBasedNdetaTask::Sum::Init(TList* list, const TH2D* data, Int_t col)
1241 {
1242   DGUARD(fDebug,1,"Initializing sums with %s", data->GetName());
1243   TString n(GetHistName(0));
1244   TString n0(GetHistName(1));
1245   const char* postfix = GetTitle();
1246
1247   fSum = static_cast<TH2D*>(data->Clone(n));
1248   if (postfix) fSum->SetTitle(Form("%s (%s)", data->GetTitle(), postfix));
1249   fSum->SetDirectory(0);
1250   fSum->SetMarkerColor(col);
1251   fSum->SetMarkerStyle(GetMarkerStyle(kCircle|kSolid));
1252   fSum->Reset();
1253   list->Add(fSum);
1254
1255   fSum0 = static_cast<TH2D*>(data->Clone(n0));
1256   if (postfix)
1257     fSum0->SetTitle(Form("%s 0-bin (%s)", data->GetTitle(), postfix));
1258   else
1259     fSum0->SetTitle(Form("%s 0-bin", data->GetTitle()));
1260   fSum0->SetDirectory(0);
1261   fSum0->SetMarkerColor(col);
1262   fSum0->SetMarkerStyle(GetMarkerStyle(kCross|kHollow));
1263   fSum0->Reset();
1264   list->Add(fSum0);
1265
1266   fEvents = new TH1I(GetHistName(2), "Event types", 2, -.5, 1.5);
1267   fEvents->SetDirectory(0);
1268   fEvents->GetXaxis()->SetBinLabel(1, "Non-zero");
1269   fEvents->GetXaxis()->SetBinLabel(2, "Zero");
1270   list->Add(fEvents);
1271 }
1272
1273 //____________________________________________________________________
1274 TString
1275 AliBasedNdetaTask::Sum::GetHistName(const char* name, Int_t what, const char* post)
1276 {
1277   TString n(name);
1278   if      (what == 1) n.Append("0");
1279   else if (what == 2) n.Append("Events");
1280   if (post && post[0] != '\0')  n.Append(post);
1281   return n;
1282 }
1283
1284 //____________________________________________________________________
1285 TString
1286 AliBasedNdetaTask::Sum::GetHistName(Int_t what) const
1287 {
1288   return GetHistName(GetName(), what, GetTitle());
1289 }
1290
1291 //____________________________________________________________________
1292 void
1293 AliBasedNdetaTask::Sum::Add(const TH2D* data, Bool_t isZero)
1294 {
1295   DGUARD(fDebug,2,"Adding %s to sums", data->GetName());
1296   if (isZero) fSum0->Add(data);
1297   else        fSum->Add(data);
1298   fEvents->Fill(isZero ? 1 : 0);
1299 }
1300
1301 //____________________________________________________________________
1302 TH2D*
1303 AliBasedNdetaTask::Sum::CalcSum(TList*       output,
1304                                 Double_t&    ntotal,
1305                                 Double_t     epsilon0,
1306                                 Double_t     epsilon,
1307                                 Int_t        marker,
1308                                 Bool_t       rootProj,
1309                                 Bool_t       corrEmpty) const
1310 {
1311   DGUARD(fDebug,2,"Calculating final summed histogram %s", fSum->GetName());
1312
1313   // The return value `ret' is not scaled in anyway
1314   TH2D* ret      = static_cast<TH2D*>(fSum->Clone(fSum->GetName()));
1315   ret->SetDirectory(0);
1316   ret->Reset();
1317   Int_t n        = Int_t(fEvents->GetBinContent(1));
1318   Int_t n0       = Int_t(fEvents->GetBinContent(2));
1319
1320   AliInfoF("Adding histograms %s(%d) and %s(%d) with weights %f and %f resp.",
1321            fSum0->GetName(), n, fSum->GetName(), n0, 1./epsilon, 1./epsilon0);
1322   DMSG(fDebug,2, "Adding histograms %s and %s with weights %f and %f resp.",
1323        fSum0->GetName(), fSum->GetName(), 1./epsilon, 1./epsilon0);
1324   // Generate merged histogram
1325   ret->Add(fSum0, fSum, 1. / epsilon0, 1. / epsilon);
1326   ntotal = n / epsilon + n0 / epsilon0;
1327
1328   TList* out = new TList;
1329   out->SetOwner();
1330   const char* postfix = GetTitle();
1331   if (!postfix) postfix = "";
1332   out->SetName(Form("partial%s", postfix));
1333   output->Add(out);
1334
1335   // Now make copies, normalize them, and store in output list
1336   // Note, these are the only ones normalized here
1337   // These are mainly for diagnostics
1338   TH2D* sumCopy  = static_cast<TH2D*>(fSum->Clone("sum"));
1339   TH2D* sum0Copy = static_cast<TH2D*>(fSum0->Clone("sum0"));
1340   TH2D* retCopy  = static_cast<TH2D*>(ret->Clone("sumAll"));
1341   sumCopy->SetMarkerStyle(FlipHollowStyle(marker));
1342   sumCopy->SetDirectory(0);
1343   sum0Copy->SetMarkerStyle(GetMarkerStyle(GetMarkerBits(marker)+4));
1344   sum0Copy->SetDirectory(0);
1345   retCopy->SetMarkerStyle(marker);
1346   retCopy->SetDirectory(0);
1347
1348   Int_t nY      = fSum->GetNbinsY();
1349   Int_t o       = 0; // nY+1;
1350   TH1D* norm    = ProjectX(fSum,  "norm",    o, o, rootProj, corrEmpty, false);
1351   TH1D* norm0   = ProjectX(fSum0, "norm0",   o, o, rootProj, corrEmpty, false);
1352   TH1D* normAll = ProjectX(ret,   "normAll", o, o, rootProj, corrEmpty, false);
1353   norm->SetTitle("#eta coverage - >0-bin");
1354   norm0->SetTitle("#eta coverage - 0-bin");
1355   normAll->SetTitle("#eta coverage");
1356   norm->SetDirectory(0);
1357   norm0->SetDirectory(0);
1358   normAll->SetDirectory(0);
1359
1360   TH1D* sumCopyPx  = ProjectX(sumCopy,  "average",    1, nY,rootProj,corrEmpty);
1361   TH1D* sum0CopyPx = ProjectX(sum0Copy, "average0",   1, nY,rootProj,corrEmpty);
1362   TH1D* retCopyPx  = ProjectX(retCopy,  "averageAll", 1, nY,rootProj,corrEmpty);
1363   sumCopyPx->SetTitle(Form("#sum_{i}^{N_{#phi}}%s", sumCopy->GetTitle()));
1364   sum0CopyPx->SetTitle(Form("#sum_{i}^{N_{#phi}}%s", sum0Copy->GetTitle()));
1365   retCopyPx->SetTitle(Form("#sum_{i}^{N_{#phi}}%s", retCopy->GetTitle()));
1366   sumCopyPx->SetDirectory(0);
1367   sum0CopyPx->SetDirectory(0);
1368   retCopyPx->SetDirectory(0);
1369
1370   TH1D* phi    = ProjectX(fSum,  "phi",    nY+1, nY+1,rootProj,corrEmpty,false);
1371   TH1D* phi0   = ProjectX(fSum0, "phi0",   nY+1, nY+1,rootProj,corrEmpty,false);
1372   TH1D* phiAll = ProjectX(ret,   "phiAll", nY+1, nY+1,rootProj,corrEmpty,false);
1373   phi->SetTitle("#phi acceptance from dead strips - >0-bin");
1374   phi0->SetTitle("#phi acceptance from dead strips - 0-bin");
1375   phiAll->SetTitle("#phi acceptance from dead strips");
1376   phi->SetDirectory(0);
1377   phi0->SetDirectory(0);
1378   phiAll->SetDirectory(0);
1379
1380   const TH1D* cov    = (corrEmpty ? norm    : phi);
1381   const TH1D* cov0   = (corrEmpty ? norm0   : phi0);
1382   const TH1D* covAll = (corrEmpty ? normAll : phiAll);
1383
1384   // Here, we scale to the coverage (or phi acceptance)
1385   ScaleToCoverage(sumCopy,  cov);
1386   ScaleToCoverage(sum0Copy, cov0);
1387   ScaleToCoverage(retCopy,  covAll);
1388
1389   // Scale our 1D histograms
1390   sumCopyPx->Scale(1., "width");
1391   sum0CopyPx->Scale(1., "width");
1392   retCopyPx->Scale(1., "width");
1393
1394   AliInfo(Form("Maximum %f,%f changed to %f", sumCopyPx->GetMaximum(),
1395                sum0CopyPx->GetMaximum(), retCopyPx->GetMaximum()));
1396
1397   // Scale the normalization - they should be 1 at the maximum
1398   norm->Scale(n > 0   ? 1. / n  : 1);
1399   norm0->Scale(n0 > 0 ? 1. / n0 : 1);
1400   normAll->Scale(ntotal > 0 ? 1. / ntotal : 1);
1401
1402   // Scale the normalization - they should be 1 at the maximum
1403   phi->Scale(n > 0   ? 1. / n  : 1);
1404   phi0->Scale(n0 > 0 ? 1. / n0 : 1);
1405   phiAll->Scale(ntotal > 0 ? 1. / ntotal : 1);
1406
1407   out->Add(sumCopy);
1408   out->Add(sum0Copy);
1409   out->Add(retCopy);
1410   out->Add(sumCopyPx);
1411   out->Add(sum0CopyPx);
1412   out->Add(retCopyPx);
1413   out->Add(norm);
1414   out->Add(norm0);
1415   out->Add(normAll);
1416   out->Add(phi);
1417   out->Add(phi0);
1418   out->Add(phiAll);
1419
1420   AliInfoF("Returning  (1/%f * %s + 1/%f * %s), "
1421            "1/%f * %d + 1/%f * %d = %d",
1422            epsilon0, fSum0->GetName(), epsilon, fSum->GetName(),
1423            epsilon0, n0, epsilon, n, int(ntotal));
1424 #if 0
1425   for (Int_t i = 1; i <= ret->GetNbinsX(); i++) {
1426     Double_t nc  = sum->GetBinContent(i, 0);
1427     Double_t nc0 = sum0->GetBinContent(i, 0);
1428     ret->SetBinContent(i, 0, nc + nc0); // Just count events
1429   }
1430 #endif
1431
1432   return ret;
1433 }
1434
1435 //====================================================================
1436 AliBasedNdetaTask::CentralityBin::CentralityBin()
1437   : TNamed("", ""),
1438     fSums(0),
1439     fOutput(0),
1440     fSum(0),
1441     fSumMC(0),
1442     fTriggers(0),
1443     fStatus(0),
1444     fLow(0),
1445     fHigh(0),
1446     fDoFinalMCCorrection(false),
1447     fSatelliteVertices(false),
1448     fDebug(0)
1449 {
1450   //
1451   // Constructor
1452   //
1453   DGUARD(fDebug,3,"Default CTOR of AliBasedNdeta::CentralityBin");
1454 }
1455 //____________________________________________________________________
1456 AliBasedNdetaTask::CentralityBin::CentralityBin(const char* name,
1457                                                 Short_t low, Short_t high)
1458   : TNamed(name, ""),
1459     fSums(0),
1460     fOutput(0),
1461     fSum(0),
1462     fSumMC(0),
1463     fTriggers(0),
1464     fStatus(0),
1465     fLow(low),
1466     fHigh(high),
1467     fDoFinalMCCorrection(false),
1468     fSatelliteVertices(false),
1469     fDebug(0)
1470 {
1471   //
1472   // Constructor
1473   //
1474   // Parameters:
1475   //    name Name used for histograms (e.g., Forward)
1476   //    low  Lower centrality cut in percent
1477   //    high Upper centrality cut in percent
1478   //
1479   DGUARD(fDebug,3,"Named CTOR of AliBasedNdeta::CentralityBin: %s [%3d,%3d]",
1480          name,low,high);
1481   if (low <= 0 && high <= 0) {
1482     fLow  = 0;
1483     fHigh = 0;
1484     SetTitle("All centralities");
1485   }
1486   else {
1487     fLow  = low;
1488     fHigh = high;
1489     SetTitle(Form("Centrality bin from %3d%% to %3d%%", low, high));
1490   }
1491 }
1492 //____________________________________________________________________
1493 AliBasedNdetaTask::CentralityBin::CentralityBin(const CentralityBin& o)
1494   : TNamed(o),
1495     fSums(o.fSums),
1496     fOutput(o.fOutput),
1497     fSum(o.fSum),
1498     fSumMC(o.fSumMC),
1499     fTriggers(o.fTriggers),
1500     fStatus(o.fStatus),
1501     fLow(o.fLow),
1502     fHigh(o.fHigh),
1503     fDoFinalMCCorrection(o.fDoFinalMCCorrection),
1504     fSatelliteVertices(o.fSatelliteVertices),
1505     fDebug(o.fDebug)
1506 {
1507   //
1508   // Copy constructor
1509   //
1510   // Parameters:
1511   //    other Object to copy from
1512   //
1513   DGUARD(fDebug,3,"Copy CTOR of AliBasedNdeta::CentralityBin");
1514 }
1515 //____________________________________________________________________
1516 AliBasedNdetaTask::CentralityBin::~CentralityBin()
1517 {
1518   //
1519   // Destructor
1520   //
1521   DGUARD(fDebug,3,"DTOR of AliBasedNdeta::CentralityBin");
1522
1523   // if (fSums) fSums->Delete();
1524   // if (fOutput) fOutput->Delete();
1525 }
1526
1527 //____________________________________________________________________
1528 AliBasedNdetaTask::CentralityBin&
1529 AliBasedNdetaTask::CentralityBin::operator=(const CentralityBin& o)
1530 {
1531   //
1532   // Assignment operator
1533   //
1534   // Parameters:
1535   //    other Object to assign from
1536   //
1537   // Return:
1538   //    Reference to this
1539   //
1540   DGUARD(fDebug,3,"Centrality bin assignment");
1541   if (&o == this) return *this;
1542   SetName(o.GetName());
1543   SetTitle(o.GetTitle());
1544   fSums      = o.fSums;
1545   fOutput    = o.fOutput;
1546   fSum       = o.fSum;
1547   fSumMC     = o.fSumMC;
1548   fTriggers  = o.fTriggers;
1549   fStatus    = o.fStatus;
1550   fLow       = o.fLow;
1551   fHigh      = o.fHigh;
1552   fDoFinalMCCorrection = o.fDoFinalMCCorrection;
1553   fSatelliteVertices = o.fSatelliteVertices;
1554
1555   return *this;
1556 }
1557 //____________________________________________________________________
1558 Int_t
1559 AliBasedNdetaTask::CentralityBin::GetColor(Int_t fallback) const
1560 {
1561   if (IsAllBin()) return fallback;
1562   Float_t  fc       = (fLow+double(fHigh-fLow)/2) / 100;
1563   Int_t    nCol     = gStyle->GetNumberOfColors();
1564   Int_t    icol     = TMath::Min(nCol-1,int(fc * nCol + .5));
1565   Int_t    col      = gStyle->GetColorPalette(icol);
1566   return col;
1567 }
1568 //____________________________________________________________________
1569 const char*
1570 AliBasedNdetaTask::CentralityBin::GetListName() const
1571 {
1572   //
1573   // Get the list name
1574   //
1575   // Return:
1576   //    List Name
1577   //
1578   if (IsAllBin()) return "all";
1579   return Form("cent%03d_%03d", fLow, fHigh);
1580 }
1581 //____________________________________________________________________
1582 void
1583 AliBasedNdetaTask::CentralityBin::CreateOutputObjects(TList* dir, Int_t mask)
1584 {
1585   //
1586   // Create output objects
1587   //
1588   // Parameters:
1589   //    dir   Parent list
1590   //
1591   DGUARD(fDebug,1,"Create centrality bin output objects");
1592   fSums = new TList;
1593   fSums->SetName(GetListName());
1594   fSums->SetOwner();
1595   dir->Add(fSums);
1596
1597   fTriggers = AliAODForwardMult::MakeTriggerHistogram("triggers", mask);
1598   fTriggers->SetDirectory(0);
1599
1600   fStatus = AliAODForwardMult::MakeStatusHistogram("status");
1601   fStatus->SetDirectory(0);
1602
1603   fSums->Add(fTriggers);
1604   fSums->Add(fStatus);
1605 }
1606 //____________________________________________________________________
1607 void
1608 AliBasedNdetaTask::CentralityBin::SetDebugLevel(Int_t lvl)
1609 {
1610   fDebug = lvl;
1611   if (fSum) fSum->fDebug = lvl;
1612   if (fSumMC) fSumMC->fDebug = lvl;
1613 }
1614
1615 //____________________________________________________________________
1616 Bool_t
1617 AliBasedNdetaTask::CentralityBin::ReadSum(TList* list, bool mc)
1618 {
1619   const char* post = (mc ? "MC" : "");
1620   TString     sn   = Sum::GetHistName(GetName(),0,post);
1621   TString     sn0  = Sum::GetHistName(GetName(),1,post);
1622   TString     ev   = Sum::GetHistName(GetName(),2,post);
1623   TH2D* sum      = static_cast<TH2D*>(list->FindObject(sn));
1624   TH2D* sum0     = static_cast<TH2D*>(list->FindObject(sn0));
1625   TH1I* events   = static_cast<TH1I*>(list->FindObject(ev));
1626   if (!sum || !sum0 || !events) {
1627     AliWarningF("Failed to find one or more histograms: "
1628                 "%s (%p) %s (%p) %s (%p)",
1629                 sn.Data(), sum,
1630                 sn0.Data(), sum0,
1631                 ev.Data(), events);
1632     return false;
1633   }
1634   Sum* ret     = new Sum(GetName(), post);
1635   ret->fSum    = sum;
1636   ret->fSum0   = sum0;
1637   ret->fEvents = events;
1638   ret->fDebug  = fDebug;
1639   if (mc) fSumMC = ret;
1640   else    fSum   = ret;
1641
1642   return true;
1643 }
1644
1645 //____________________________________________________________________
1646 void
1647 AliBasedNdetaTask::CentralityBin::CreateSums(const TH2D* data, const TH2D* mc)
1648 {
1649   //
1650   // Create sum histogram
1651   //
1652   // Parameters:
1653   //    data  Data histogram to clone
1654   //    mc    (optional) MC histogram to clone
1655   //
1656   DGUARD(fDebug,1,"Create centrality bin sums from %s",
1657          data ? data->GetName() : "(null)");
1658   if (data) {
1659     fSum = new Sum(GetName(),"");
1660     fSum->Init(fSums, data, GetColor());
1661     fSum->fDebug = fDebug;
1662   }
1663
1664   // If no MC data is given, then do not create MC sum histogram
1665   if (!mc) return;
1666
1667   fSumMC = new Sum(GetName(), "MC");
1668   fSumMC->Init(fSums, mc, GetColor());
1669   fSumMC->fDebug = fDebug;
1670 }
1671
1672 //____________________________________________________________________
1673 Bool_t
1674 AliBasedNdetaTask::CentralityBin::CheckEvent(const AliAODForwardMult* forward,
1675                                              Int_t triggerMask,
1676                                              Double_t vzMin, Double_t vzMax)
1677 {
1678   //
1679   // Check the trigger, vertex, and centrality
1680   //
1681   // Parameters:
1682   //    aod Event input
1683   //
1684   // Return:
1685   //    true if the event is to be used
1686   //
1687   if (!forward) return false;
1688
1689   DGUARD(fDebug,2,"Check the event");
1690   // We do not check for centrality here - it's already done
1691   return forward->CheckEvent(triggerMask, vzMin, vzMax, 0, 0,
1692                              fTriggers, fStatus);
1693 }
1694
1695
1696 //____________________________________________________________________
1697 Bool_t
1698 AliBasedNdetaTask::CentralityBin::ProcessEvent(const AliAODForwardMult* forward,
1699                                                Int_t triggerMask, Bool_t isZero,
1700                                                Double_t vzMin, Double_t vzMax,
1701                                                const TH2D* data, const TH2D* mc)
1702 {
1703   //
1704   // Process an event
1705   //
1706   // Parameters:
1707   //    forward     Forward data (for trigger, vertex, & centrality)
1708   //    triggerMask Trigger mask
1709   //    vzMin       Minimum IP z coordinate
1710   //    vzMax       Maximum IP z coordinate
1711   //    data        Data histogram
1712   //    mc          MC histogram
1713   //
1714   DGUARD(fDebug,1,"Process one event for %s a given centrality bin",
1715          data ? data->GetName() : "(null)");
1716   if (!CheckEvent(forward, triggerMask, vzMin, vzMax)) return false;
1717   if (!data) return false;
1718   if (!fSum) CreateSums(data, mc);
1719
1720   fSum->Add(data, isZero);
1721   if (mc) fSumMC->Add(mc, isZero);
1722
1723   return true;
1724 }
1725
1726 //________________________________________________________________________
1727 Double_t
1728 AliBasedNdetaTask::CentralityBin::Normalization(const TH1I& t,
1729                                                 UShort_t    scheme,
1730                                                 Double_t    trigEff,
1731                                                 Double_t&   ntotal,
1732                                                 TString*    text) const
1733 {
1734   //
1735   // Calculate normalization
1736   //
1737   // Parameters:
1738   //    t       Trigger histogram
1739   //    scheme  Normaliztion scheme
1740   //    trigEff From MC
1741   //    ntotal  On return, contains the number of events.
1742   //
1743   DGUARD(fDebug,1,"Normalize centrality bin %s [%3d-%3d%%] with %s",
1744          GetName(), fLow, fHigh, t.GetName());
1745   Double_t nAll        = t.GetBinContent(AliAODForwardMult::kBinAll);
1746   Double_t nB          = t.GetBinContent(AliAODForwardMult::kBinB);
1747   Double_t nA          = t.GetBinContent(AliAODForwardMult::kBinA);
1748   Double_t nC          = t.GetBinContent(AliAODForwardMult::kBinC);
1749   Double_t nE          = t.GetBinContent(AliAODForwardMult::kBinE);
1750   Double_t nOffline    = t.GetBinContent(AliAODForwardMult::kBinOffline);
1751   Double_t nTriggered  = t.GetBinContent(AliAODForwardMult::kWithTrigger);
1752   Double_t nWithVertex = t.GetBinContent(AliAODForwardMult::kWithVertex);
1753   Double_t nAccepted   = ntotal;
1754   ntotal               = 0;
1755
1756   if (nTriggered <= 0.1) {
1757     AliError("Number of triggered events <= 0");
1758     return -1;
1759   }
1760   if (nWithVertex <= 0.1) {
1761     AliError("Number of events with vertex <= 0");
1762     return -1;
1763   }
1764   ntotal          = nAccepted;
1765   Double_t vtxEff = nWithVertex / nTriggered;
1766   Double_t scaler = 1;
1767   Double_t beta   = nA + nC - 2*nE;
1768
1769
1770   TString rhs("N = N_acc");
1771   if (!(scheme & kZeroBin)) {
1772     if (scheme & kEventLevel) {
1773       ntotal = nAccepted / vtxEff;
1774       scaler = vtxEff;
1775       AliInfoF("Calculating event normalisation as\n"
1776                " N = N_A * N_T / N_V = %d * %d / %d = %f (%f)",
1777                Int_t(nAccepted), Int_t(nTriggered), Int_t(nWithVertex),
1778                ntotal, scaler);
1779       if (scheme & kBackground) {
1780         //          1            E_V             E_V
1781         //   s = --------- = ------------- = ------------
1782         //        1 - beta   1 - beta E_V    1 - beta N_V
1783         //       ---  ----       --------        ---- ---
1784         //       E_V  N_V          N_V           N_V  N_T
1785         //
1786         //          E_V
1787         //     = ------------
1788         //        1 - beta
1789         //            ----
1790         //             N_T
1791         //
1792         ntotal -= nAccepted * beta / nWithVertex;
1793         // This one is direct and correct.
1794         // scaler = 1. / (1. / vtxEff - beta / nWithVertex);
1795         // A simpler expresion
1796         scaler /= (1 - beta / nTriggered); // 0.831631 -> 0.780689
1797         AliInfo(Form("Calculating event normalisation as\n"
1798                      " beta = N_a + N_c + 2 N_e = %d + %d - 2 * %d = %d\n"
1799                      " N = N - N_A * beta / N_V = %f - %d * %d / %d = %f (%f)",
1800                      Int_t(nA), Int_t(nC), Int_t(nE), Int_t(beta),
1801                      nAccepted / vtxEff, Int_t(nAccepted), Int_t(beta),
1802                      Int_t(nWithVertex), ntotal, scaler));
1803         rhs.Append("(1/eps_V - beta/N_vtx)");
1804       } // Background
1805       else
1806         rhs.Append("/eps_V");
1807     } // Event-level
1808     if (scheme & kTriggerEfficiency) {
1809       ntotal /= trigEff;
1810       scaler *= trigEff;
1811       AliInfo(Form("Correcting for trigger efficiency:\n"
1812                    " N = 1 / E_X * N = 1 / %f * %d = %f (%f)",
1813                    trigEff, Int_t(ntotal), ntotal / trigEff, scaler));
1814       rhs.Append("/eps_T");
1815     } // Trigger efficiency
1816   }
1817   else  {
1818     // Calculate as
1819     //
1820     //  N = N_A + 1/E_X * N_A / N_V (N_T - N_V - beta)
1821     //    = N_A (1 + 1/E_X (N_T/N_V - 1 - beta / N_V))
1822     //    = N_A (1 + 1/E_X (1/E_V - 1 - beta / N_V))
1823     //
1824     //  s = N_A/N = 1 / (1 + 1/E_X (N_T/N_V - 1 - beta / N_V))
1825     //    = N_V / (N_V + 1/E_X (N_T - N_V - beta))
1826     //
1827     if (!(scheme & kBackground)) beta = 0;
1828     ntotal = nAccepted * (1 + 1/trigEff * (nTriggered / nWithVertex - 1
1829                                          - beta / nWithVertex));
1830     scaler = nWithVertex / (nWithVertex +
1831                             1/trigEff * (nTriggered-nWithVertex-beta));
1832     AliInfo(Form("Calculating event normalisation as\n"
1833                  "  beta = N_a + N_c + 2 N_e = %d + %d - 2 * %d = %d\n"
1834                  "  N = N_A (1 + 1/E_X (N_T/N_V - 1 - beta / N_V)) = "
1835                  "%d (1 + 1 / %f (%d / %d - 1 - %d / %d)) = %f (%f)",
1836                  Int_t(nA), Int_t(nC), Int_t(nE), Int_t(beta),
1837                  Int_t(nAccepted), trigEff, Int_t(nTriggered),
1838                  Int_t(nWithVertex), Int_t(beta), Int_t(nWithVertex),
1839                  ntotal, scaler));
1840     rhs.Append("(1+1/eps_T(1/eps_V-1-beta/N_vtx))");
1841   }
1842
1843   if (text) {
1844     text->Append(Form("%-40s = %d\n", "N_all",            UInt_t(nAll)));
1845     text->Append(Form("%-40s = %d\n", "N_acc",            UInt_t(nAccepted)));
1846     text->Append(Form("%-40s = %d\n", "N_trg",            UInt_t(nTriggered)));
1847     text->Append(Form("%-40s = %d\n", "N_vtx",            UInt_t(nWithVertex)));
1848     text->Append(Form("%-40s = %d\n", "N_B",              UInt_t(nB)));
1849     text->Append(Form("%-40s = %d\n", "N_A",              UInt_t(nA)));
1850     text->Append(Form("%-40s = %d\n", "N_C",              UInt_t(nC)));
1851     text->Append(Form("%-40s = %d\n", "N_E",              UInt_t(nE)));
1852     text->Append(Form("%-40s = %d\n", "beta = N_A + N_C - 2N_E",UInt_t(beta)));
1853     text->Append(Form("%-40s = %f\n", "eps_V = N_vtx/N_trg",vtxEff));
1854     text->Append(Form("%-40s = %f\n", "eps_T",            trigEff));
1855     text->Append(Form("%-40s = %f\n", rhs.Data(),         ntotal));
1856   }
1857   AliInfo(Form("\n"
1858                " Total of        %9d events for %s\n"
1859                "  of these       %9d have an offline trigger\n"
1860                "  of these N_T = %9d has the selected trigger\n"
1861                "  of these N_V = %9d has a vertex\n"
1862                "  of these N_A = %9d were in the selected range\n"
1863                "  Triggers by hardware type:\n"
1864                "    N_b   =  %9d\n"
1865                "    N_ac  =  %9d (%d+%d)\n"
1866                "    N_e   =  %9d\n"
1867                "  Vertex efficiency:          %f\n"
1868                "  Trigger efficiency:         %f\n"
1869                "  Total number of events: N = %f\n"
1870                "  Scaler (N_A/N):             %f\n"
1871                "  %25s = %f",
1872                Int_t(nAll), GetTitle(), Int_t(nOffline),
1873                Int_t(nTriggered), Int_t(nWithVertex), Int_t(nAccepted),
1874                Int_t(nB), Int_t(nA+nC), Int_t(nA), Int_t(nC), Int_t(nE),
1875                vtxEff, trigEff, ntotal, scaler, rhs.Data(), ntotal));
1876   return scaler;
1877 }
1878
1879 //________________________________________________________________________
1880 const char*
1881 AliBasedNdetaTask::CentralityBin::GetResultName(Int_t rebin,
1882                                                 Bool_t sym,
1883                                                 const char* postfix) const
1884 {
1885   static TString n;
1886   n = Form("dndeta%s%s",GetName(), postfix);
1887   if (rebin > 1) n.Append(Form("_rebin%02d", rebin));
1888   if (sym)       n.Append("_mirror");
1889   return n.Data();
1890 }
1891 //________________________________________________________________________
1892 TH1*
1893 AliBasedNdetaTask::CentralityBin::GetResult(Int_t rebin,
1894                                             Bool_t sym,
1895                                             const char* postfix) const
1896 {
1897   if (!fOutput) {
1898     AliWarning(Form("No output list defined in %s [%3d,%3d]", GetName(),
1899                     fLow, fHigh));
1900     return 0;
1901   }
1902   TString  n = GetResultName(rebin, sym, postfix);
1903   TObject* o = fOutput->FindObject(n.Data());
1904   if (!o) {
1905     // AliWarning(Form("Object %s not found in output list", n.Data()));
1906     return 0;
1907   }
1908   return static_cast<TH1*>(o);
1909 }
1910
1911 //________________________________________________________________________
1912 void
1913 AliBasedNdetaTask::CentralityBin::MakeResult(const TH2D* sum,
1914                                              const char* postfix,
1915                                              bool        rootProj,
1916                                              bool        corrEmpty,
1917                                              const TH2F* shapeCorr,
1918                                              Double_t    scaler,
1919                                              bool        symmetrice,
1920                                              Int_t       rebin,
1921                                              bool        cutEdges,
1922                                              Int_t       marker,
1923                                              Int_t       color,
1924                                              TList*      mclist,
1925                                              TList*      truthlist)
1926 {
1927   //
1928   // Generate the dN/deta result from input
1929   //
1930   // Parameters:
1931   //     sum        Sum of 2D hists
1932   //     postfix    Post fix on names
1933   //     rootProj   Whether to use ROOT TH2::ProjectionX
1934   //     corrEmpty  Correct for empty bins
1935   //     shapeCorr  Shape correction to use
1936   //     scaler     Event-level normalization scaler
1937   //     symmetrice Whether to make symmetric extensions
1938   //     rebin      Whether to rebin
1939   //     cutEdges   Whether to cut edges when rebinning
1940   //
1941   DGUARD(fDebug,1,"Make centrality bin result from %s", sum->GetName());
1942   TH2D* copy    = static_cast<TH2D*>(sum->Clone(Form("d2Ndetadphi%s%s",
1943                                                      GetName(), postfix)));
1944
1945   TH1D* accNorm = 0;
1946   Int_t nY      = sum->GetNbinsY();
1947   // Hack HHD Hans test code to manually remove FMD2 dead channel (but
1948   // it is on outer?)
1949   //
1950   // cholm comment: The original hack has been moved to
1951   // AliForwarddNdetaTask::CheckEventData - this simplifies things a
1952   // great deal, and we could in principle use the new phi acceptance.
1953   //
1954   // However, since we may have filtered out the dead sectors in the
1955   // AOD production already, we can't be sure we can recalculate the
1956   // phi acceptance correctly, so for now, we rely on fCorrEmpty being set.
1957   Int_t o       = (corrEmpty ? 0 : nY+1);
1958   accNorm = ProjectX(sum, Form("norm%s%s",GetName(), postfix), o, o,
1959                      rootProj, corrEmpty, false);
1960   accNorm->SetDirectory(0);
1961
1962   // ---- Scale by shape correction ----------------------------------
1963   if (shapeCorr) copy->Divide(shapeCorr);
1964   else AliInfo("No shape correction specified, or disabled");
1965
1966   // --- Normalize to the coverage -----------------------------------
1967   if (corrEmpty) {
1968     ScaleToCoverage(copy, accNorm);
1969     // --- Event-level normalization ---------------------------------
1970     copy->Scale(scaler);
1971   }
1972
1973   // --- Project on X axis -------------------------------------------
1974   TH1D* dndeta = ProjectX(copy, Form("dndeta%s%s",GetName(), postfix),
1975                           1, nY, rootProj, corrEmpty);
1976   dndeta->SetDirectory(0);
1977   // Event-level normalization
1978   if (!corrEmpty) {
1979     ScaleToCoverage(dndeta, accNorm);
1980     dndeta->Scale(scaler);
1981   }
1982   dndeta->Scale(1., "width");
1983   copy->Scale(1., "width");
1984
1985   TH1D* dndetaMCCorrection = 0;
1986   TList* centlist          = 0;
1987   TH1D* dndetaMCtruth      = 0;
1988   TList* truthcentlist     = 0;
1989
1990   // --- Possible final correction to <MC analysis> / <MC truth> -----
1991   // we get the rebinned distribution for satellite to make the correction
1992   TString rebinSuf(fSatelliteVertices ? "_rebin05" : "");
1993   if(mclist)
1994     centlist = static_cast<TList*> (mclist->FindObject(GetListName()));
1995   if(centlist)
1996     dndetaMCCorrection =
1997       static_cast<TH1D*>(centlist->FindObject(Form("dndeta%s%s%s",
1998                                                    GetName(), postfix,
1999                                                    rebinSuf.Data())));
2000   if (truthlist)
2001     truthcentlist =
2002       static_cast<TList*>(truthlist->FindObject(GetListName()));
2003   if (truthcentlist)
2004     // TODO here new is "dndetaTruth"
2005     dndetaMCtruth =
2006       static_cast<TH1D*>(truthcentlist->FindObject(Form("dndetaMCTruth%s",
2007                                                         rebinSuf.Data())));
2008
2009   if (dndetaMCCorrection && dndetaMCtruth) {
2010     AliInfo("Correcting with final MC correction");
2011     TString testString(dndetaMCCorrection->GetName());
2012
2013     // We take the measured MC dN/deta and divide with truth
2014     dndetaMCCorrection->Divide(dndetaMCtruth);
2015     dndetaMCCorrection->SetTitle("Final MC correction");
2016     dndetaMCCorrection->SetName("finalMCCorr");
2017     for(Int_t m = 1; m <= dndetaMCCorrection->GetNbinsX(); m++) {
2018       if(dndetaMCCorrection->GetBinContent(m) < 0.5 ||
2019          dndetaMCCorrection->GetBinContent(m) > 1.75) {
2020         dndetaMCCorrection->SetBinContent(m,1.);
2021         dndetaMCCorrection->SetBinError(m,0.1);
2022       }
2023     }
2024     // Applying the correction
2025     if (!fSatelliteVertices)
2026       // For non-satellites we took the same binning, so we do a straight
2027       // division
2028       dndeta->Divide(dndetaMCCorrection);
2029     else {
2030       // For satelitte events, we took coarser binned histograms, so
2031       // we need to do a bit more
2032       for(Int_t m = 1; m <= dndeta->GetNbinsX(); m++) {
2033         if(dndeta->GetBinContent(m) <= 0.01 ) continue;
2034
2035         Double_t eta     = dndeta->GetXaxis()->GetBinCenter(m);
2036         Int_t    bin     = dndetaMCCorrection->GetXaxis()->FindBin(eta);
2037         Double_t mccorr  = dndetaMCCorrection->GetBinContent(bin);
2038         Double_t mcerror = dndetaMCCorrection->GetBinError(bin);
2039         if (mccorr < 1e-6) {
2040           dndeta->SetBinContent(m, 0);
2041           dndeta->SetBinError(m, 0);
2042         }
2043         Double_t value   = dndeta->GetBinContent(m);
2044         Double_t error   = dndeta->GetBinError(m);
2045         Double_t sumw2   = (error   * error   * mccorr * mccorr +
2046                             mcerror * mcerror * value  * value);
2047         dndeta->SetBinContent(m,value/mccorr) ;
2048         dndeta->SetBinError(m,TMath::Sqrt(sumw2)/mccorr/mccorr);
2049       }
2050     }
2051   }
2052   else
2053     AliInfo("No final MC correction applied");
2054
2055   // --- Set some histogram attributes -------------------------------
2056   TString post;
2057   Int_t rColor = GetColor(color);
2058   if (postfix && postfix[0] != '\0') post = Form(" (%s)", postfix);
2059   SetHistogramAttributes(dndeta,  rColor, marker,
2060                          Form("ALICE %s%s", GetName(), post.Data()));
2061   SetHistogramAttributes(accNorm, rColor, marker,
2062                          Form("ALICE %s normalisation%s",
2063                               GetName(), post.Data()));
2064
2065   // --- Make symmetric extensions and rebinnings --------------------
2066   if (symmetrice) fOutput->Add(Symmetrice(dndeta));
2067   fOutput->Add(dndeta);
2068   fOutput->Add(accNorm);
2069   fOutput->Add(copy);
2070   fOutput->Add(Rebin(dndeta, rebin, cutEdges));
2071   if (symmetrice) fOutput->Add(Symmetrice(Rebin(dndeta, rebin, cutEdges)));
2072   if (dndetaMCCorrection) fOutput->Add(dndetaMCCorrection);
2073
2074   // HHD Test of dN/deta in phi bins add flag later?
2075   //
2076   // cholm comment: We disable this for now
2077 #if 0
2078   for (Int_t nn=1; nn <= sum->GetNbinsY(); nn++) {
2079     TH1D* dndeta_phi = ProjectX(copy, Form("dndeta%s%s_phibin%d",
2080                                            GetName(), postfix, nn),
2081                                 nn, nn, rootProj, corrEmpty);
2082     dndeta_phi->SetDirectory(0);
2083     // Event-level normalization
2084     dndeta_phi->Scale(TMath::Pi()/10., "width");
2085
2086     if(centlist)
2087       dndetaMCCorrection =
2088         static_cast<TH1D*>(centlist->FindObject(Form("dndeta%s%s_phibin%d",
2089                                                      GetName(), postfix,nn)));
2090     if(truthcentlist)
2091       dndetaMCtruth
2092         = static_cast<TH1D*>(truthcentlist->FindObject("dndetaMCTruth"));
2093
2094     if (dndetaMCCorrection && dndetaMCtruth) {
2095       AliInfo("Correcting with final MC correction");
2096       TString testString(dndetaMCCorrection->GetName());
2097       dndetaMCCorrection->Divide(dndetaMCtruth);
2098       dndetaMCCorrection->SetTitle(Form("Final_MC_correction_phibin%d",nn));
2099       dndetaMCCorrection->SetName(Form("Final_MC_correction_phibin%d",nn));
2100       for(Int_t m = 1; m <= dndetaMCCorrection->GetNbinsX(); m++) {
2101         if(dndetaMCCorrection->GetBinContent(m) < 0.25 ||
2102            dndetaMCCorrection->GetBinContent(m) > 1.75) {
2103           dndetaMCCorrection->SetBinContent(m,1.);
2104           dndetaMCCorrection->SetBinError(m,0.1);
2105         }
2106       }
2107       //Applying the correction
2108       dndeta_phi->Divide(dndetaMCCorrection);
2109     }
2110     fOutput->Add(dndeta_phi);
2111     fOutput->Add(Rebin(dndeta_phi, rebin, cutEdges));
2112     if(dndetaMCCorrection) fOutput->Add(dndetaMCCorrection);
2113   } // End of phi
2114 #endif
2115 }
2116
2117 //________________________________________________________________________
2118 void
2119 AliBasedNdetaTask::CentralityBin::End(TList*      sums,
2120                                       TList*      results,
2121                                       UShort_t    scheme,
2122                                       const TH2F* shapeCorr,
2123                                       Double_t    trigEff,
2124                                       Double_t    trigEff0,
2125                                       Bool_t      symmetrice,
2126                                       Int_t       rebin,
2127                                       Bool_t      rootProj,
2128                                       Bool_t      corrEmpty,
2129                                       Bool_t      cutEdges,
2130                                       Int_t       triggerMask,
2131                                       Int_t       marker,
2132                                       Int_t       color,
2133                                       TList*      mclist,
2134                                       TList*      truthlist)
2135 {
2136   //
2137   // End of processing
2138   //
2139   // Parameters:
2140   //    sums        List of sums
2141   //    results     Output list of results
2142   //    shapeCorr   Shape correction or nil
2143   //    trigEff     Trigger efficiency
2144   //    symmetrice  Whether to symmetrice the results
2145   //    rebin       Whether to rebin the results
2146   //    corrEmpty   Whether to correct for empty bins
2147   //    cutEdges    Whether to cut edges when rebinning
2148   //    triggerMask Trigger mask
2149   //
2150   DGUARD(fDebug,1,"End centrality bin procesing");
2151
2152   fSums = dynamic_cast<TList*>(sums->FindObject(GetListName()));
2153   if(!fSums) {
2154     AliError("Could not retrieve TList fSums");
2155     return;
2156   }
2157
2158   fOutput = new TList;
2159   fOutput->SetName(GetListName());
2160   fOutput->SetOwner();
2161   results->Add(fOutput);
2162
2163   if (!fSum) {
2164     if (!ReadSum(fSums, false)) {
2165       AliInfo("This task did not produce any output");
2166       return;
2167     }
2168   }
2169   if (!fSumMC) ReadSum(fSums, true);
2170
2171   fTriggers = static_cast<TH1I*>(fSums->FindObject("triggers"));
2172
2173   if (!fTriggers) {
2174     AliError("Couldn't find histogram 'triggers' in list");
2175     return;
2176   }
2177
2178   // --- Get normalization scaler ------------------------------------
2179   Double_t epsilonT  = trigEff;
2180   Double_t epsilonT0 = trigEff0;
2181   AliInfoF("Using epsilonT=%f, epsilonT0=%f for 0x%x",
2182            epsilonT, epsilonT0, triggerMask);
2183 #if 0
2184   // These hard-coded trigger efficiencies are not used anymore, and
2185   // are only left in the code for reference.  We should remove this
2186   // soon.
2187   if (triggerMask == AliAODForwardMult::kNSD) {
2188     // This is a local change
2189     epsilonT = 0.96; // New value since HHD code 29/08/2012, why?
2190     //epsilonT = 0.92; //First paper...
2191     //epsilonT = 0.954; //First paper...
2192     //epsilonT = 1.03; //phojet
2193     AliWarning(Form("Using hard-coded NSD trigger efficiency of %f",epsilonT));
2194   }
2195   else if (triggerMask == AliAODForwardMult::kInel) {
2196     // This is a local change
2197     epsilonT = 0.934; // New value since HHD code 29/08/2012, why?
2198     // 900 GeV Inel eff from Martin
2199     //epsilonT = 0.916;
2200     //epsilonT = 0.97;  //phojet
2201
2202     AliWarning(Form("Using hard-coded Inel trigger efficiency of %f",epsilonT));
2203   }
2204   if (scheme & kZeroBin) {
2205     if (triggerMask==AliAODForwardMult::kInel)
2206       epsilonT0 = 0.785021; // 0.100240;
2207     else if (triggerMask==AliAODForwardMult::kInelGt0)
2208       epsilonT0 = 0;
2209     else if (triggerMask==AliAODForwardMult::kNSD)
2210       epsilonT0 = .706587;
2211     epsilonT = 1;
2212     AliWarning(Form("Using hard-coded NCluster>0 trigger efficiency of %f",
2213                     epsilonT0));
2214   }
2215 #endif
2216
2217   // Get our histograms
2218   Double_t nSum   = 0;
2219   TH2D*    sum    = fSum->CalcSum(fOutput, nSum, epsilonT0, epsilonT,
2220                                  marker, rootProj, corrEmpty);
2221   Double_t nSumMC = 0;
2222   TH2D*    sumMC  = 0;
2223   if (fSumMC) sumMC = fSumMC->CalcSum(fOutput, nSumMC,
2224                                       epsilonT0, epsilonT, marker,
2225                                       rootProj, corrEmpty);
2226   if (!sum) {
2227     AliError("Failed to get sum from summer - bailing out");
2228     return;
2229   }
2230
2231   TString  text;
2232   Double_t ntotal = nSum;
2233   Double_t scaler = Normalization(*fTriggers, scheme, epsilonT, ntotal, &text);
2234   if (scaler < 0) {
2235     AliError("Failed to calculate normalization - bailing out");
2236     return;
2237   }
2238   fOutput->Add(fTriggers->Clone());
2239   fOutput->Add(new TNamed("normCalc", text.Data()));
2240
2241   // --- Make result and store ---------------------------------------
2242   MakeResult(sum, "", rootProj, corrEmpty, (scheme & kShape) ? shapeCorr : 0,
2243              scaler, symmetrice, rebin, cutEdges, marker, color,
2244              mclist, truthlist);
2245
2246   // --- Process result from TrackRefs -------------------------------
2247   if (sumMC)
2248     MakeResult(sumMC, "MC", rootProj, corrEmpty,
2249                (scheme & kShape) ? shapeCorr : 0,
2250                scaler, symmetrice, rebin, cutEdges,
2251                GetMarkerStyle(GetMarkerBits(marker)+4), color,
2252                mclist, truthlist);
2253
2254   // Temporary stuff
2255   // if (!IsAllBin()) return;
2256
2257 }
2258 //____________________________________________________________________
2259 Bool_t
2260 AliBasedNdetaTask::ApplyEmpiricalCorrection(const AliAODForwardMult* aod,
2261                                             TH2D* data)
2262 {
2263   if (!fglobalempiricalcorrection || !data)
2264     return true;
2265   Float_t zvertex=aod->GetIpZ();
2266   Int_t binzvertex=fglobalempiricalcorrection->GetXaxis()->FindBin(zvertex);
2267   if(binzvertex<1||binzvertex>fglobalempiricalcorrection->GetNbinsX())
2268     return false;
2269   for (int i=1;i<=data->GetNbinsX();i++) {
2270     Int_t bincorrection=fglobalempiricalcorrection->GetYaxis()
2271       ->FindBin(data->GetXaxis()->GetBinCenter(i));
2272     if(bincorrection<1||bincorrection>fglobalempiricalcorrection->GetNbinsY())
2273       return false;
2274     Float_t correction=fglobalempiricalcorrection
2275       ->GetBinContent(binzvertex,bincorrection);
2276     if(correction<0.001) {
2277       data->SetBinContent(i,0,0);
2278       data->SetBinContent(i,data->GetNbinsY()+1,0);
2279     }
2280     for(int j=1;j<=data->GetNbinsY();j++) {
2281       if (data->GetBinContent(i,j)>0.0) {
2282         data->SetBinContent(i,j,data->GetBinContent(i,j)*correction);
2283         data->SetBinError(i,j,data->GetBinError(i,j)*correction);
2284       }
2285     }
2286   }
2287   return true;
2288 }
2289
2290 //
2291 // EOF
2292 //