PWG2/FORWARD/analysis2/AliFMDEnergyFitterTask.cxx

   1 //
   2 // Histogram and fit the energy loss distributions for the FMD
   3 //
   4 // Inputs:
   5 //   - AliESDEvent
   6 //
   7 // Outputs:
   8 //   - None
   9 //
  10 // Histograms:
  11 //
  12 // Corrections used:
  13 //   - None
  14 //
  15 //
  16 //
  17 #include "AliFMDEnergyFitterTask.h"
  18 #include "AliLog.h"
  19 #include "AliESDEvent.h"
  20 #include "AliAODForwardMult.h"
  21 #include "AliForwardCorrectionManager.h"
  22 #include "AliAnalysisManager.h"
  23 #include "AliAnalysisDataSlot.h"
  24 #include "AliAnalysisDataContainer.h"
  25 #include <TH1.h>
  26 #include <TDirectory.h>
  27 #include <TTree.h>
  28 #include <TFile.h>
  29 #include "AliMCEvent.h"
  30 #include "AliGenHijingEventHeader.h"
  31 #include "AliHeader.h"
  32 #include <iostream>
  33
  34 //====================================================================
  35 AliFMDEnergyFitterTask::AliFMDEnergyFitterTask()
  36   : AliAnalysisTaskSE(),
  37     fFirstEvent(true),
  38     fEventInspector(),
  39     fEnergyFitter(),
  40     fList(0),
  41     fbLow(0),
  42     fbHigh(100)
  43 {
  44   //
  45   // Constructor
  46   //
  47 }
  48
  49 //____________________________________________________________________
  50 AliFMDEnergyFitterTask::AliFMDEnergyFitterTask(const char* name)
  51   : AliAnalysisTaskSE(name),
  52     fFirstEvent(true),
  53     fEventInspector("event"),
  54     fEnergyFitter("energy"),
  55     fList(0),
  56     fbLow(0),
  57     fbHigh(100)
  58 {
  59   //
  60   // Constructor
  61   //
  62   // Parameters:
  63   //    name Name of task
  64   //
  65   DefineOutput(1, TList::Class());
  66   DefineOutput(2, TList::Class());
  67 }
  68
  69 //____________________________________________________________________
  70 AliFMDEnergyFitterTask::AliFMDEnergyFitterTask(const AliFMDEnergyFitterTask& o)
  71   : AliAnalysisTaskSE(o),
  72     fFirstEvent(o.fFirstEvent),
  73     fEventInspector(o.fEventInspector),
  74     fEnergyFitter(o.fEnergyFitter),
  75     fList(o.fList),
  76     fbLow(o.fbLow),
  77     fbHigh(o.fbHigh)
  78 {
  79   //
  80   // Copy constructor
  81   //
  82   // Parameters:
  83   //    o Object to copy from
  84   //
  85   DefineOutput(1, TList::Class());
  86   DefineOutput(2, TList::Class());
  87 }
  88
  89 //____________________________________________________________________
  90 AliFMDEnergyFitterTask&
  91 AliFMDEnergyFitterTask::operator=(const AliFMDEnergyFitterTask& o)
  92 {
  93   //
  94   // Assignment operator
  95   //
  96   // Parameters:
  97   //    o Object to assign from
  98   //
  99   // Return:
 100   //    Reference to this object
 101   //
 102   AliAnalysisTaskSE::operator=(o);
 103
 104   fFirstEvent        = o.fFirstEvent;
 105   fEventInspector    = o.fEventInspector;
 106   fEnergyFitter      = o.fEnergyFitter;
 107   fList              = o.fList;
 108   fbLow              = o.fbLow;
 109   fbHigh             = o.fbHigh;
 110
 111   return *this;
 112 }
 113
 114 //____________________________________________________________________
 115 void
 116 AliFMDEnergyFitterTask::SetDebug(Int_t dbg)
 117 {
 118   //
 119   // Set the debug level
 120   //
 121   // Parameters:
 122   //    dbg Debug level
 123   //
 124   fEventInspector.SetDebug(dbg);
 125   fEnergyFitter.SetDebug(dbg);
 126 }
 127
 128 //____________________________________________________________________
 129 void
 130 AliFMDEnergyFitterTask::Init()
 131 {
 132   //
 133   // Initialize the task
 134   //
 135   //
 136   fFirstEvent = true;
 137 }
 138
 139 //____________________________________________________________________
 140 void
 141 AliFMDEnergyFitterTask::InitializeSubs()
 142 {
 143   //
 144   // Initialise the sub objects and stuff.  Called on first event
 145   //
 146   //
 147   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
 148   fcm.Init(fEventInspector.GetCollisionSystem(),
 149            fEventInspector.GetEnergy(),
 150            fEventInspector.GetField(), 0);
 151   TAxis eAxis(0,0,0);
 152   TAxis vAxis(10,-10,10);
 153   fEnergyFitter.Init(eAxis);
 154   fEventInspector.Init(vAxis);
 155 }
 156
 157 //____________________________________________________________________
 158 void
 159 AliFMDEnergyFitterTask::UserCreateOutputObjects()
 160 {
 161   //
 162   // Create output objects
 163   //
 164   //
 165   fList = new TList;
 166
 167   fEventInspector.DefineOutput(fList);
 168   fEnergyFitter.DefineOutput(fList);
 169
 170   PostData(1, fList);
 171 }
 172 //____________________________________________________________________
 173 void
 174 AliFMDEnergyFitterTask::UserExec(Option_t*)
 175 {
 176   //
 177   // Process each event
 178   //
 179   // Parameters:
 180   //    option Not used
 181   //
 182
 183   // static Int_t cnt = 0;
 184   // cnt++;
 185   // Get the input data
 186
 187   AliMCEvent* mcevent = MCEvent();
 188   if(mcevent) {
 189     AliHeader* header            = mcevent->Header();
 190     AliGenHijingEventHeader* hijingHeader = dynamic_cast<AliGenHijingEventHeader*>(header->GenEventHeader());
 191     if(hijingHeader) {
 192       Float_t b = hijingHeader->ImpactParameter();
 193       if(b<fbLow || b>fbHigh) return;
 194       else
 195         std::cout<<"Selecting event with impact parameter "<<b<<std::endl;
 196     }
 197
 198   }
 199
 200   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 201   // AliInfo(Form("Event # %6d (esd=%p)", cnt, esd));
 202   if (!esd) {
 203     AliWarning("No ESD event found for input event");
 204     return;
 205   }
 206
 207   // On the first event, initialize the parameters
 208   if (fFirstEvent && esd->GetESDRun()) {
 209     fEventInspector.ReadRunDetails(esd);
 210
 211     AliInfo(Form("Initializing with parameters from the ESD:\n"
 212                  "         AliESDEvent::GetBeamEnergy()   ->%f\n"
 213                  "         AliESDEvent::GetBeamType()     ->%s\n"
 214                  "         AliESDEvent::GetCurrentL3()    ->%f\n"
 215                  "         AliESDEvent::GetMagneticField()->%f\n"
 216                  "         AliESDEvent::GetRunNumber()    ->%d\n",
 217                  esd->GetBeamEnergy(),
 218                  esd->GetBeamType(),
 219                  esd->GetCurrentL3(),
 220                  esd->GetMagneticField(),
 221                  esd->GetRunNumber()));
 222
 223
 224
 225     // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 226     // pars->SetParametersFromESD(esd);
 227     // pars->PrintStatus();
 228     fFirstEvent = false;
 229
 230     InitializeSubs();
 231   }
 232   Bool_t   lowFlux  = kFALSE;
 233   UInt_t   triggers = 0;
 234   UShort_t ivz      = 0;
 235   Double_t vz       = 0;
 236   Double_t cent     = 0;
 237   UInt_t   found    = fEventInspector.Process(esd, triggers, lowFlux,
 238                                               ivz, vz, cent);
 239   if (found & AliFMDEventInspector::kNoEvent)    return;
 240   if (found & AliFMDEventInspector::kNoTriggers) return;
 241   if (found & AliFMDEventInspector::kNoSPD)     return;
 242   if (found & AliFMDEventInspector::kNoFMD)     return;
 243   if (found & AliFMDEventInspector::kNoVertex)  return;
 244   if (found & AliFMDEventInspector::kBadVertex) return;
 245
 246   //  if(cent > 0) {
 247   //  if( cent < 40 || cent >50 ) return;
 248   //  else std::cout<<"selecting event with cent "<<cent<<std::endl;
 249   // }
 250
 251   // Get FMD data
 252   AliESDFMD* esdFMD = esd->GetFMDData();
 253   // Do the energy stuff
 254   if (!fEnergyFitter.Accumulate(*esdFMD, cent,
 255                                 triggers & AliAODForwardMult::kEmpty)){
 256     AliWarning("Energy fitter failed");
 257     return;
 258   }
 259   PostData(1, fList);
 260 }
 261
 262 //____________________________________________________________________
 263 void
 264 AliFMDEnergyFitterTask::Terminate(Option_t*)
 265 {
 266   //
 267   // End of job
 268   //
 269   // Parameters:
 270   //    option Not used
 271   //
 272   AliInfo(Form("Running terminate of %s", GetName()));
 273   TList* list = dynamic_cast<TList*>(GetOutputData(1));
 274   if (!list) {
 275     AliError(Form("No output list defined (%p)", GetOutputData(1)));
 276     if (GetOutputData(1)) GetOutputData(1)->Print();
 277     return;
 278   }
 279
 280 #if 0
 281   // Get our histograms from the container
 282   TH1I* hEventsTr    = 0;//static_cast<TH1I*>(list->FindObject("nEventsTr"));
 283   TH1I* hEventsTrVtx = 0;//static_cast<TH1I*>(list->FindObject("nEventsTrVtx"));
 284   TH1I* hTriggers    = 0;
 285   if (!fEventInspector.FetchHistograms(list, hEventsTr,
 286                                        hEventsTrVtx, hTriggers)) {
 287     AliError(Form("Didn't get histograms from event selector "
 288                   "(hEventsTr=%p,hEventsTrVtx=%p)",
 289                   hEventsTr, hEventsTrVtx));
 290     list->ls();
 291     return;
 292   }
 293 #endif
 294   AliInfo("Fitting energy loss spectra");
 295   fEnergyFitter.Fit(list);
 296
 297   // Make a deep copy and post that as output 2
 298   TList* list2 = static_cast<TList*>(list->Clone(Form("%sResults",
 299                                                       list->GetName())));
 300   PostData(2, list2);
 301 #if 0
 302   // Temporary code to save output to a file - should be disabled once
 303   // the plugin stuff works
 304   list->ls();
 305   TDirectory* savdir = gDirectory;
 306   TFile* tmp = TFile::Open("elossfits.root", "RECREATE");
 307   list->Write(list->GetName(), TObject::kSingleKey);
 308   tmp->Write();
 309   tmp->Close();
 310   savdir->cd();
 311 #endif
 312 }
 313
 314 //____________________________________________________________________
 315 void
 316 AliFMDEnergyFitterTask::Print(Option_t*) const
 317 {
 318   //
 319   // Print information
 320   //
 321   // Parameters:
 322   //    option Not used
 323   //
 324 }
 325
 326 //
 327 // EOF
 328 //