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