PWG2/FORWARD/analysis2/AliForwardQATask.cxx

   1 //
   2 // Calculate the multiplicity in the forward regions event-by-event
   3 //
   4 // Inputs:
   5 //   - AliESDEvent
   6 //
   7 // Outputs:
   8 //   - AliAODForwardMult
   9 //
  10 // Histograms
  11 //
  12 // Corrections used
  13 //
  14 #include "AliForwardQATask.h"
  15 #include "AliForwardUtil.h"
  16 #include "AliTriggerAnalysis.h"
  17 #include "AliPhysicsSelection.h"
  18 #include "AliLog.h"
  19 #include "AliESDEvent.h"
  20 #include "AliAODHandler.h"
  21 #include "AliMultiplicity.h"
  22 #include "AliInputEventHandler.h"
  23 #include "AliForwardCorrectionManager.h"
  24 #include "AliAnalysisManager.h"
  25 #include "AliAODForwardMult.h"
  26 #include <TH1.h>
  27 #include <TDirectory.h>
  28 #include <TTree.h>
  29 #include <TROOT.h>
  30 #include <TStopwatch.h>
  31
  32 //====================================================================
  33 AliForwardQATask::AliForwardQATask()
  34   : AliAnalysisTaskSE(),
  35     fEnableLowFlux(false),
  36     fFirstEvent(true),
  37     fCorrManager(0),
  38     fESDFMD(),
  39     fHistos(),
  40     fEventInspector(),
  41     fEnergyFitter(),
  42     fSharingFilter(),
  43     fDensityCalculator(),
  44     fList(0),
  45     fDebug(0)
  46 {
  47   //
  48   // Constructor
  49   //
  50 }
  51
  52 //____________________________________________________________________
  53 AliForwardQATask::AliForwardQATask(const char* name)
  54   : AliAnalysisTaskSE(name),
  55     fEnableLowFlux(false),
  56     fFirstEvent(true),
  57     fCorrManager(0),
  58     fESDFMD(),
  59     fHistos(),
  60     fEventInspector("event"),
  61     fEnergyFitter("energy"),
  62     fSharingFilter("sharing"),
  63     fDensityCalculator("density"),
  64     fList(0),
  65     fDebug(0)
  66 {
  67   //
  68   // Constructor
  69   //
  70   // Parameters:
  71   //    name Name of task
  72   //
  73   DefineOutput(1, TList::Class());
  74   DefineOutput(2, TList::Class());
  75   fCorrManager = &AliForwardCorrectionManager::Instance();
  76   fEnergyFitter.SetNParticles(1); // Just find the 1st peak
  77   fEnergyFitter.SetDoMakeObject(false);
  78   fEnergyFitter.SetUseIncreasingBins(true);
  79   fEnergyFitter.SetDoFits(kTRUE);
  80   fEnergyFitter.SetLowCut(0.4);
  81   fEnergyFitter.SetFitRangeBinWidth(4);
  82   fEnergyFitter.SetMinEntries(1000);
  83 }
  84
  85 //____________________________________________________________________
  86 AliForwardQATask::AliForwardQATask(const AliForwardQATask& o)
  87   : AliAnalysisTaskSE(o),
  88     fEnableLowFlux(o.fEnableLowFlux),
  89     fFirstEvent(o.fFirstEvent),
  90     fCorrManager(o.fCorrManager),
  91     fESDFMD(o.fESDFMD),
  92     fHistos(o.fHistos),
  93     fEventInspector(o.fEventInspector),
  94     fEnergyFitter(o.fEnergyFitter),
  95     fSharingFilter(o.fSharingFilter),
  96     fDensityCalculator(o.fDensityCalculator),
  97     fList(o.fList),
  98     fDebug(o.fDebug)
  99 {
 100   //
 101   // Copy constructor
 102   //
 103   // Parameters:
 104   //    o Object to copy from
 105   //
 106   DefineOutput(1, TList::Class());
 107   DefineOutput(2, TList::Class());
 108 }
 109
 110 //____________________________________________________________________
 111 AliForwardQATask&
 112 AliForwardQATask::operator=(const AliForwardQATask& o)
 113 {
 114   //
 115   // Assignment operator
 116   //
 117   // Parameters:
 118   //    o Object to assign from
 119   //
 120   // Return:
 121   //    Reference to this object
 122   //
 123   AliAnalysisTaskSE::operator=(o);
 124
 125   fEnableLowFlux     = o.fEnableLowFlux;
 126   fFirstEvent        = o.fFirstEvent;
 127   fCorrManager       = o.fCorrManager;
 128   fEventInspector    = o.fEventInspector;
 129   fEnergyFitter      = o.fEnergyFitter;
 130   fSharingFilter     = o.fSharingFilter;
 131   fDensityCalculator = o.fDensityCalculator;
 132   fHistos            = o.fHistos;
 133   fList              = o.fList;
 134   fDebug             = o.fDebug;
 135
 136   return *this;
 137 }
 138
 139 //____________________________________________________________________
 140 void
 141 AliForwardQATask::SetDebug(Int_t dbg)
 142 {
 143   //
 144   // Set debug level
 145   //
 146   // Parameters:
 147   //    dbg Debug level
 148   //
 149   fDebug = dbg;
 150   fEventInspector.SetDebug(dbg);
 151   fEnergyFitter.SetDebug(dbg);
 152   fSharingFilter.SetDebug(dbg);
 153   fDensityCalculator.SetDebug(dbg);
 154 }
 155
 156 //____________________________________________________________________
 157 Bool_t
 158 AliForwardQATask::CheckCorrections(UInt_t what) const
 159 {
 160   //
 161   // Check if all needed corrections are there and accounted for.  If not,
 162   // do a Fatal exit
 163   //
 164   // Parameters:
 165   //    what Which corrections is needed
 166   //
 167   // Return:
 168   //    true if all present, false otherwise
 169   //
 170
 171   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
 172   // Check that we have the energy loss fits, needed by
 173   //   AliFMDSharingFilter
 174   //   AliFMDDensityCalculator
 175   if (what & AliForwardCorrectionManager::kELossFits && !fcm.GetELossFit()) {
 176     AliWarning("No energy loss fits");
 177     return false;
 178   }
 179   return true;
 180 }
 181
 182 //____________________________________________________________________
 183 Bool_t
 184 AliForwardQATask::ReadCorrections(const TAxis*& pe,
 185                                   const TAxis*& pv,
 186                                   Bool_t        mc)
 187 {
 188   //
 189   // Read corrections
 190   //
 191   //
 192   UInt_t what = AliForwardCorrectionManager::kAll;
 193   what ^= AliForwardCorrectionManager::kDoubleHit;
 194   what ^= AliForwardCorrectionManager::kVertexBias;
 195   what ^= AliForwardCorrectionManager::kAcceptance;
 196   what ^= AliForwardCorrectionManager::kMergingEfficiency;
 197
 198   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
 199   if (!fcm.Init(GetEventInspector().GetCollisionSystem(),
 200                 GetEventInspector().GetEnergy(),
 201                 GetEventInspector().GetField(),
 202                 mc,
 203                 what)) return false;
 204   if (!CheckCorrections(what)) {
 205     return false;
 206   }
 207
 208   // Sett our persistency pointer
 209   // fCorrManager = &fcm;
 210
 211   // Get the eta axis from the secondary maps - if read in
 212   if (!pe) {
 213     pe = fcm.GetEtaAxis();
 214     if (!pe) AliFatal("No eta axis defined");
 215   }
 216   // Get the vertex axis from the secondary maps - if read in
 217   if (!pv) {
 218     pv = fcm.GetVertexAxis();
 219     if (!pv) AliFatal("No vertex axis defined");
 220   }
 221
 222   return true;
 223 }
 224
 225 //____________________________________________________________________
 226 AliESDEvent*
 227 AliForwardQATask::GetESDEvent()
 228 {
 229   //
 230   // Get the ESD event. IF this is the first event, initialise
 231   //
 232   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 233   if (!esd) {
 234     AliWarning("No ESD event found for input event");
 235     return 0;
 236   }
 237
 238   // On the first event, initialize the parameters
 239   if (fFirstEvent && esd->GetESDRun()) {
 240     GetEventInspector().ReadRunDetails(esd);
 241
 242     AliInfoF("Initializing with parameters from the ESD:\n"
 243              "         AliESDEvent::GetBeamEnergy()   ->%f\n"
 244              "         AliESDEvent::GetBeamType()     ->%s\n"
 245              "         AliESDEvent::GetCurrentL3()    ->%f\n"
 246              "         AliESDEvent::GetMagneticField()->%f\n"
 247              "         AliESDEvent::GetRunNumber()    ->%d\n",
 248              esd->GetBeamEnergy(),
 249              esd->GetBeamType(),
 250              esd->GetCurrentL3(),
 251              esd->GetMagneticField(),
 252              esd->GetRunNumber());
 253
 254
 255     if (!InitializeSubs()) {
 256       AliWarning("Initialisation of sub algorithms failed!");
 257       return 0;
 258     }
 259     AliInfoF("Clearing first event flag from %s to false",
 260              fFirstEvent ? "true" : "false");
 261     fFirstEvent = false;
 262   }
 263   return esd;
 264 }
 265 //____________________________________________________________________
 266 Bool_t
 267 AliForwardQATask::InitializeSubs()
 268 {
 269   //
 270   // Initialise the sub objects and stuff.  Called on first event
 271   //
 272   //
 273   const TAxis* pe = 0;
 274   const TAxis* pv = 0;
 275
 276   if (!ReadCorrections(pe,pv)) {
 277     AliWarning("Failed to read corrections");
 278     return false;
 279   }
 280
 281   fHistos.Init(*pe);
 282
 283   fEventInspector.Init(*pv);
 284   fEnergyFitter.Init(*pe);
 285   fSharingFilter.Init();
 286   fDensityCalculator.Init(*pe);
 287
 288   this->Print();
 289
 290   return true;
 291 }
 292
 293 //____________________________________________________________________
 294 void
 295 AliForwardQATask::UserCreateOutputObjects()
 296 {
 297   //
 298   // Create output objects
 299   //
 300   //
 301   fList = new TList;
 302   fList->SetOwner();
 303
 304   fEventInspector.DefineOutput(fList);
 305   fEnergyFitter.DefineOutput(fList);
 306   fSharingFilter.DefineOutput(fList);
 307   fDensityCalculator.DefineOutput(fList);
 308
 309   PostData(1, fList);
 310 }
 311 //____________________________________________________________________
 312 void
 313 AliForwardQATask::UserExec(Option_t*)
 314 {
 315   //
 316   // Process each event
 317   //
 318   // Parameters:
 319   //    option Not used
 320   //
 321
 322   // static Int_t cnt = 0;
 323   // cnt++;
 324   // Get the input data
 325   AliESDEvent* esd = GetESDEvent();
 326   if (!esd) {
 327     AliWarning("Got no ESD event");
 328     return;
 329   }
 330   if (fFirstEvent) {
 331     // If the first event flag wasn't cleared in the above call to
 332     // GetESDEvent, we should not do anything, since nothing has been
 333     // initialised yet, so we opt out here (with a warning)
 334     AliWarning("Nothing has been initialized yet, opt'ing out");
 335     return;
 336   }
 337
 338   // Clear stuff
 339   fHistos.Clear();
 340   fESDFMD.Clear();
 341
 342   Bool_t   lowFlux   = kFALSE;
 343   UInt_t   triggers  = 0;
 344   UShort_t ivz       = 0;
 345   Double_t vz        = 0;
 346   Double_t cent      = -1;
 347   UShort_t nClusters = 0;
 348   UInt_t   found     = fEventInspector.Process(esd, triggers, lowFlux,
 349                                                ivz, vz, cent, nClusters);
 350
 351   if (found & AliFMDEventInspector::kNoEvent)    return;
 352   if (found & AliFMDEventInspector::kNoTriggers) return;
 353   if (found & AliFMDEventInspector::kNoSPD)      return;
 354   if (found & AliFMDEventInspector::kNoFMD)      return;
 355   if (found & AliFMDEventInspector::kNoVertex)   return;
 356   if (triggers & AliAODForwardMult::kPileUp)     return;
 357   if (found & AliFMDEventInspector::kBadVertex)  return;
 358
 359   // We we do not want to use low flux specific code, we disable it here.
 360   if (!fEnableLowFlux) lowFlux = false;
 361
 362   // Get FMD data
 363   AliESDFMD* esdFMD = esd->GetFMDData();
 364
 365   // Run the energy loss fitter
 366   if (!fEnergyFitter.Accumulate(*esdFMD, cent,
 367                                 triggers & AliAODForwardMult::kEmpty)) {
 368     AliWarning("Energy fitter failed");
 369     return;
 370   }
 371
 372   //  // Apply the sharing filter (or hit merging or clustering if you like)
 373   if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD)) {
 374     AliWarning("Sharing filter failed!");
 375     return;
 376   }
 377
 378   // Calculate the inclusive charged particle density
 379   if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) {
 380     // if (!fDensityCalculator.Calculate(*esdFMD, fHistos, ivz, lowFlux)) {
 381     AliWarning("Density calculator failed!");
 382     return;
 383   }
 384   PostData(1, fList);
 385 }
 386
 387 //____________________________________________________________________
 388 void
 389 AliForwardQATask::Terminate(Option_t*)
 390 {
 391   //
 392   // End of job
 393   //
 394   // Parameters:
 395   //    option Not used
 396   //
 397   if (fDebug) AliInfo("In Forwards terminate");
 398   TStopwatch swt;
 399   swt.Start();
 400
 401   TList* list = dynamic_cast<TList*>(GetOutputData(1));
 402   if (!list) {
 403     AliError(Form("No output list defined (%p)", GetOutputData(1)));
 404     if (GetOutputData(1)) GetOutputData(1)->Print();
 405     return;
 406   }
 407
 408   // Get our histograms from the container
 409   TH1I* hEventsTr    = 0;//static_cast<TH1I*>(list->FindObject("nEventsTr"));
 410   TH1I* hEventsTrVtx = 0;//static_cast<TH1I*>(list->FindObject("nEventsTrVtx"));
 411   TH1I* hTriggers    = 0;
 412   if (!fEventInspector.FetchHistograms(list, hEventsTr,
 413                                        hEventsTrVtx, hTriggers)) {
 414     AliError(Form("Didn't get histograms from event selector "
 415                   "(hEventsTr=%p,hEventsTrVtx=%p)",
 416                   hEventsTr, hEventsTrVtx));
 417     return;
 418   }
 419
 420   TStopwatch swf;
 421   swf.Start();
 422   fEnergyFitter.Fit(list);
 423   swf.Stop();
 424   AliInfoF("Fitting took %d real-time seconds, and %f CPU seconds",
 425            Int_t(swf.RealTime()), swf.CpuTime());
 426
 427   fSharingFilter.ScaleHistograms(list,Int_t(hEventsTr->Integral()));
 428   fDensityCalculator.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
 429
 430   // Make a deep copy and post that as output 2
 431   TList* list2 = static_cast<TList*>(list->Clone(Form("%sResults",
 432                                                       list->GetName())));
 433   if (fDebug) AliInfoF("Posting post processing results to %s",
 434                        list2->GetName());
 435   list2->SetOwner();
 436   PostData(2, list2);
 437
 438   swt.Stop();
 439   AliInfoF("Terminate took %d real-time seconds, and %f CPU seconds",
 440            Int_t(swt.RealTime()), swt.CpuTime());
 441
 442 }
 443
 444 //____________________________________________________________________
 445 void
 446 AliForwardQATask::Print(Option_t* option) const
 447 {
 448   //
 449   // Print information
 450   //
 451   // Parameters:
 452   //    option Not used
 453   //
 454
 455   std::cout << ClassName() << ": " << GetName() << "\n"
 456             << "  Enable low flux code:   " << (fEnableLowFlux ? "yes" : "no")
 457             << "\n"
 458             << "  Off-line trigger mask:  0x"
 459             << std::hex     << std::setfill('0')
 460             << std::setw (8) << fOfflineTriggerMask
 461             << std::dec     << std::setfill (' ') << std::endl;
 462   gROOT->IncreaseDirLevel();
 463   if (fCorrManager) fCorrManager->Print();
 464   else
 465     std::cout << "  Correction manager not set yet" << std::endl;
 466   GetEventInspector()   .Print(option);
 467   GetEnergyFitter()     .Print(option);
 468   GetSharingFilter()    .Print(option);
 469   gROOT->DecreaseDirLevel();
 470 }
 471
 472 //
 473 // EOF
 474 //