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 "AliAnalysisManager.h"
  22 #include "AliAnalysisDataSlot.h"
  23 #include "AliAnalysisDataContainer.h"
  24 #include <TH1.h>
  25 #include <TDirectory.h>
  26 #include <TTree.h>
  27 #include <TFile.h>
  28 #include "AliMCEvent.h"
  29 #include "AliGenHijingEventHeader.h"
  30 #include "AliHeader.h"
  31 #include <iostream>
  32
  33 //====================================================================
  34 AliFMDEnergyFitterTask::AliFMDEnergyFitterTask()
  35   : AliAnalysisTaskSE(),
  36     fFirstEvent(true),
  37     fEventInspector(),
  38     fEnergyFitter(),
  39     fList(0),
  40     fbLow(0),
  41     fbHigh(100)
  42 {
  43   //
  44   // Constructor
  45   //
  46   DGUARD(fDebug, 3,"Default CTOR of AliFMDEnergyFitterTask");
  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   DGUARD(fDebug, 3,"Named CTOR of AliFMDEnergyFitterTask: %s", name);
  66   DefineOutput(1, TList::Class());
  67   DefineOutput(2, TList::Class());
  68   fBranchNames =
  69     "ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.,"
  70     "AliESDFMD.,SPDVertex.,PrimaryVertex.";
  71 }
  72
  73 //____________________________________________________________________
  74 AliFMDEnergyFitterTask::AliFMDEnergyFitterTask(const AliFMDEnergyFitterTask& o)
  75   : AliAnalysisTaskSE(o),
  76     fFirstEvent(o.fFirstEvent),
  77     fEventInspector(o.fEventInspector),
  78     fEnergyFitter(o.fEnergyFitter),
  79     fList(o.fList),
  80     fbLow(o.fbLow),
  81     fbHigh(o.fbHigh)
  82 {
  83   //
  84   // Copy constructor
  85   //
  86   // Parameters:
  87   //    o Object to copy from
  88   //
  89   DGUARD(fDebug, 3,"COPY CTOR of AliFMDEnergyFitterTask");
  90   DefineOutput(1, TList::Class());
  91   DefineOutput(2, TList::Class());
  92 }
  93
  94 //____________________________________________________________________
  95 AliFMDEnergyFitterTask&
  96 AliFMDEnergyFitterTask::operator=(const AliFMDEnergyFitterTask& o)
  97 {
  98   //
  99   // Assignment operator
 100   //
 101   // Parameters:
 102   //    o Object to assign from
 103   //
 104   // Return:
 105   //    Reference to this object
 106   //
 107   fDebug = o.fDebug;
 108   DGUARD(fDebug,3,"Assignment of AliFMDEnergyFitterTask");
 109   if (&o == this) return *this;
 110   AliAnalysisTaskSE::operator=(o);
 111
 112   fFirstEvent        = o.fFirstEvent;
 113   fEventInspector    = o.fEventInspector;
 114   fEnergyFitter      = o.fEnergyFitter;
 115   fList              = o.fList;
 116   fbLow              = o.fbLow;
 117   fbHigh             = o.fbHigh;
 118
 119   return *this;
 120 }
 121
 122 //____________________________________________________________________
 123 void
 124 AliFMDEnergyFitterTask::SetDebug(Int_t dbg)
 125 {
 126   //
 127   // Set the debug level
 128   //
 129   // Parameters:
 130   //    dbg Debug level
 131   //
 132   fDebug = dbg;
 133   fEventInspector.SetDebug(dbg);
 134   fEnergyFitter.SetDebug(dbg);
 135 }
 136
 137 //____________________________________________________________________
 138 void
 139 AliFMDEnergyFitterTask::Init()
 140 {
 141   //
 142   // Initialize the task
 143   //
 144   //
 145   DGUARD(fDebug,1,"Initialize of AliFMDEnergyFitterTask");
 146   fFirstEvent = true;
 147 }
 148
 149 //____________________________________________________________________
 150 void
 151 AliFMDEnergyFitterTask::SetupForData()
 152 {
 153   //
 154   // Initialise the sub objects and stuff.  Called on first event
 155   //
 156   //
 157   DGUARD(fDebug,1,"Initialize subs of AliFMDEnergyFitterTask");
 158   TAxis eAxis(0,0,0); // Default only
 159   TAxis vAxis(10,-10,10); // Default only
 160   fEnergyFitter.SetupForData(eAxis);
 161   fEventInspector.SetupForData(vAxis);
 162
 163 }
 164
 165 //____________________________________________________________________
 166 void
 167 AliFMDEnergyFitterTask::UserCreateOutputObjects()
 168 {
 169   //
 170   // Create output objects
 171   //
 172   //
 173   DGUARD(fDebug,1,"Create output objects of AliFMDEnergyFitterTask");
 174   fList = new TList;
 175   fList->SetOwner();
 176
 177   fEventInspector.CreateOutputObjects(fList);
 178   fEnergyFitter.CreateOutputObjects(fList);
 179
 180   PostData(1, fList);
 181 }
 182 //____________________________________________________________________
 183 void
 184 AliFMDEnergyFitterTask::UserExec(Option_t*)
 185 {
 186   //
 187   // Process each event
 188   //
 189   // Parameters:
 190   //    option Not used
 191   //
 192
 193   // static Int_t cnt = 0;
 194   // cnt++;
 195   // Get the input data
 196   DGUARD(fDebug,3,"Analyse event of AliFMDEnergyFitterTask");
 197
 198   AliMCEvent* mcevent = MCEvent();
 199   if(mcevent) {
 200     AliHeader* header            = mcevent->Header();
 201     AliGenHijingEventHeader* hijingHeader =
 202       dynamic_cast<AliGenHijingEventHeader*>(header->GenEventHeader());
 203     if(hijingHeader) {
 204       Float_t b = hijingHeader->ImpactParameter();
 205       if(b<fbLow || b>fbHigh) return;
 206       else
 207         std::cout<<"Selecting event with impact parameter "<<b<<std::endl;
 208     }
 209
 210   }
 211
 212   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 213   // AliInfo(Form("Event # %6d (esd=%p)", cnt, esd));
 214   if (!esd) {
 215     AliWarning("No ESD event found for input event");
 216     return;
 217   }
 218
 219   // --- Read in the data --------------------------------------------
 220   LoadBranches();
 221
 222   // On the first event, initialize the parameters
 223   if (fFirstEvent && esd->GetESDRun()) {
 224     fEventInspector.ReadRunDetails(esd);
 225
 226     AliInfo(Form("Initializing with parameters from the ESD:\n"
 227                  "         AliESDEvent::GetBeamEnergy()   ->%f\n"
 228                  "         AliESDEvent::GetBeamType()     ->%s\n"
 229                  "         AliESDEvent::GetCurrentL3()    ->%f\n"
 230                  "         AliESDEvent::GetMagneticField()->%f\n"
 231                  "         AliESDEvent::GetRunNumber()    ->%d\n",
 232                  esd->GetBeamEnergy(),
 233                  esd->GetBeamType(),
 234                  esd->GetCurrentL3(),
 235                  esd->GetMagneticField(),
 236                  esd->GetRunNumber()));
 237
 238
 239
 240     // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 241     // pars->SetParametersFromESD(esd);
 242     // pars->PrintStatus();
 243     fFirstEvent = false;
 244
 245     SetupForData();
 246   }
 247   Bool_t   lowFlux   = kFALSE;
 248   UInt_t   triggers  = 0;
 249   UShort_t ivz       = 0;
 250   TVector3 ip;
 251   Double_t cent      = 0;
 252   UShort_t nClusters = 0;
 253   UInt_t   found     = fEventInspector.Process(esd, triggers, lowFlux,
 254                                                ivz, ip, cent, nClusters);
 255   if (found & AliFMDEventInspector::kNoEvent)    return;
 256   if (found & AliFMDEventInspector::kNoTriggers) return;
 257   if (found & AliFMDEventInspector::kNoSPD)     return;
 258   if (found & AliFMDEventInspector::kNoFMD)     return;
 259   if (found & AliFMDEventInspector::kNoVertex)  return;
 260   if (found & AliFMDEventInspector::kBadVertex) return;
 261
 262   //  if(cent > 0) {
 263   //  if( cent < 40 || cent >50 ) return;
 264   //  else std::cout<<"selecting event with cent "<<cent<<std::endl;
 265   // }
 266
 267   // Get FMD data
 268   AliESDFMD* esdFMD = esd->GetFMDData();
 269   // Do the energy stuff
 270   if (!fEnergyFitter.Accumulate(*esdFMD, cent,
 271                                 triggers & AliAODForwardMult::kEmpty)){
 272     AliWarning("Energy fitter failed");
 273     return;
 274   }
 275   PostData(1, fList);
 276 }
 277
 278 //____________________________________________________________________
 279 void
 280 AliFMDEnergyFitterTask::FinishTaskOutput()
 281 {
 282   if (!fList) {
 283     Warning("FinishTaskOutput", "No list defined");
 284   }
 285   // else {
 286   //   fList->ls();
 287   // }
 288   // if (fDebug)
 289   // gDebug = 1;
 290   AliAnalysisTaskSE::FinishTaskOutput();
 291 }
 292
 293 //____________________________________________________________________
 294 void
 295 AliFMDEnergyFitterTask::Terminate(Option_t*)
 296 {
 297   //
 298   // End of job
 299   //
 300   // Parameters:
 301   //    option Not used
 302   //
 303   DGUARD(fDebug,1,"Processing merged output of AliFMDEnergyFitterTask");
 304
 305   TList* list = dynamic_cast<TList*>(GetOutputData(1));
 306   if (!list) {
 307     AliError(Form("No output list defined (%p)", GetOutputData(1)));
 308     if (GetOutputData(1)) GetOutputData(1)->Print();
 309     return;
 310   }
 311
 312   AliInfo("Fitting energy loss spectra");
 313   fEnergyFitter.Fit(list);
 314
 315   // Make a deep copy and post that as output 2
 316   TList* list2 = static_cast<TList*>(list->Clone(Form("%sResults",
 317                                                       list->GetName())));
 318   list2->SetOwner();
 319   PostData(2, list2);
 320 }
 321
 322 //____________________________________________________________________
 323 void
 324 AliFMDEnergyFitterTask::Print(Option_t*) const
 325 {
 326   //
 327   // Print information
 328   //
 329   // Parameters:
 330   //    option Not used
 331   //
 332 }
 333
 334 //
 335 // EOF
 336 //