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 "AliFMD.h"                             // ALIFMD_H
  47 #include "AliFMDDigit.h"                        // ALIFMDDIGIT_H
  48 #include "AliFMDParticles.h"                    // ALIFMDPARTICLES_H
  49 #include "AliFMDReconstructor.h"                // ALIFMDRECONSTRUCTOR_H
  50 #include "AliAltroBuffer.h"                     // ALIALTROBUFFER_H
  51 #include "AliLog.h"                             // ALILOG_H
  52 #include "AliRun.h"                             // ALIRUN_H
  53 #include "AliRunLoader.h"                       // ALIRUNLOADER_H
  54 #include "AliLoader.h"                          // ALILOADER_H
  55 #include "AliHeader.h"                          // ALIHEADER_H
  56 #include "AliGenEventHeader.h"                  // ALIGENEVENTHEADER_H
  57 #include "AliFMDRawStream.h"                    // ALIFMDRAWSTREAM_H
  58 #include "AliFMDRawReader.h"                    // ALIFMDRAWREADER_H
  59 #include "AliRawReader.h"                       // ALIRAWREADER_H
  60 #include "AliFMDReconstructionAlgorithm.h"      // ALIFMDRECONSTRUCTIONALGORITHM_H
  61
  62 //____________________________________________________________________
  63 ClassImp(AliFMDReconstructor);
  64
  65 //____________________________________________________________________
  66 AliFMDReconstructor::AliFMDReconstructor()
  67   : AliReconstructor(),
  68     fDeltaEta(0),
  69     fDeltaPhi(0),
  70     fThreshold(0),
  71     fPedestal(0),
  72     fPedestalWidth(0),
  73     fPedestalFactor(0)
  74 {
  75   // Make a new FMD reconstructor object - default CTOR.
  76   SetDeltaEta();
  77   SetDeltaPhi();
  78   SetThreshold();
  79   SetPedestal();
  80
  81   fParticles = new TClonesArray("AliFMDParticles", 1000);
  82   fFMDLoader = 0;
  83   fRunLoader = 0;
  84   fFMD       = 0;
  85 }
  86
  87
  88 //____________________________________________________________________
  89 AliFMDReconstructor::AliFMDReconstructor(const AliFMDReconstructor& other)
  90   : AliReconstructor(),
  91     fDeltaEta(0),
  92     fDeltaPhi(0),
  93     fThreshold(0),
  94     fPedestal(0),
  95     fPedestalWidth(0),
  96     fPedestalFactor(0)
  97 {
  98   // Make a new FMD reconstructor object - default CTOR.
  99   SetDeltaEta(other.fDeltaEta);
 100   SetDeltaPhi(other.fDeltaPhi);
 101   SetThreshold(other.fThreshold);
 102   SetPedestal(other.fPedestal, other.fPedestalWidth, other.fPedestalFactor);
 103
 104   // fParticles = new TClonesArray("AliFMDParticles", 1000);
 105   fFMDLoader = other.fFMDLoader;
 106   fRunLoader = other.fRunLoader;
 107   fFMD       = other.fFMD;
 108 }
 109
 110
 111 //____________________________________________________________________
 112 AliFMDReconstructor&
 113 AliFMDReconstructor::operator=(const AliFMDReconstructor& other)
 114 {
 115   // Make a new FMD reconstructor object - default CTOR.
 116   SetDeltaEta(other.fDeltaEta);
 117   SetDeltaPhi(other.fDeltaPhi);
 118   SetThreshold(other.fThreshold);
 119   SetPedestal(other.fPedestal, other.fPedestalWidth);
 120
 121   // fParticles = new TClonesArray("AliFMDParticles", 1000);
 122   fFMDLoader = other.fFMDLoader;
 123   fRunLoader = other.fRunLoader;
 124   fFMD       = other.fFMD;
 125
 126   return *this;
 127 }
 128
 129 //____________________________________________________________________
 130 void
 131 AliFMDReconstructor::SetPedestal(Float_t mean, Float_t width, Float_t factor)
 132 {
 133   // Set the pedestal, and pedestal width
 134   fPedestal       = mean;
 135   fPedestalWidth  = width;
 136   fPedestalFactor = factor;
 137 }
 138
 139 //____________________________________________________________________
 140 void
 141 AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader,
 142                                  AliRawReader* rawReader) const
 143 {
 144   // Collects all digits in the same active volume into number of
 145   // particles
 146   //
 147   // Reconstruct number of particles in given group of pads for given
 148   // FMDvolume determined by numberOfVolume,
 149   // numberOfMinSector, numberOfMaxSector, numberOfMinRing,
 150   // numberOgMaxRing
 151   //
 152   // The reconstruction method is choosen based on the number of empty
 153   // strips.
 154   fParticles->Clear();
 155   if (!runLoader) {
 156     Error("Exec","Run Loader loader is NULL - Session not opened");
 157     return;
 158   }
 159   fRunLoader = runLoader;
 160   fFMDLoader = runLoader->GetLoader("FMDLoader");
 161   if (!fFMDLoader)
 162     Fatal("AliFMDReconstructor","Can not find FMD (loader) "
 163           "in specified event");
 164
 165   // Get the AliRun object
 166   if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
 167
 168   // Get the AliFMD object
 169   fFMD = static_cast<AliFMD*>(fRunLoader->GetAliRun()->GetDetector("FMD"));
 170   if (!fFMD) {
 171     AliError("Can not get FMD from gAlice");
 172     return;
 173   }
 174   fFMDLoader->LoadRecPoints("RECREATE");
 175
 176   if (!fRunLoader->TreeE())     fRunLoader->LoadHeader();
 177
 178   if (rawReader) {
 179     Int_t event = 0;
 180     while (rawReader->NextEvent()) {
 181       ProcessEvent(event, rawReader);
 182       event++;
 183     }
 184   }
 185   else {
 186     Int_t nEvents= Int_t(fRunLoader->TreeE()->GetEntries());
 187     for(Int_t event = 0; event < nEvents; event++)
 188       ProcessEvent(event, 0);
 189   }
 190
 191
 192   fFMDLoader->UnloadRecPoints();
 193   fFMDLoader = 0;
 194   fRunLoader = 0;
 195   fFMD       = 0;
 196 }
 197
 198 //____________________________________________________________________
 199 void
 200 AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader) const
 201 {
 202   // Collects all digits in the same active volume into number of
 203   // particles
 204   //
 205   // Reconstruct number of particles in given group of pads for given
 206   // FMDvolume determined by numberOfVolume,
 207   // numberOfMinSector, numberOfMaxSector, numberOfMinRing,
 208   // numberOgMaxRing
 209   //
 210   // The reconstruction method is choosen based on the number of empty
 211   // strips.
 212   Reconstruct(runLoader, 0);
 213 }
 214
 215
 216 //____________________________________________________________________
 217 void
 218 AliFMDReconstructor::ProcessEvent(Int_t event,
 219                                   AliRawReader* reader) const
 220 {
 221   // Process one event read from either a clones array or from a a raw
 222   // data reader.
 223   fRunLoader->GetEvent(event) ;
 224   //event z-vertex for correction eta-rad dependence
 225   AliHeader *header            = fRunLoader->GetHeader();
 226   if (!header) Warning("ProcessEvent", "no AliHeader found!");
 227   AliGenEventHeader* genHeader = (header ? header->GenEventHeader() : 0);
 228
 229   // Get the Z--coordinate from the event header
 230   TArrayF o(3);
 231   if (genHeader) genHeader->PrimaryVertex(o);
 232   else           Warning("ProcessEvent", "No GenEventHeader Found");
 233   fCurrentVertex = o.At(2);
 234
 235   // If the recontruction tree isn't loaded, load it
 236   if(fFMDLoader->TreeR()==0) fFMDLoader->MakeTree("R");
 237
 238   //Make branches to hold the reconstructed particles
 239   const Int_t kBufferSize = 16000;
 240   fFMDLoader->TreeR()->Branch("FMD", &fParticles, kBufferSize);
 241
 242   // Load or recreate the digits
 243   if (fFMDLoader->LoadDigits((reader ? "UPDATE" : "READ"))) {
 244     if (!reader) {
 245       Error("Exec","Error occured while loading digits. Exiting.");
 246       return;
 247     }
 248   }
 249   // Get the digits tree
 250   TTree* digitTree = fFMDLoader->TreeD();
 251   if (!digitTree) {
 252     if (!reader) {
 253       Error("Exec","Can not get Tree with Digits. "
 254             "Nothing to reconstruct - Exiting");
 255       return;
 256     }
 257     fFMDLoader->MakeTree("D");
 258     digitTree = fFMDLoader->TreeD();
 259
 260   }
 261   // Get the FMD branch holding the digits.
 262   TBranch *digitBranch = digitTree->GetBranch("FMD");
 263   TClonesArray* digits = fFMD->Digits();
 264   if (!digitBranch) {
 265     if (!reader) {
 266       Error("Exec", "No digit branch for the FMD found");
 267       return;
 268     }
 269     fFMD->MakeBranchInTree(digitTree, fFMD->GetName(), &(digits), 4000, 0);
 270   }
 271   if (!reader) digitBranch->SetAddress(&digits);
 272
 273   if  (reader) {
 274     AliFMDRawReader rawRead(fFMD, reader);
 275     // rawRead->SetSampleRate(fSampleRate);
 276     rawRead.Exec();
 277   }
 278   else {
 279     // Read the ADC values from a clones array.
 280     AliDebug(10, "Reading ADCs from Digits array");
 281     // read Digits, and reconstruct the particles
 282     if (!fFMDLoader->TreeD()->GetEvent(0)) return;
 283   }
 284
 285   TIter next(&fAlgorithms);
 286   AliFMDReconstructionAlgorithm* algorithm = 0;
 287   while ((algorithm = static_cast<AliFMDReconstructionAlgorithm*>(next())))
 288     algorithm->PreEvent();
 289
 290   ProcessDigits(digits);
 291
 292   next.Reset();
 293   algorithm = 0;
 294   while ((algorithm = static_cast<AliFMDReconstructionAlgorithm*>(next())))
 295     algorithm->PostEvent();
 296
 297   if (reader) {
 298     digitTree->Fill();
 299     fFMDLoader->WriteDigits("OVERWRITE");
 300   }
 301   fFMDLoader->UnloadDigits();
 302   fFMDLoader->TreeR()->Reset();
 303   fFMDLoader->TreeR()->Fill();
 304   fFMDLoader->WriteRecPoints("OVERWRITE");
 305 }
 306
 307 //____________________________________________________________________
 308 UShort_t
 309 AliFMDReconstructor::SubtractPedestal(AliFMDDigit* digit) const
 310 {
 311   // Member function to subtract the pedestal from a digit
 312   // This implementation does nothing, but a derived class could over
 313   // load this to subtract a pedestal that was given in a database or
 314   // something like that.
 315
 316   Int_t counts = 0;
 317   Float_t ped = fPedestal * fPedestalFactor * fPedestalWidth;
 318   if (digit->Count3() >= 0)      counts = digit->Count3();
 319   else if (digit->Count2() >= 0) counts = digit->Count2();
 320   else                           counts = digit->Count2();
 321   counts = TMath::Max(Int_t(counts - ped), 0);
 322   return  UShort_t(counts);
 323 }
 324
 325 //____________________________________________________________________
 326 void
 327 AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
 328 {
 329   Int_t nDigits = digits->GetEntries();
 330   for (Int_t i = 0; i < nDigits; i++) {
 331     AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
 332     AliFMDSubDetector* subDetector = 0;
 333     switch (digit->Detector()) {
 334     case 1: subDetector = fFMD->GetFMD1(); break;
 335     case 2: subDetector = fFMD->GetFMD2(); break;
 336     case 3: subDetector = fFMD->GetFMD3(); break;
 337     }
 338     if (!subDetector) {
 339       Warning("ProcessDigits", "Unknown detector: FMD%d" , digit->Detector());
 340       continue;
 341     }
 342
 343     AliFMDRing* ring  = 0;
 344     Float_t     ringZ = 0;
 345     switch(digit->Ring()) {
 346     case 'i':
 347     case 'I':
 348       ring  = subDetector->GetInner();
 349       ringZ = subDetector->GetInnerZ();
 350       break;
 351     case 'o':
 352     case 'O':
 353       ring  = subDetector->GetOuter();
 354       ringZ = subDetector->GetOuterZ();
 355       break;
 356     }
 357     if (!ring) {
 358       Warning("ProcessDigits", "Unknown ring: FMD%d%c", digit->Detector(),
 359               digit->Ring());
 360       break;
 361     }
 362
 363     Float_t  realZ    = fCurrentVertex + ringZ;
 364     Float_t  stripR   = ((ring->GetHighR() - ring->GetLowR()) / ring->GetNStrips()
 365                          * (digit->Strip() + .5) + ring->GetLowR());
 366     Float_t  theta    = TMath::ATan2(stripR, realZ);
 367     Float_t  phi      = (2 * TMath::Pi() / ring->GetNSectors()
 368                          * (digit->Sector() + .5));
 369     Float_t  eta      = -TMath::Log(TMath::Tan(theta / 2));
 370     UShort_t counts   = SubtractPedestal(digit);
 371
 372     TIter next(&fAlgorithms);
 373     AliFMDReconstructionAlgorithm* algorithm = 0;
 374     while ((algorithm = static_cast<AliFMDReconstructionAlgorithm*>(next())))
 375       algorithm->ProcessDigit(digit, eta, phi, counts);
 376   }
 377 }
 378
 379
 380 //____________________________________________________________________
 381 void
 382 AliFMDReconstructor::ReconstructFromCache() const
 383 {
 384   // Based on the information in the cache, do the reconstruction.
 385   Int_t nRecon = 0;
 386   // Loop over the detectors
 387   for (Int_t i = 1; i <= 3; i++) {
 388     AliFMDSubDetector* sub = 0;
 389     switch (i) {
 390     case 1: sub = fFMD->GetFMD1(); break;
 391     case 2: sub = fFMD->GetFMD2(); break;
 392     case 3: sub = fFMD->GetFMD3(); break;
 393     }
 394     if (!sub) continue;
 395
 396     // Loop over the rings in the detector
 397     for (Int_t j = 0; j < 2; j++) {
 398       Float_t     rZ = 0;
 399       AliFMDRing* r  = 0;
 400       switch (j) {
 401       case 0: r  = sub->GetInner(); rZ = sub->GetInnerZ(); break;
 402       case 1: r  = sub->GetOuter(); rZ = sub->GetOuterZ(); break;
 403       }
 404       if (!r) continue;
 405
 406       // Calculate low/high theta and eta
 407       // FIXME: Is this right?
 408       Float_t realZ    = fCurrentVertex + rZ;
 409       Float_t thetaOut = TMath::ATan2(r->GetHighR(), realZ);
 410       Float_t thetaIn  = TMath::ATan2(r->GetLowR(), realZ);
 411       Float_t etaOut   = - TMath::Log(TMath::Tan(thetaOut / 2));
 412       Float_t etaIn    = - TMath::Log(TMath::Tan(thetaIn / 2));
 413       if (TMath::Abs(etaOut) > TMath::Abs(etaIn)) {
 414         Float_t tmp = etaIn;
 415         etaIn       = etaOut;
 416         etaOut      = tmp;
 417       }
 418
 419       //-------------------------------------------------------------
 420       //
 421       // Here starts poisson method
 422       //
 423       // Calculate eta step per strip, number of eta steps, number of
 424       // phi steps, and check the sign of the eta increment
 425       Float_t stripEta = (Float_t(r->GetNStrips()) / (etaIn - etaOut));
 426       Int_t   nEta     = Int_t(TMath::Abs(etaIn - etaOut) / fDeltaEta);
 427       Int_t   nPhi     = Int_t(360. / fDeltaPhi);
 428       Float_t sign     = TMath::Sign(Float_t(1.), etaIn);
 429
 430       AliDebug(10, Form("FMD%d%c Eta range: %f, %f %d Phi steps",
 431                         sub->GetId(), r->GetId(), etaOut, etaIn, nPhi));
 432
 433       // Loop over relevant phi values
 434       for (Int_t p = 0; p < nPhi; p++) {
 435         Float_t  minPhi    = p * fDeltaPhi;
 436         Float_t  maxPhi    = minPhi + fDeltaPhi;
 437         UShort_t minSector = UShort_t(minPhi / 360) * r->GetNSectors();
 438         UShort_t maxSector = UShort_t(maxPhi / 360) * r->GetNSectors();
 439
 440         AliDebug(10, Form(" Now in phi range %f, %f (sectors %d,%d)",
 441                           minPhi, maxPhi, minSector, maxSector));
 442         // Loop over relevant eta values
 443         for (Int_t e = nEta; e >= 0; --e) {
 444           Float_t  maxEta   = etaIn  - sign * e * fDeltaEta;
 445           Float_t  minEta   = maxEta - sign * fDeltaEta;
 446           if (sign > 0)  minEta = TMath::Max(minEta, etaOut);
 447           else           minEta = TMath::Min(minEta, etaOut);
 448           Float_t  theta1   = 2 * TMath::ATan(TMath::Exp(-minEta));
 449           Float_t  theta2   = 2 * TMath::ATan(TMath::Exp(-maxEta));
 450           Float_t  minR     = TMath::Abs(realZ * TMath::Tan(theta2));
 451           Float_t  maxR     = TMath::Abs(realZ * TMath::Tan(theta1));
 452           UShort_t minStrip = UShort_t((etaIn - maxEta) * stripEta + 0.5);
 453           UShort_t maxStrip = UShort_t((etaIn - minEta) * stripEta + 0.5);
 454
 455           AliDebug(10, Form("  Now in eta range %f, %f (strips %d, %d)\n"
 456                             "    [radii %f, %f, thetas %f, %f, sign %d]",
 457                             minEta, maxEta, minStrip, maxStrip,
 458                             minR, maxR, theta1, theta2, sign));
 459
 460           // Count number of empty strips
 461           Int_t   emptyStrips = 0;
 462           for (Int_t sector = minSector; sector < maxSector; sector++)
 463             for (Int_t strip = minStrip; strip < maxStrip; strip++) emptyStrips++;
 464           // if (fAdcs(sub->GetId() - 1, r->GetId(), sector, strip)
 465           //     < fThreshold) emptyStrips++;
 466
 467           // The total number of strips
 468           Float_t nTotal = (maxSector - minSector) * (maxStrip - minStrip);
 469
 470           // Log ratio of empty to total number of strips
 471           AliDebug(10, Form("Lambda= %d / %d = %f",
 472                             emptyStrips, nTotal,
 473                             Float_t(emptyStrips) / nTotal));
 474
 475           Double_t lambda = (emptyStrips > 0 ?
 476                              - TMath::Log(Double_t(emptyStrips) / nTotal) :
 477                              1);
 478
 479           // The reconstructed number of particles is then given by
 480           Int_t reconstructed = Int_t(lambda * nTotal + 0.5);
 481
 482           // Add a AliFMDParticles to the reconstruction tree.
 483           new((*fParticles)[nRecon])
 484             AliFMDParticles(sub->GetId(), r->GetId(),
 485                             minSector, maxSector, minStrip, maxStrip,
 486                             minEta, maxEta, minPhi, maxPhi,
 487                             reconstructed, AliFMDParticles::kPoission);
 488           nRecon++;
 489         } // phi
 490       } // eta
 491     } // ring
 492   } // detector
 493 }
 494
 495
 496 //____________________________________________________________________
 497 void
 498 AliFMDReconstructor::FillESD(AliRunLoader* /*fRunLoader*/,
 499                              AliESD* /*esd*/) const
 500 {
 501   // nothing to be done
 502
 503 }
 504
 505 //____________________________________________________________________
 506 //
 507 // EOF
 508 //