FMD/AliFMDReconstructor.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //____________________________________________________________________
  19 //
  20 // This is a class that constructs AliFMDMult (reconstructed
  21 // multiplicity) from of Digits
  22 //
  23 // This class reads either digits from a TClonesArray or raw data from
  24 // a DDL file (or similar), and stores the read ADC counts in an
  25 // internal cache (fAdcs).
  26 //
  27 // From the cached values it then calculates the number of particles
  28 // that hit a region of the FMDs, as specified by the user.
  29 //
  30 // The reconstruction can be done in two ways: Either via counting the
  31 // number of empty strips (Poisson method), or by converting the ADC
  32 // signal to an energy deposition, and then dividing by the typical
  33 // energy loss of a particle.
  34 //
  35 //      +---------------------+       +---------------------+
  36 //      | AliFMDReconstructor |<>-----| AliFMDMultAlgorithm |
  37 //      +---------------------+       +---------------------+
  38 //                                               ^
  39 //                                               |
  40 //                                   +-----------+---------+
  41 //                                   |                     |
  42 //                         +-------------------+   +------------------+
  43 //                         | AliFMDMultPoisson |   | AliFMDMultNaiive |
  44 //                         +-------------------+   +------------------+
  45 //
  46 // AliFMDReconstructor acts as a manager class.  It contains a list of
  47 // AliFMDMultAlgorithm objects.  The call graph looks something like
  48 //
  49 //
  50 //       +----------------------+            +----------------------+
  51 //       | :AliFMDReconstructor |            | :AliFMDMultAlgorithm |
  52 //       +----------------------+            +----------------------+
  53 //                  |                                  |
  54 //    Reconstruct  +-+                                 |
  55 //    ------------>| |                         PreRun +-+
  56 //                 | |------------------------------->| |
  57 //                 | |                                +-+
  58 //                 | |-----+ (for each event)          |
  59 //                 | |     | *ProcessEvent             |
  60 //                 |+-+    |                           |
  61 //                 || |<---+                 PreEvent +-+
  62 //                 || |------------------------------>| |
  63 //                 || |                               +-+
  64 //                 || |-----+                          |
  65 //                 || |     | ProcessDigits            |
  66 //                 ||+-+    |                          |
  67 //                 ||| |<---+                          |
  68 //                 ||| |         *ProcessDigit(digit) +-+
  69 //                 ||| |----------------------------->| |
  70 //                 ||| |                              +-+
  71 //                 ||+-+                               |
  72 //                 || |                     PostEvent +-+
  73 //                 || |------------------------------>| |
  74 //                 || |                               +-+
  75 //                 |+-+                                |
  76 //                 | |                        PostRun +-+
  77 //                 | |------------------------------->| |
  78 //                 | |                                +-+
  79 //                 +-+                                 |
  80 //                  |                                  |
  81 //
  82 //
  83 //
  84 //-- Authors: Evgeny Karpechev(INR) and Alla Maevsksia
  85 //  Latest changes by Christian Holm Christensen <cholm@nbi.dk>
  86 //
  87 //
  88 //____________________________________________________________________
  89
  90 #include <AliLog.h>                        // ALILOG_H
  91 #include <AliRun.h>                        // ALIRUN_H
  92 #include <AliRunLoader.h>                  // ALIRUNLOADER_H
  93 #include <AliLoader.h>                     // ALILOADER_H
  94 #include <AliHeader.h>                     // ALIHEADER_H
  95 #include <AliRawReader.h>                  // ALIRAWREADER_H
  96 #include <AliGenEventHeader.h>             // ALIGENEVENTHEADER_H
  97 #include "AliFMD.h"                        // ALIFMD_H
  98 #include "AliFMDDigit.h"                   // ALIFMDDIGIT_H
  99 #include "AliFMDReconstructor.h"           // ALIFMDRECONSTRUCTOR_H
 100 #include "AliFMDRawStream.h"               // ALIFMDRAWSTREAM_H
 101 #include "AliFMDRawReader.h"               // ALIFMDRAWREADER_H
 102 #include "AliFMDMultAlgorithm.h"           // ALIFMDMULTALGORITHM_H
 103 #include "AliFMDMultPoisson.h"             // ALIFMDMULTPOISSON_H
 104 #include "AliFMDMultNaiive.h"              // ALIFMDMULTNAIIVE_H
 105
 106 //____________________________________________________________________
 107 ClassImp(AliFMDReconstructor);
 108
 109 //____________________________________________________________________
 110 AliFMDReconstructor::AliFMDReconstructor()
 111   : AliReconstructor(),
 112     fPedestal(0),
 113     fPedestalWidth(0),
 114     fPedestalFactor(0)
 115 {
 116   // Make a new FMD reconstructor object - default CTOR.
 117   SetPedestal();
 118
 119   fFMDLoader = 0;
 120   fRunLoader = 0;
 121   fFMD       = 0;
 122   fAlgorithms.Add(new AliFMDMultNaiive);
 123   fAlgorithms.Add(new AliFMDMultPoisson);
 124 }
 125
 126
 127 //____________________________________________________________________
 128 AliFMDReconstructor::AliFMDReconstructor(const AliFMDReconstructor& other)
 129   : AliReconstructor(),
 130     fPedestal(0),
 131     fPedestalWidth(0),
 132     fPedestalFactor(0)
 133 {
 134   // Copy constructor
 135   SetPedestal(other.fPedestal, other.fPedestalWidth, other.fPedestalFactor);
 136
 137   fFMDLoader = other.fFMDLoader;
 138   fRunLoader = other.fRunLoader;
 139   fFMD       = other.fFMD;
 140
 141   fAlgorithms.Delete();
 142   TIter next(&(other.fAlgorithms));
 143   AliFMDMultAlgorithm* algorithm = 0;
 144   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
 145     fAlgorithms.Add(algorithm);
 146   fAlgorithms.SetOwner(kFALSE);
 147 }
 148
 149
 150 //____________________________________________________________________
 151 AliFMDReconstructor&
 152 AliFMDReconstructor::operator=(const AliFMDReconstructor& other)
 153 {
 154   // Assignment operator
 155   SetPedestal(other.fPedestal, other.fPedestalWidth, other.fPedestalFactor);
 156
 157   fFMDLoader = other.fFMDLoader;
 158   fRunLoader = other.fRunLoader;
 159   fFMD       = other.fFMD;
 160
 161   fAlgorithms.Delete();
 162   TIter next(&(other.fAlgorithms));
 163   AliFMDMultAlgorithm* algorithm = 0;
 164   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
 165     fAlgorithms.Add(algorithm);
 166   fAlgorithms.SetOwner(kFALSE);
 167
 168   return *this;
 169 }
 170
 171 //____________________________________________________________________
 172 AliFMDReconstructor::~AliFMDReconstructor()
 173 {
 174   // Destructor
 175   fAlgorithms.Delete();
 176 }
 177
 178 //____________________________________________________________________
 179 void
 180 AliFMDReconstructor::SetPedestal(Float_t mean, Float_t width, Float_t factor)
 181 {
 182   // Set the pedestal, and pedestal width
 183   fPedestal       = mean;
 184   fPedestalWidth  = width;
 185   fPedestalFactor = factor;
 186 }
 187
 188 //____________________________________________________________________
 189 void
 190 AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader,
 191                                  AliRawReader* rawReader) const
 192 {
 193   // Collects all digits in the same active volume into number of
 194   // particles
 195   //
 196   // Reconstruct number of particles in given group of pads for given
 197   // FMDvolume determined by numberOfVolume,
 198   // numberOfMinSector, numberOfMaxSector, numberOfMinRing,
 199   // numberOgMaxRing
 200   //
 201   // The reconstruction method is choosen based on the number of empty
 202   // strips.
 203   if (!runLoader) {
 204     Error("Exec","Run Loader loader is NULL - Session not opened");
 205     return;
 206   }
 207   fRunLoader = runLoader;
 208   fFMDLoader = runLoader->GetLoader("FMDLoader");
 209   if (!fFMDLoader)
 210     Fatal("AliFMDReconstructor","Can not find FMD (loader) "
 211           "in specified event");
 212
 213   // Get the AliRun object
 214   if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
 215
 216   // Get the AliFMD object
 217   fFMD = static_cast<AliFMD*>(fRunLoader->GetAliRun()->GetDetector("FMD"));
 218   if (!fFMD) {
 219     AliError("Can not get FMD from gAlice");
 220     return;
 221   }
 222   fFMDLoader->LoadRecPoints("RECREATE");
 223
 224   if (!fRunLoader->TreeE())     fRunLoader->LoadHeader();
 225
 226   TIter next(&fAlgorithms);
 227   AliFMDMultAlgorithm* algorithm = 0;
 228   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
 229     algorithm->PreRun(fFMD);
 230
 231   if (rawReader) {
 232     Int_t event = 0;
 233     while (rawReader->NextEvent()) {
 234       ProcessEvent(event, rawReader);
 235       event++;
 236     }
 237   }
 238   else {
 239     Int_t nEvents= Int_t(fRunLoader->TreeE()->GetEntries());
 240     for(Int_t event = 0; event < nEvents; event++)
 241       ProcessEvent(event, 0);
 242   }
 243
 244   next.Reset();
 245   algorithm = 0;
 246   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
 247     algorithm->PostRun();
 248
 249   fFMDLoader->UnloadRecPoints();
 250   fFMDLoader = 0;
 251   fRunLoader = 0;
 252   fFMD       = 0;
 253 }
 254
 255 //____________________________________________________________________
 256 void
 257 AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader) const
 258 {
 259   // Collects all digits in the same active volume into number of
 260   // particles
 261   //
 262   // Reconstruct number of particles in given group of pads for given
 263   // FMDvolume determined by numberOfVolume,
 264   // numberOfMinSector, numberOfMaxSector, numberOfMinRing,
 265   // numberOgMaxRing
 266   //
 267   // The reconstruction method is choosen based on the number of empty
 268   // strips.
 269   Reconstruct(runLoader, 0);
 270 }
 271
 272
 273 //____________________________________________________________________
 274 void
 275 AliFMDReconstructor::ProcessEvent(Int_t event,
 276                                   AliRawReader* reader) const
 277 {
 278   // Process one event read from either a clones array or from a a raw
 279   // data reader.
 280   fRunLoader->GetEvent(event) ;
 281   //event z-vertex for correction eta-rad dependence
 282   AliHeader *header            = fRunLoader->GetHeader();
 283   if (!header) Warning("ProcessEvent", "no AliHeader found!");
 284   AliGenEventHeader* genHeader = (header ? header->GenEventHeader() : 0);
 285
 286   // Get the Z--coordinate from the event header
 287   TArrayF o(3);
 288   if (genHeader) genHeader->PrimaryVertex(o);
 289   else           Warning("ProcessEvent", "No GenEventHeader Found");
 290   fCurrentVertex = o.At(2);
 291
 292   // If the recontruction tree isn't loaded, load it
 293   if(fFMDLoader->TreeR()==0) fFMDLoader->MakeTree("R");
 294
 295   // Load or recreate the digits
 296   if (fFMDLoader->LoadDigits((reader ? "UPDATE" : "READ"))) {
 297     if (!reader) {
 298       Error("Exec","Error occured while loading digits. Exiting.");
 299       return;
 300     }
 301   }
 302   // Get the digits tree
 303   TTree* digitTree = fFMDLoader->TreeD();
 304   if (!digitTree) {
 305     if (!reader) {
 306       Error("Exec","Can not get Tree with Digits. "
 307             "Nothing to reconstruct - Exiting");
 308       return;
 309     }
 310     fFMDLoader->MakeTree("D");
 311     digitTree = fFMDLoader->TreeD();
 312
 313   }
 314   // Get the FMD branch holding the digits.
 315   TBranch *digitBranch = digitTree->GetBranch("FMD");
 316   TClonesArray* digits = fFMD->Digits();
 317   if (!digitBranch) {
 318     if (!reader) {
 319       Error("Exec", "No digit branch for the FMD found");
 320       return;
 321     }
 322     fFMD->MakeBranchInTree(digitTree, fFMD->GetName(), &(digits), 4000, 0);
 323   }
 324   if (!reader) digitBranch->SetAddress(&digits);
 325
 326   if  (reader) {
 327     AliFMDRawReader rawRead(fFMD, reader);
 328     AliDebug(10, Form("Making raw reader with sample rate: %d",
 329                       fFMD->GetSampleRate()));
 330     rawRead.SetSampleRate(fFMD->GetSampleRate());
 331     rawRead.Exec();
 332   }
 333   else {
 334     // Read the ADC values from a clones array.
 335     AliDebug(10, "Reading ADCs from Digits array");
 336     // read Digits, and reconstruct the particles
 337     if (!fFMDLoader->TreeD()->GetEvent(0)) return;
 338   }
 339
 340   TIter next(&fAlgorithms);
 341   AliFMDMultAlgorithm* algorithm = 0;
 342   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
 343     algorithm->PreEvent(fFMDLoader->TreeR(), fCurrentVertex);
 344
 345   ProcessDigits(digits);
 346
 347   next.Reset();
 348   algorithm = 0;
 349   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
 350     algorithm->PostEvent();
 351
 352   if (reader) {
 353     digitTree->Fill();
 354     fFMDLoader->WriteDigits("OVERWRITE");
 355   }
 356   fFMDLoader->UnloadDigits();
 357   fFMDLoader->TreeR()->Reset();
 358   fFMDLoader->TreeR()->Fill();
 359   fFMDLoader->WriteRecPoints("OVERWRITE");
 360 }
 361
 362 //____________________________________________________________________
 363 UShort_t
 364 AliFMDReconstructor::SubtractPedestal(AliFMDDigit* digit) const
 365 {
 366   // Member function to subtract the pedestal from a digit
 367   // This implementation does nothing, but a derived class could over
 368   // load this to subtract a pedestal that was given in a database or
 369   // something like that.
 370
 371   Int_t counts = 0;
 372   Float_t ped = fPedestal + fPedestalFactor * fPedestalWidth;
 373   if (digit->Count3() > 0)      counts = digit->Count3();
 374   else if (digit->Count2() > 0) counts = digit->Count2();
 375   else                          counts = digit->Count1();
 376   counts = TMath::Max(Int_t(counts - ped), 0);
 377   return  UShort_t(counts);
 378 }
 379
 380 //____________________________________________________________________
 381 void
 382 AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
 383 {
 384   Int_t nDigits = digits->GetEntries();
 385   for (Int_t i = 0; i < nDigits; i++) {
 386     AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
 387     AliFMDSubDetector* subDetector = 0;
 388     switch (digit->Detector()) {
 389     case 1: subDetector = fFMD->GetFMD1(); break;
 390     case 2: subDetector = fFMD->GetFMD2(); break;
 391     case 3: subDetector = fFMD->GetFMD3(); break;
 392     }
 393     if (!subDetector) {
 394       Warning("ProcessDigits", "Unknown detector: FMD%d" , digit->Detector());
 395       continue;
 396     }
 397
 398     AliFMDRing* ring  = 0;
 399     Float_t     ringZ = 0;
 400     switch(digit->Ring()) {
 401     case 'i':
 402     case 'I':
 403       ring  = subDetector->GetInner();
 404       ringZ = subDetector->GetInnerZ();
 405       break;
 406     case 'o':
 407     case 'O':
 408       ring  = subDetector->GetOuter();
 409       ringZ = subDetector->GetOuterZ();
 410       break;
 411     }
 412     if (!ring) {
 413       Warning("ProcessDigits", "Unknown ring: FMD%d%c", digit->Detector(),
 414               digit->Ring());
 415       break;
 416     }
 417
 418     Float_t  realZ    = fCurrentVertex + ringZ;
 419     Float_t  stripR   = ((ring->GetHighR() - ring->GetLowR())
 420                          / ring->GetNStrips() * (digit->Strip() + .5)
 421                          + ring->GetLowR());
 422     Float_t  theta    = TMath::ATan2(stripR, realZ);
 423     Float_t  phi      = (2 * TMath::Pi() / ring->GetNSectors()
 424                          * (digit->Sector() + .5));
 425     Float_t  eta      = -TMath::Log(TMath::Tan(theta / 2));
 426     UShort_t counts   = SubtractPedestal(digit);
 427
 428     TIter next(&fAlgorithms);
 429     AliFMDMultAlgorithm* algorithm = 0;
 430     while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
 431       algorithm->ProcessDigit(digit, eta, phi, counts);
 432   }
 433 }
 434
 435
 436 //____________________________________________________________________
 437 void
 438 AliFMDReconstructor::FillESD(AliRunLoader* /*fRunLoader*/,
 439                              AliESD* /*esd*/) const
 440 {
 441   // nothing to be done
 442
 443 }
 444
 445 //____________________________________________________________________
 446 //
 447 // EOF
 448 //