PWGLF/FORWARD/analysis2/trains/ELossTimeTask.C

   1 #ifndef __CINT__
   2 # include <AliFMDEventInspector.h>
   3 # include <AliForwardCorrectionManager.h>
   4 # include <TH2.h>
   5 # include <AliESDFMD.h>
   6 # include <AliAODForwardMult.h>
   7 # include <AliESDEvent.h>
   8 # include <TFile.h>
   9 # include <TSystem.h>
  10 #else
  11 class AliFMDEventInspector;
  12 class TH2;
  13 class AliESDFMD;
  14 #endif
  15 #include <AliBaseESDTask.h>
  16 #include <AliForwardUtil.h>
  17 #define NO_TASK
  18 #define NO_SORTER
  19 #include "EventTimeTask.C"
  20
  21 /**
  22  * @defgroup pwglf_forward_eventtime Investigation of time-between-events
  23  * @ingroup pwglf_forward_trains_specific
  24  */
  25 /**
  26  * Task to analyse the energy loss in the FMD rings as function of the
  27  * time to previous event.
  28  *
  29  * @ingroup pwglf_forward_eventtime
  30  */
  31 struct ELossTimeTask : public AliBaseESDTask
  32 {
  33   /**
  34    * Constructor - for I/O only
  35    */
  36   ELossTimeTask()
  37     : AliBaseESDTask(),
  38       fEventInspector(),
  39       fFMD1i(),
  40       fFMD2i(),
  41       fFMD2o(),
  42       fFMD3i(),
  43       fFMD3o(),
  44       fMap(0),
  45       fDt(0)
  46   {}
  47   /**
  48    * Constructor
  49    *
  50    * @param name Name of task
  51    */
  52   ELossTimeTask(const char* name)
  53     : AliBaseESDTask(name, "ELossTimeTask",
  54                      &(AliForwardCorrectionManager::Instance())),
  55       fEventInspector("fmdEventInspector"),
  56       fFMD1i(1, 'i'),
  57       fFMD2i(2, 'i'),
  58       fFMD2o(2, 'o'),
  59       fFMD3i(3, 'i'),
  60       fFMD3o(3, 'o'),
  61       fMap(0),
  62       fDt(0)
  63   {
  64   }
  65   /**
  66    * Book output objects. Derived class should define this to book
  67    * output objects on the processing output list @c fList before the
  68    * actual event processing.  This is called on the master and on
  69    * each slave.
  70    *
  71    * If this member function returns false, the execution is stopped
  72    * with a fatal signal.
  73    *
  74    * @return true on success.
  75    */
  76   virtual Bool_t Book()
  77   {
  78     fNeededCorrections = 0;
  79     fExtraCorrections  = 0;
  80
  81     fDt = new TH1D("dt", "Time-to-last event (PS-triggered)", 60, 0, 15);
  82     fDt->SetXTitle("log_{10}(#Deltat)");
  83     fDt->SetFillColor(kYellow+2);
  84     fDt->SetLineColor(kYellow+2);
  85     fDt->SetFillStyle(3001);
  86     fDt->SetDirectory(0);
  87     fList->Add(fDt);
  88
  89     fFMD1i.Book(fList, fDt);
  90     fFMD2i.Book(fList, fDt);
  91     fFMD2o.Book(fList, fDt);
  92     fFMD3i.Book(fList, fDt);
  93     fFMD3o.Book(fList, fDt);
  94
  95     // Possibly re-read map
  96     ReadMap("map.root");
  97
  98     return true;
  99   }
 100   /**
 101    * Process a single event
 102    *
 103    * @param esd Input event
 104    *
 105    * @return true on success
 106    */
 107   virtual Bool_t Event(AliESDEvent& esd)
 108   {
 109     Bool_t   lowFlux   = kFALSE;
 110     UInt_t   triggers  = 0;
 111     UShort_t ivz       = 0;
 112     TVector3 ip;
 113     Double_t cent      = 0;
 114     UShort_t nClusters = 0;
 115     UInt_t   found     = fEventInspector.Process(&esd, triggers, lowFlux,
 116                                                  ivz, ip, cent, nClusters);
 117     if (found & AliFMDEventInspector::kNoEvent)    return false;
 118     if (found & AliFMDEventInspector::kNoTriggers) return false;
 119     if (found & AliFMDEventInspector::kNoSPD)      return false;
 120     if (found & AliFMDEventInspector::kNoFMD)      return false;
 121     if (found & AliFMDEventInspector::kNoVertex)   return false;
 122     if (found & AliFMDEventInspector::kBadVertex)  return false;
 123
 124     // do not process pile-up, A, C, and E events
 125     if (triggers & AliAODForwardMult::kPileUp)     return false;
 126     if (triggers & AliAODForwardMult::kA)          return false;
 127     if (triggers & AliAODForwardMult::kC)          return false;
 128     if (triggers & AliAODForwardMult::kE)          return false;
 129
 130     // We want only the events found by off-line
 131     if (!(triggers & AliAODForwardMult::kOffline)) return false;
 132
 133     // Perhaps we should also insist on MB only
 134     // if (fOnlyMB && (!(triggers & AliAODForwardMult::kInel))) return false;
 135
 136     // Get FMD data
 137     AliESDFMD* esdFMD = esd.GetFMDData();
 138     ULong64_t  period = esd.GetPeriodNumber();
 139     ULong64_t  orbit  = esd.GetOrbitNumber();
 140     ULong64_t  bc     = esd.GetBunchCrossNumber();
 141     ULong64_t  full   = EventTimeData::EncodeFull(period, orbit, bc);
 142     ULong64_t  dt     = fMap->Get(full);
 143     Double_t   logDt  = TMath::Log10(25. * dt);
 144     if (dt == EventTimeMap::kInvalidTime) {
 145       logDt = 0;
 146       Printf("!!! Failed to find dT for 0x%016llu", full);
 147     }
 148     // else
 149     //   Printf("=== 0x%016llu -> 0x%016llu", full, dt);
 150     fDt->Fill(logDt);
 151
 152     for (UShort_t d = 1; d <= 3; d++) {
 153       UShort_t nQ = d == 1 ? 1 : 2;
 154       for (UShort_t q = 0; q < nQ; q++) {
 155         RingHistos* r = 0;
 156         switch (d) {
 157         case 1: r = &fFMD1i; break;
 158         case 2: r = (q == 0 ? &fFMD2i : &fFMD2o); break;
 159         case 3: r = (q == 0 ? &fFMD3i : &fFMD3o); break;
 160         }
 161         r->Event(*esdFMD, logDt);
 162       }
 163     }
 164     return true;
 165   }
 166   /**
 167    * Do the final analysis on the merged output.
 168    *
 169    * @return true on success
 170    */
 171   virtual Bool_t Finalize()
 172   {
 173     fDt = static_cast<TH1*>(fList->FindObject("dt"));
 174
 175     fFMD1i.Finalize(fList, fResults, fDt);
 176     fFMD2i.Finalize(fList, fResults, fDt);
 177     fFMD2o.Finalize(fList, fResults, fDt);
 178     fFMD3i.Finalize(fList, fResults, fDt);
 179     fFMD3o.Finalize(fList, fResults, fDt);
 180     return true;
 181   }
 182   /**
 183    * Get a reference to the event inspector. User must override this
 184    * to return proper object
 185    *
 186    * @return Reference to the event inspector
 187    */
 188   virtual AliFMDEventInspector& GetEventInspector() { return fEventInspector; }
 189   /**
 190    * Get a reference to the event inspector. User must override this
 191    * to return proper object
 192    *
 193    * @return Reference to the event inspector
 194    */
 195   virtual const AliFMDEventInspector& GetEventInspector() const
 196   {
 197     return fEventInspector;
 198   }
 199   /**
 200    * Read the map from timestamp to time-to-previous event
 201    *
 202    * @param filename File to read the map from
 203    *
 204    * @return true on success, false otherwise
 205    */
 206   Bool_t ReadMap(const char* filename)
 207   {
 208     if (gSystem->AccessPathName(filename, kReadPermission)) {
 209       // TSystem::AccessPathName returns false if file can be accessed!
 210       Error("ReadMap", "File \"%s\" cannot be open for reading", filename);
 211       return false;
 212     }
 213     Printf("Opening \"%s\" ...", filename);
 214     TFile* file = TFile::Open(filename, "READ");
 215     if (!file) {
 216       Error("ReadMap", "Failed to open file \"%s\"", filename);
 217       return false;
 218     }
 219     Printf("Opened \"%s\" ...", filename);
 220     TObject* o = file->Get("eventTimeMap");
 221     if (!o) {
 222       Error("ReadMap", "Failed to get \"eventTimeMap\" from %s", filename);
 223       return false;
 224     }
 225     Printf("Got object \"eventTimeMap\" from \"%s\" ...", filename);
 226     if (!o->IsA()->InheritsFrom(EventTimeMap::Class())) {
 227       Error("ReadMap", "Object \"%s\" is not an EventTimeMap, but a %s",
 228             o->GetName(), o->ClassName());
 229       return false;
 230     }
 231     Printf("Set the \"eventTimeMap\" to use");
 232     if (fMap) {
 233       delete fMap;
 234       fMap = 0;
 235     }
 236     fMap = static_cast<EventTimeMap*>(o);
 237     file->Close();
 238     return true;
 239   }
 240   /**
 241    * Create and connect the task
 242    *
 243    * @param mapfile File name of file containing timestamp map
 244    *
 245    * @return true on connect
 246    */
 247   static Bool_t Create(const char* mapfile)
 248   {
 249     ELossTimeTask* task = new ELossTimeTask("elossTime");
 250     if (!task->ReadMap(mapfile)) return false;
 251     task->Connect();
 252     return true;
 253   }
 254 protected:
 255   /**
 256    * Dummy copy constructor
 257    *
 258    * @param o Object to copy from
 259    */
 260   ELossTimeTask(const ELossTimeTask& o) : AliBaseESDTask(o) {}
 261   /** Our event inspector */
 262   AliFMDEventInspector fEventInspector;
 263   /**
 264    * Structure to hold per-ring histograms
 265    */
 266   struct RingHistos : public AliForwardUtil::RingHistos
 267   {
 268     /**
 269      * Constructor - for I/O only
 270      */
 271     RingHistos()
 272       : AliForwardUtil::RingHistos(),
 273         fDtVsELoss(0)
 274     {}
 275     /**
 276      * Constructor
 277      *
 278      * @param d detector number
 279      * @param r ring identifier
 280      */
 281     RingHistos(UShort_t d, Char_t r)
 282       : AliForwardUtil::RingHistos(d,r),
 283         fDtVsELoss(0)
 284     {}
 285     /**
 286      * Book histograms
 287      *
 288      * @param dir Parent list to add our list to
 289      * @param dt  Histogram of time differences
 290      *
 291      * @return true on success
 292      */
 293     Bool_t Book(TList* dir, TH1* dt)
 294     {
 295       TList* l = DefineOutputList(dir);
 296       // Double_t dtBins[] = { 0, 1, 2e3, 5e5, 1.1e6, 5e6, 1e12, 1e20 };
 297       fDtVsELoss = new TH2D("dtVsELoss",
 298                             Form("#Deltat vs #Delta/#Delta_{mip} - %s",
 299                                  GetName()), 450, 0, 15,
 300                             dt->GetXaxis()->GetNbins(),
 301                             dt->GetXaxis()->GetXmin(),
 302                             dt->GetXaxis()->GetXmax());
 303       fDtVsELoss->SetXTitle("#Delta/#Delta_{mip}");
 304       fDtVsELoss->SetYTitle("log_{10}(#Deltat) [ns]");
 305       fDtVsELoss->SetMarkerColor(Color());
 306       fDtVsELoss->SetDirectory(0);
 307       l->Add(fDtVsELoss);
 308
 309       return true;
 310     }
 311     /**
 312      * Process a single event
 313      *
 314      * @param fmd   FMD ESD data
 315      * @param logDt Logarithm (base 10) of Time-to-previous event
 316      */
 317     void Event(AliESDFMD& fmd, Double_t logDt)
 318     {
 319       const UShort_t nSec = NSector();
 320       const UShort_t nStr = NStrip();
 321       for (UShort_t s = 0; s < nSec; s++) {
 322         for (UShort_t t = 0; t < nStr; t++) {
 323           Float_t mult = fmd.Multiplicity(fDet,fRing,s,t);
 324           if(mult == AliESDFMD::kInvalidMult || mult <= 0)
 325             continue;
 326           fDtVsELoss->Fill(mult, logDt);
 327         }
 328       }
 329     }
 330     void Finalize(const TList* input, TList* output, TH1* dt)
 331     {
 332       TList* in  = GetOutputList(input);
 333       TList* out = DefineOutputList(output);
 334
 335       fDtVsELoss = static_cast<TH2*>(in->FindObject("dtVsELoss"));
 336
 337       TH2* dtVsELoss = static_cast<TH2*>(fDtVsELoss->Clone());
 338       dtVsELoss->SetDirectory(0);
 339       dtVsELoss->SetZTitle("1/N_{events}");
 340       dtVsELoss->Reset();
 341       for (UShort_t j = 1; j <= dt->GetNbinsX(); j++) {
 342         Double_t norm = dt->GetBinContent(j);
 343         if (norm <= 0) {
 344           // Warning("Finalize", "Bin %d empty in dT", j);
 345           continue;
 346         }
 347         for (UShort_t i = 1; i <= fDtVsELoss->GetNbinsX(); i++) {
 348           Double_t c = fDtVsELoss->GetBinContent(i, j);
 349           Double_t e = fDtVsELoss->GetBinError(i, j);
 350           dtVsELoss->SetBinContent(i,j,c/norm);
 351           dtVsELoss->SetBinError(i,j,e/norm);
 352         }
 353       }
 354       out->Add(dtVsELoss);
 355     }
 356     /** Our histogram */
 357     TH2* fDtVsELoss;
 358     ClassDef(RingHistos,1)
 359   };
 360   /** Container of FMD1i histograms */
 361   RingHistos    fFMD1i;
 362   /** Container of FMD2i histograms */
 363   RingHistos    fFMD2i;
 364   /** Container of FMD2o histograms */
 365   RingHistos    fFMD2o;
 366   /** Container of FMD3i histograms */
 367   RingHistos    fFMD3i;
 368   /** Container of FMD3o histograms */
 369   RingHistos    fFMD3o;
 370   /** Map from timestamp to time-to-previous event*/
 371   EventTimeMap* fMap;
 372   /** Distribution of log10(dt) */
 373   TH1* fDt;
 374   ClassDef(ELossTimeTask,2)
 375 };
 376 //
 377 // EOF
 378 //