PWGLF/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   if (&o == this) return *this;
 103   AliAnalysisTaskSE::operator=(o);
 104
 105   fFirstEvent        = o.fFirstEvent;
 106   fEventInspector    = o.fEventInspector;
 107   fEnergyFitter      = o.fEnergyFitter;
 108   fList              = o.fList;
 109   fbLow              = o.fbLow;
 110   fbHigh             = o.fbHigh;
 111
 112   return *this;
 113 }
 114
 115 //____________________________________________________________________
 116 void
 117 AliFMDEnergyFitterTask::SetDebug(Int_t dbg)
 118 {
 119   //
 120   // Set the debug level
 121   //
 122   // Parameters:
 123   //    dbg Debug level
 124   //
 125   fEventInspector.SetDebug(dbg);
 126   fEnergyFitter.SetDebug(dbg);
 127 }
 128
 129 //____________________________________________________________________
 130 void
 131 AliFMDEnergyFitterTask::Init()
 132 {
 133   //
 134   // Initialize the task
 135   //
 136   //
 137   fFirstEvent = true;
 138 }
 139
 140 //____________________________________________________________________
 141 void
 142 AliFMDEnergyFitterTask::InitializeSubs()
 143 {
 144   //
 145   // Initialise the sub objects and stuff.  Called on first event
 146   //
 147   //
 148   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
 149   UShort_t sys = fEventInspector.GetCollisionSystem();
 150   UShort_t sNN = fEventInspector.GetEnergy();
 151   Short_t  fld = fEventInspector.GetField();
 152   fcm.Init(sys, sNN, fld, 0);
 153   TAxis eAxis(0,0,0);
 154   TAxis vAxis(10,-10,10);
 155   fEnergyFitter.Init(eAxis);
 156   fEventInspector.Init(vAxis);
 157
 158 }
 159
 160 //____________________________________________________________________
 161 void
 162 AliFMDEnergyFitterTask::UserCreateOutputObjects()
 163 {
 164   //
 165   // Create output objects
 166   //
 167   //
 168   fList = new TList;
 169   fList->SetOwner();
 170
 171   fEventInspector.DefineOutput(fList);
 172   fEnergyFitter.DefineOutput(fList);
 173
 174   PostData(1, fList);
 175 }
 176 //____________________________________________________________________
 177 void
 178 AliFMDEnergyFitterTask::UserExec(Option_t*)
 179 {
 180   //
 181   // Process each event
 182   //
 183   // Parameters:
 184   //    option Not used
 185   //
 186
 187   // static Int_t cnt = 0;
 188   // cnt++;
 189   // Get the input data
 190
 191   AliMCEvent* mcevent = MCEvent();
 192   if(mcevent) {
 193     AliHeader* header            = mcevent->Header();
 194     AliGenHijingEventHeader* hijingHeader =
 195       dynamic_cast<AliGenHijingEventHeader*>(header->GenEventHeader());
 196     if(hijingHeader) {
 197       Float_t b = hijingHeader->ImpactParameter();
 198       if(b<fbLow || b>fbHigh) return;
 199       else
 200         std::cout<<"Selecting event with impact parameter "<<b<<std::endl;
 201     }
 202
 203   }
 204
 205   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 206   // AliInfo(Form("Event # %6d (esd=%p)", cnt, esd));
 207   if (!esd) {
 208     AliWarning("No ESD event found for input event");
 209     return;
 210   }
 211
 212   // On the first event, initialize the parameters
 213   if (fFirstEvent && esd->GetESDRun()) {
 214     fEventInspector.ReadRunDetails(esd);
 215
 216     AliInfo(Form("Initializing with parameters from the ESD:\n"
 217                  "         AliESDEvent::GetBeamEnergy()   ->%f\n"
 218                  "         AliESDEvent::GetBeamType()     ->%s\n"
 219                  "         AliESDEvent::GetCurrentL3()    ->%f\n"
 220                  "         AliESDEvent::GetMagneticField()->%f\n"
 221                  "         AliESDEvent::GetRunNumber()    ->%d\n",
 222                  esd->GetBeamEnergy(),
 223                  esd->GetBeamType(),
 224                  esd->GetCurrentL3(),
 225                  esd->GetMagneticField(),
 226                  esd->GetRunNumber()));
 227
 228
 229
 230     // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 231     // pars->SetParametersFromESD(esd);
 232     // pars->PrintStatus();
 233     fFirstEvent = false;
 234
 235     InitializeSubs();
 236   }
 237   Bool_t   lowFlux   = kFALSE;
 238   UInt_t   triggers  = 0;
 239   UShort_t ivz       = 0;
 240   Double_t vz        = 0;
 241   Double_t cent      = 0;
 242   UShort_t nClusters = 0;
 243   UInt_t   found     = fEventInspector.Process(esd, triggers, lowFlux,
 244                                                ivz, vz, cent, nClusters);
 245   if (found & AliFMDEventInspector::kNoEvent)    return;
 246   if (found & AliFMDEventInspector::kNoTriggers) return;
 247   if (found & AliFMDEventInspector::kNoSPD)     return;
 248   if (found & AliFMDEventInspector::kNoFMD)     return;
 249   if (found & AliFMDEventInspector::kNoVertex)  return;
 250   if (found & AliFMDEventInspector::kBadVertex) return;
 251
 252   //  if(cent > 0) {
 253   //  if( cent < 40 || cent >50 ) return;
 254   //  else std::cout<<"selecting event with cent "<<cent<<std::endl;
 255   // }
 256
 257   // Get FMD data
 258   AliESDFMD* esdFMD = esd->GetFMDData();
 259   // Do the energy stuff
 260   if (!fEnergyFitter.Accumulate(*esdFMD, cent,
 261                                 triggers & AliAODForwardMult::kEmpty)){
 262     AliWarning("Energy fitter failed");
 263     return;
 264   }
 265   PostData(1, fList);
 266 }
 267
 268 //____________________________________________________________________
 269 void
 270 AliFMDEnergyFitterTask::Terminate(Option_t*)
 271 {
 272   //
 273   // End of job
 274   //
 275   // Parameters:
 276   //    option Not used
 277   //
 278   AliInfo(Form("Running terminate of %s", GetName()));
 279   TList* list = dynamic_cast<TList*>(GetOutputData(1));
 280   if (!list) {
 281     AliError(Form("No output list defined (%p)", GetOutputData(1)));
 282     if (GetOutputData(1)) GetOutputData(1)->Print();
 283     return;
 284   }
 285
 286   AliInfo("Fitting energy loss spectra");
 287   fEnergyFitter.Fit(list);
 288
 289   // Make a deep copy and post that as output 2
 290   TList* list2 = static_cast<TList*>(list->Clone(Form("%sResults",
 291                                                       list->GetName())));
 292   list2->SetOwner();
 293   PostData(2, list2);
 294 }
 295
 296 //____________________________________________________________________
 297 void
 298 AliFMDEnergyFitterTask::Print(Option_t*) const
 299 {
 300   //
 301   // Print information
 302   //
 303   // Parameters:
 304   //    option Not used
 305   //
 306 }
 307
 308 //
 309 // EOF
 310 //