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 "AliFMDGeometry.h"                // ALIFMDGEOMETRY_H
  99 #include "AliFMDParameters.h"              // ALIFMDPARAMETERS_H
 100 #include "AliFMDDetector.h"                // ALIFMDDETECTOR_H
 101 #include "AliFMDRing.h"                    // ALIFMDRING_H
 102 #include "AliFMDDigit.h"                   // ALIFMDDIGIT_H
 103 #include "AliFMDReconstructor.h"           // ALIFMDRECONSTRUCTOR_H
 104 #include "AliFMDRawStream.h"               // ALIFMDRAWSTREAM_H
 105 #include "AliFMDRawReader.h"               // ALIFMDRAWREADER_H
 106 #include "AliFMDMultAlgorithm.h"           // ALIFMDMULTALGORITHM_H
 107 #include "AliFMDMultPoisson.h"             // ALIFMDMULTPOISSON_H
 108 #include "AliFMDMultNaiive.h"              // ALIFMDMULTNAIIVE_H
 109 #include "AliESD.h"                        // ALIESD_H
 110
 111 //____________________________________________________________________
 112 ClassImp(AliFMDReconstructor)
 113 #if 0
 114   ; // This is here to keep Emacs for indenting the next line
 115 #endif
 116
 117 //____________________________________________________________________
 118 AliFMDReconstructor::AliFMDReconstructor()
 119   : AliReconstructor(),
 120     fPedestal(0),
 121     fPedestalWidth(0),
 122     fPedestalFactor(0)
 123 {
 124   // Make a new FMD reconstructor object - default CTOR.
 125   AliFMDParameters* pars = AliFMDParameters::Instance();
 126   fPedestal       = pars->GetPedestal();
 127   fPedestalWidth  = pars->GetPedestalWidth();
 128   fPedestalFactor = pars->GetPedestalFactor();
 129
 130   fAlgorithms.Add(new AliFMDMultNaiive);
 131   fAlgorithms.Add(new AliFMDMultPoisson);
 132
 133   TIter next(&fAlgorithms);
 134   AliFMDMultAlgorithm* algorithm = 0;
 135   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
 136     algorithm->PreRun(0);
 137 }
 138
 139
 140 //____________________________________________________________________
 141 AliFMDReconstructor::AliFMDReconstructor(const AliFMDReconstructor& other)
 142   : AliReconstructor(),
 143     fPedestal(0),
 144     fPedestalWidth(0),
 145     fPedestalFactor(0)
 146 {
 147   // Copy constructor
 148   fPedestal       = other.fPedestal;
 149   fPedestalWidth  = other.fPedestalWidth;
 150   fPedestalFactor = other.fPedestalFactor;
 151
 152
 153   fAlgorithms.Delete();
 154   TIter next(&(other.fAlgorithms));
 155   AliFMDMultAlgorithm* algorithm = 0;
 156   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
 157     fAlgorithms.Add(algorithm);
 158   fAlgorithms.SetOwner(kFALSE);
 159 }
 160
 161
 162 //____________________________________________________________________
 163 AliFMDReconstructor&
 164 AliFMDReconstructor::operator=(const AliFMDReconstructor& other)
 165 {
 166   // Assignment operator
 167   fPedestal       = other.fPedestal;
 168   fPedestalWidth  = other.fPedestalWidth;
 169   fPedestalFactor = other.fPedestalFactor;
 170
 171   fAlgorithms.Delete();
 172   TIter next(&(other.fAlgorithms));
 173   AliFMDMultAlgorithm* algorithm = 0;
 174   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
 175     fAlgorithms.Add(algorithm);
 176   fAlgorithms.SetOwner(kFALSE);
 177
 178   return *this;
 179 }
 180
 181 //____________________________________________________________________
 182 AliFMDReconstructor::~AliFMDReconstructor()
 183 {
 184   // Destructor
 185   fAlgorithms.Delete();
 186 }
 187
 188 //____________________________________________________________________
 189 void
 190 AliFMDReconstructor::Init(AliRunLoader* runLoader)
 191 {
 192   // Initialize the reconstructor
 193   AliDebug(1, Form("Init called with runloader 0x%x", runLoader));
 194   fCurrentVertex = 0;
 195   if (!runLoader) {
 196     Warning("Init", "No run loader");
 197     return;
 198   }
 199   AliHeader* header = runLoader->GetHeader();
 200   if (!header) {
 201     Warning("Init", "No header");
 202     return;
 203   }
 204   AliGenEventHeader* eventHeader = header->GenEventHeader();
 205   if (eventHeader) {
 206     TArrayF vtx;
 207     eventHeader->PrimaryVertex(vtx);
 208     fCurrentVertex = vtx[2];
 209     AliDebug(1, Form("Primary vertex Z coordinate for event # %d/%d is %f",
 210                      header->GetRun(), header->GetEvent(), fCurrentVertex));
 211     Warning("Init", "no generator event header");
 212   }
 213   else {
 214     Warning("Init", "No generator event header - "
 215             "perhaps we get the vertex from ESD?");
 216   }
 217   // Get the ESD tree
 218   SetESD(new AliESD);
 219 }
 220
 221 //____________________________________________________________________
 222 void
 223 AliFMDReconstructor::ConvertDigits(AliRawReader* reader,
 224                                    TTree* digitsTree) const
 225 {
 226   // Convert Raw digits to AliFMDDigit's in a tree
 227   AliFMDRawReader rawRead(reader, digitsTree);
 228   // rawRead.SetSampleRate(fFMD->GetSampleRate());
 229   rawRead.Exec();
 230 }
 231
 232 //____________________________________________________________________
 233 void
 234 AliFMDReconstructor::Reconstruct(TTree* digitsTree,
 235                                  TTree* clusterTree) const
 236 {
 237   // Reconstruct event from digits in tree
 238   // Get the FMD branch holding the digits.
 239   // FIXME: The vertex may not be known yet, so we may have to move
 240   // some of this to FillESD.
 241   AliDebug(1, "Reconstructing from digits in a tree");
 242 #if 0
 243   if (fESD) {
 244     // const AliESDVertex* vertex = fESD->GetVertex();
 245     // if (vertex) {
 246     //   AliDebug(1, Form("Got vertex from ESD: %f", vertex->GetZv()));
 247     //   fCurrentVertex = vertex->GetZv();
 248     // }
 249   }
 250 #endif
 251   TBranch *digitBranch = digitsTree->GetBranch("FMD");
 252   if (!digitBranch) {
 253     Error("Exec", "No digit branch for the FMD found");
 254     return;
 255   }
 256   TClonesArray* digits = new TClonesArray("AliFMDDigit");
 257   digitBranch->SetAddress(&digits);
 258
 259   TIter next(&fAlgorithms);
 260   AliFMDMultAlgorithm* algorithm = 0;
 261   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
 262     algorithm->PreEvent(clusterTree, fCurrentVertex);
 263   digitBranch->GetEntry(0);
 264
 265   ProcessDigits(digits);
 266
 267   next.Reset();
 268   algorithm = 0;
 269   while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
 270     algorithm->PostEvent();
 271   clusterTree->Fill();
 272   digits->Delete();
 273   delete digits;
 274 }
 275
 276 //____________________________________________________________________
 277 void
 278 AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
 279 {
 280   // For each digit, find the pseudo rapdity, azimuthal angle, and
 281   // number of corrected ADC counts, and pass it on to the algorithms
 282   // used.
 283   Int_t nDigits = digits->GetEntries();
 284   AliDebug(1, Form("Got %d digits", nDigits));
 285   for (Int_t i = 0; i < nDigits; i++) {
 286     AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
 287     AliFMDGeometry* fmd = AliFMDGeometry::Instance();
 288     AliFMDDetector* subDetector = fmd->GetDetector(digit->Detector());
 289     if (!subDetector) {
 290       Warning("ProcessDigits", "Unknown detector: FMD%d" , digit->Detector());
 291       continue;
 292     }
 293
 294     AliFMDRing* ring  = subDetector->GetRing(digit->Ring());
 295     Float_t     ringZ = subDetector->GetRingZ(digit->Ring());
 296     if (!ring) {
 297       Warning("ProcessDigits", "Unknown ring: FMD%d%c", digit->Detector(),
 298               digit->Ring());
 299       break;
 300     }
 301
 302     Float_t  realZ    = fCurrentVertex + ringZ;
 303     Float_t  stripR   = ((ring->GetHighR() - ring->GetLowR())
 304                          / ring->GetNStrips() * (digit->Strip() + .5)
 305                          + ring->GetLowR());
 306     Float_t  theta    = TMath::ATan2(stripR, realZ);
 307     Float_t  phi      = (2 * TMath::Pi() / ring->GetNSectors()
 308                          * (digit->Sector() + .5));
 309     Float_t  eta      = -TMath::Log(TMath::Tan(theta / 2));
 310     UShort_t counts   = SubtractPedestal(digit);
 311
 312     TIter next(&fAlgorithms);
 313     AliFMDMultAlgorithm* algorithm = 0;
 314     while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next())))
 315       algorithm->ProcessDigit(digit, eta, phi, counts);
 316   }
 317 }
 318
 319 //____________________________________________________________________
 320 UShort_t
 321 AliFMDReconstructor::SubtractPedestal(AliFMDDigit* digit) const
 322 {
 323   // Member function to subtract the pedestal from a digit
 324   // This implementation does nothing, but a derived class could over
 325   // load this to subtract a pedestal that was given in a database or
 326   // something like that.
 327
 328   Int_t counts = 0;
 329   Float_t ped = fPedestal + fPedestalFactor * fPedestalWidth;
 330   if (digit->Count3() > 0)      counts = digit->Count3();
 331   else if (digit->Count2() > 0) counts = digit->Count2();
 332   else                          counts = digit->Count1();
 333   counts = TMath::Max(Int_t(counts - ped), 0);
 334   return  UShort_t(counts);
 335 }
 336
 337 //____________________________________________________________________
 338 void
 339 AliFMDReconstructor::FillESD(TTree*  /* digitsTree */,
 340                              TTree*  /* clusterTree */,
 341                              AliESD* /* esd*/) const
 342 {
 343   // nothing to be done
 344   // FIXME: The vertex may not be known when Reconstruct is executed,
 345   // so we may have to move some of that member function here.
 346 #if 0
 347   TClonesArray* multStrips  = 0;
 348   TClonesArray* multRegions = 0;
 349   TTree*        treeR  = fmdLoader->TreeR();
 350   TBranch*      branchRegions = treeR->GetBranch("FMDPoisson");
 351   TBranch*      branchStrips  = treeR->GetBranch("FMDNaiive");
 352   branchRegions->SetAddress(&multRegions);
 353   branchStrips->SetAddress(&multStrips);
 354
 355   Int_t total = 0;
 356   Int_t nEntries  = clusterTree->GetEntries();
 357   for (Int_t entry = 0; entry < nEntries; entry++) {
 358     AliDebug(5, Form("Entry # %d in cluster tree", entry));
 359     treeR->GetEntry(entry);
 360
 361
 362     Int_t nMults = multRegions->GetLast();
 363     for (Int_t i = 0; i <= nMults; i++) {
 364       AliFMDMultRegion* multR =
 365         static_cast<AliFMDMultRegion*>(multRegions->UncheckedAt(i));
 366       Int_t nParticles=multR->Particles();
 367       if (i>=0 && i<=13)   hEtaPoissonI1->AddBinContent(i+1,nParticles);
 368       if (i>=14 && i<=27 ) hEtaPoissonI2->AddBinContent(i-13,nParticles);
 369       if (i>=28 && i<=33 );
 370       if (i>=34 && i<=47 ) hEtaPoissonI3->AddBinContent(48-i,nParticles);
 371       if (i>=48 && i<=53)  hEtaPoissonO3->AddBinContent(54-i,nParticles);
 372     }
 373   }
 374 #endif
 375 }
 376
 377
 378 //____________________________________________________________________
 379 void
 380 AliFMDReconstructor::Reconstruct(AliRawReader*,TTree*) const
 381 {
 382   // Cannot be used.  See member function with same name but with 2
 383   // TTree arguments.   Make sure you do local reconstrucion
 384   AliError("MayNotUse");
 385 }
 386 //____________________________________________________________________
 387 void
 388 AliFMDReconstructor::Reconstruct(AliRunLoader*) const
 389 {
 390   // Cannot be used.  See member function with same name but with 2
 391   // TTree arguments.   Make sure you do local reconstrucion
 392   AliError("MayNotUse");
 393 }
 394 //____________________________________________________________________
 395 void
 396 AliFMDReconstructor::Reconstruct(AliRunLoader*, AliRawReader*) const
 397 {
 398   // Cannot be used.  See member function with same name but with 2
 399   // TTree arguments.   Make sure you do local reconstrucion
 400   AliError("MayNotUse");
 401 }
 402 //____________________________________________________________________
 403 void
 404 AliFMDReconstructor::FillESD(AliRawReader*,TTree*,AliESD*) const
 405 {
 406   // Cannot be used.  See member function with same name but with 2
 407   // TTree arguments.   Make sure you do local reconstrucion
 408   AliError("MayNotUse");
 409 }
 410 //____________________________________________________________________
 411 void
 412 AliFMDReconstructor::FillESD(AliRunLoader*,AliESD*) const
 413 {
 414   // Cannot be used.  See member function with same name but with 2
 415   // TTree arguments.   Make sure you do local reconstrucion
 416   AliError("MayNotUse");
 417 }
 418 //____________________________________________________________________
 419 void
 420 AliFMDReconstructor::FillESD(AliRunLoader*,AliRawReader*,AliESD*) const
 421 {
 422   // Cannot be used.  See member function with same name but with 2
 423   // TTree arguments.   Make sure you do local reconstrucion
 424   AliError("MayNotUse");
 425 }
 426
 427 //____________________________________________________________________
 428 //
 429 // EOF
 430 //