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"
  47 #include "AliFMDDigit.h"
  48 #include "AliFMDParticles.h"
  49 #include "AliFMDReconstructor.h"
  50 #include "AliAltroBuffer.h"
  51 #include "AliLog.h"
  52 #include "AliRun.h"
  53 #include "AliRunLoader.h"
  54 #include "AliLoader.h"
  55 #include "AliHeader.h"
  56 #include "AliGenEventHeader.h"
  57 #include "AliFMDRawStream.h"
  58 #include "AliRawReader.h"
  59
  60 //____________________________________________________________________
  61 ClassImp(AliFMDReconstructor);
  62
  63 //____________________________________________________________________
  64 AliFMDReconstructor::AliFMDReconstructor()
  65   : AliReconstructor(),
  66     fAdcs(kMaxDetectors, kMaxRings, kMaxSectors, kMaxStrips),
  67     fDeltaEta(0),
  68     fDeltaPhi(0),
  69     fThreshold(0),
  70     fPedestal(0),
  71     fPedestalWidth(0)
  72 {
  73   // Make a new FMD reconstructor object - default CTOR.
  74   SetDeltaEta();
  75   SetDeltaPhi();
  76   SetThreshold();
  77   SetPedestal();
  78
  79   fParticles = new TClonesArray("AliFMDParticles", 1000);
  80   fFMDLoader = 0;
  81   fRunLoader = 0;
  82   fFMD       = 0;
  83 }
  84
  85
  86 //____________________________________________________________________
  87 AliFMDReconstructor::AliFMDReconstructor(const AliFMDReconstructor& other)
  88   : AliReconstructor(),
  89     fAdcs(kMaxDetectors, kMaxRings, kMaxSectors, kMaxStrips),
  90     fDeltaEta(0),
  91     fDeltaPhi(0),
  92     fThreshold(0),
  93     fPedestal(0),
  94     fPedestalWidth(0)
  95 {
  96   // Make a new FMD reconstructor object - default CTOR.
  97   SetDeltaEta(other.fDeltaEta);
  98   SetDeltaPhi(other.fDeltaPhi);
  99   SetThreshold(other.fThreshold);
 100   SetPedestal(other.fPedestal, other.fPedestalWidth);
 101
 102   // fParticles = new TClonesArray("AliFMDParticles", 1000);
 103   fFMDLoader = other.fFMDLoader;
 104   fRunLoader = other.fRunLoader;
 105   fFMD       = other.fFMD;
 106 }
 107
 108
 109 //____________________________________________________________________
 110 AliFMDReconstructor&
 111 AliFMDReconstructor::operator=(const AliFMDReconstructor& other)
 112 {
 113   // Make a new FMD reconstructor object - default CTOR.
 114   SetDeltaEta(other.fDeltaEta);
 115   SetDeltaPhi(other.fDeltaPhi);
 116   SetThreshold(other.fThreshold);
 117   SetPedestal(other.fPedestal, other.fPedestalWidth);
 118
 119   // fParticles = new TClonesArray("AliFMDParticles", 1000);
 120   fFMDLoader = other.fFMDLoader;
 121   fRunLoader = other.fRunLoader;
 122   fFMD       = other.fFMD;
 123
 124   return *this;
 125 }
 126
 127 //____________________________________________________________________
 128 void
 129 AliFMDReconstructor::SetPedestal(Float_t mean, Float_t width)
 130 {
 131   // Set the pedestal, and pedestal width
 132   fPedestal      = mean;
 133   fPedestalWidth = width;
 134 }
 135
 136 //____________________________________________________________________
 137 void
 138 AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader,
 139                                  AliRawReader* rawReader) const
 140 {
 141   // Collects all digits in the same active volume into number of
 142   // particles
 143   //
 144   // Reconstruct number of particles in given group of pads for given
 145   // FMDvolume determined by numberOfVolume,
 146   // numberOfMinSector, numberOfMaxSector, numberOfMinRing,
 147   // numberOgMaxRing
 148   //
 149   // The reconstruction method is choosen based on the number of empty
 150   // strips.
 151   fParticles->Clear();
 152   if (!runLoader) {
 153     Error("Exec","Run Loader loader is NULL - Session not opened");
 154     return;
 155   }
 156   fRunLoader = runLoader;
 157   fFMDLoader = runLoader->GetLoader("FMDLoader");
 158   if (!fFMDLoader)
 159     Fatal("AliFMDReconstructor","Can not find FMD (loader) "
 160           "in specified event");
 161
 162   // Get the AliRun object
 163   if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
 164
 165   // Get the AliFMD object
 166   fFMD = static_cast<AliFMD*>(fRunLoader->GetAliRun()->GetDetector("FMD"));
 167   if (!fFMD) {
 168     AliError("Can not get FMD from gAlice");
 169     return;
 170   }
 171   fFMDLoader->LoadRecPoints("RECREATE");
 172
 173   if (!fRunLoader->TreeE())     fRunLoader->LoadHeader();
 174
 175   TClonesArray* digits = fFMD->Digits();
 176   if (rawReader) {
 177     Int_t event = 0;
 178     while (rawReader->NextEvent()) {
 179       ProcessEvent(event, rawReader, digits);
 180       event++;
 181     }
 182   }
 183   else {
 184     Int_t nEvents= Int_t(fRunLoader->TreeE()->GetEntries());
 185     for(Int_t event = 0; event < nEvents; event++)
 186       ProcessEvent(event, 0, digits);
 187   }
 188
 189
 190   fFMDLoader->UnloadRecPoints();
 191   fFMDLoader = 0;
 192   fRunLoader = 0;
 193   fFMD       = 0;
 194 }
 195
 196 //____________________________________________________________________
 197 void
 198 AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader) const
 199 {
 200   // Collects all digits in the same active volume into number of
 201   // particles
 202   //
 203   // Reconstruct number of particles in given group of pads for given
 204   // FMDvolume determined by numberOfVolume,
 205   // numberOfMinSector, numberOfMaxSector, numberOfMinRing,
 206   // numberOgMaxRing
 207   //
 208   // The reconstruction method is choosen based on the number of empty
 209   // strips.
 210   Reconstruct(runLoader, 0);
 211 }
 212
 213
 214 //____________________________________________________________________
 215 void
 216 AliFMDReconstructor::ProcessEvent(Int_t event,
 217                                   AliRawReader* reader,
 218                                   TClonesArray* digits) const
 219 {
 220   // Process one event read from either a clones array or from a a raw
 221   // data reader.
 222   fRunLoader->GetEvent(event) ;
 223   //event z-vertex for correction eta-rad dependence
 224   AliHeader *header            = fRunLoader->GetHeader();
 225   if (!header)
 226     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   Float_t zVertex = o.At(2);
 233
 234   // If the recontruction tree isn't loaded, load it
 235   if(fFMDLoader->TreeR()==0) fFMDLoader->MakeTree("R");
 236
 237   //Make branches to hold the reconstructed particles
 238   const Int_t kBufferSize = 16000;
 239   fFMDLoader->TreeR()->Branch("FMD", &fParticles, kBufferSize);
 240
 241   // Load or recreate the digits
 242   if (fFMDLoader->LoadDigits((reader ? "UPDATE" : "READ"))) {
 243     if (!reader) {
 244       Error("Exec","Error occured while loading digits. Exiting.");
 245       return;
 246     }
 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   if (!digitBranch) {
 264     if (!reader) {
 265       Error("Exec", "No digit branch for the FMD found");
 266       return;
 267     }
 268     fFMD->MakeBranchInTree(digitTree, fFMD->GetName(), &(digits), 4000, 0);
 269   }
 270   if (!reader) digitBranch->SetAddress(&digits);
 271
 272   fEmptyStrips = 0;
 273   fTotalStrips = 0;
 274   Bool_t ok = kFALSE;
 275   if      (reader) ok = ReadAdcs(reader);
 276   else if (digits) ok = ReadAdcs(digits);
 277   if (!ok) return;
 278
 279   ReconstructFromCache(zVertex);
 280
 281   if (reader) {
 282     digitTree->Fill();
 283     fFMDLoader->WriteDigits("OVERWRITE");
 284   }
 285   fFMDLoader->UnloadDigits();
 286   fFMDLoader->TreeR()->Reset();
 287   fFMDLoader->TreeR()->Fill();
 288   fFMDLoader->WriteRecPoints("OVERWRITE");
 289 }
 290
 291 //____________________________________________________________________
 292 Bool_t
 293 AliFMDReconstructor::ReadAdcs(TClonesArray* digits) const
 294 {
 295   // Read the ADC values from a clones array.
 296   AliDebug(10, "Reading ADCs from Digits array");
 297   // read Digits, and reconstruct the particles
 298   if (!fFMDLoader->TreeD()->GetEvent(0)) return kFALSE;
 299
 300   // Reads the digits from the array, and fills up the cache (fAdcs)
 301   fAdcs.Clear();
 302   Int_t nDigits = digits->GetEntries();
 303   for (Int_t digit = 0; digit < nDigits; digit++) {
 304     AliFMDDigit* fmdDigit =
 305       static_cast<AliFMDDigit*>(digits->UncheckedAt(digit));
 306
 307     ProcessDigit(fmdDigit);
 308   } //digit loop
 309   return kTRUE;
 310 }
 311
 312 //____________________________________________________________________
 313 Bool_t
 314 AliFMDReconstructor::ReadAdcs(AliRawReader* reader) const
 315 {
 316   // Read the ADC values from a raw data reader.
 317   AliDebug(10, "Reading ADCs from RawReader");
 318   // Reads the digits from a RAW data
 319   fAdcs.Clear();
 320   // reader->Reset();
 321
 322   if (!reader->ReadHeader()) {
 323     Error("ReadAdcs", "Couldn't read header");
 324     return kFALSE;
 325   }
 326
 327   // Use AliAltroRawStream to read the ALTRO format.  No need to
 328   // reinvent the wheel :-)
 329   AliFMDRawStream input(reader);
 330   // Select FMD DDL's
 331   reader->Select(AliFMD::kBaseDDL >> 8);
 332
 333   Int_t    oldDDL      = -1;
 334   Int_t    count       = 0;
 335   UShort_t detector    = 1; // Must be one here
 336   UShort_t oldDetector = 0;
 337   // Loop over data in file
 338   Bool_t   next       = kTRUE;
 339
 340   // local Cache
 341   TArrayI counts(10);
 342   counts.Reset(-1);
 343   Int_t offset = 0;
 344
 345   while (next) {
 346     next = input.Next();
 347
 348
 349     count++;
 350     Int_t ddl = reader->GetDDLID();
 351     if (ddl != oldDDL
 352         || input.IsNewStrip()
 353         || !next) {
 354       // Make a new digit, if we have some data!
 355       if (counts[0] >= 0) {
 356         // Got a new strip.
 357         AliDebug(10, Form("Add a new strip: FMD%d%c[%2d,%3d] "
 358                           "(current: FMD%d%c[%2d,%3d])",
 359                           oldDetector, input.PrevRing(),
 360                           input.PrevSector() , input.PrevStrip(),
 361                           detector , input.Ring(), input.Sector(),
 362                           input.Strip()));
 363         fFMD->AddDigit(oldDetector,
 364                        input.PrevRing(),
 365                        input.PrevSector(),
 366                        input.PrevStrip(),
 367                        counts[0], counts[1], counts[2]);
 368         AliFMDDigit* digit =
 369           static_cast<AliFMDDigit*>(fFMD->Digits()->
 370                                     UncheckedAt(fFMD->GetNdigits()-1));
 371         ProcessDigit(digit);
 372       }
 373
 374       if (!next) {
 375         AliDebug(10, Form("Read %d channels for FMD%d",
 376                           count + 1, detector));
 377         break;
 378       }
 379
 380
 381       // If we got a new DDL, it means we have a new detector.
 382       if (ddl != oldDDL) {
 383         if (detector != 0)
 384           AliDebug(10, Form("Read %d channels for FMD%d",
 385                             count + 1, detector));
 386         // Reset counts, and update the DDL cache
 387         count       = 0;
 388         oldDDL      = ddl;
 389         // Check that we're processing a FMD detector
 390         Int_t detId = reader->GetDetectorID();
 391         if (detId != (AliFMD::kBaseDDL >> 8)) {
 392           Error("ReadAdcs", "Detector ID %d != %d",
 393                 detId, (AliFMD::kBaseDDL >> 8));
 394           break;
 395         }
 396         // Figure out what detector we're deling with
 397         oldDetector = detector;
 398         switch (ddl) {
 399         case 0: detector = 1; break;
 400         case 1: detector = 2; break;
 401         case 2: detector = 3; break;
 402         default:
 403           Error("ReadAdcs", "Unknown DDL 0x%x for FMD", ddl);
 404           return kFALSE;
 405         }
 406         AliDebug(10, Form("Reading ADCs for 0x%x  - That is FMD%d",
 407                           reader->GetEquipmentId(), detector));
 408       }
 409       counts.Reset(-1);
 410       offset = 0;
 411     }
 412
 413     counts[offset] = input.Count();
 414     offset++;
 415
 416     AliDebug(10, Form("ADC of FMD%d%c[%2d,%3d] += %d",
 417                       detector, input.Ring(), input.Sector(),
 418                       input.Strip(), input.Count()));
 419     oldDetector = detector;
 420   }
 421   return kTRUE;
 422 }
 423
 424 //____________________________________________________________________
 425 void
 426 AliFMDReconstructor::ProcessDigit(AliFMDDigit* digit) const
 427 {
 428   // Process a digit.  Derived classes can overload this member
 429   // function to do stuff to the digit.  However, it should write the
 430   // ADC count to the internal cache
 431   //
 432   //    fAdcs(detector - 1, ring, sector, strip) = counts;
 433   //
 434   // In this implementation, we count the number of strips below
 435   // threshold.   This we do to later choose what kind of
 436   // reconstruction algorithm we'd like to use.
 437   //
 438   UShort_t detector = digit->Detector();
 439   Char_t   ring     = digit->Ring();
 440   UShort_t sector   = digit->Sector();
 441   UShort_t strip    = digit->Strip();
 442
 443   UShort_t counts  = SubtractPedestal(digit);
 444
 445   fAdcs(detector - 1, ring, sector, strip) = counts;
 446   if (counts < fThreshold) fEmptyStrips++;
 447   fTotalStrips++;
 448 }
 449
 450 //____________________________________________________________________
 451 UShort_t
 452 AliFMDReconstructor::SubtractPedestal(AliFMDDigit* digit) const
 453 {
 454   // Member function to subtract the pedestal from a digit
 455   // This implementation does nothing, but a derived class could over
 456   // load this to subtract a pedestal that was given in a database or
 457   // something like that.
 458
 459   Int_t counts =
 460     TMath::Max(Int_t(digit->Count1() - fPedestal - 3 * fPedestalWidth), 0);
 461   if (digit->Count2() >= 0)
 462     counts +=
 463       TMath::Max(Int_t(digit->Count2() - fPedestal - 3 * fPedestalWidth), 0);
 464   if (digit->Count3() >= 0)
 465     counts +=
 466       TMath::Max(Int_t(digit->Count3() - fPedestal - 3 * fPedestalWidth), 0);
 467
 468   if (counts < 0) counts = 0;
 469   return  UShort_t(counts);
 470 }
 471
 472 //____________________________________________________________________
 473 void
 474 AliFMDReconstructor::ReconstructFromCache(Float_t zVertex) const
 475 {
 476   // Based on the information in the cache, do the reconstruction.
 477   Int_t nRecon = 0;
 478   // Loop over the detectors
 479   for (Int_t i = 1; i <= 3; i++) {
 480     AliFMDSubDetector* sub = 0;
 481     switch (i) {
 482     case 1: sub = fFMD->GetFMD1(); break;
 483     case 2: sub = fFMD->GetFMD2(); break;
 484     case 3: sub = fFMD->GetFMD3(); break;
 485     }
 486     if (!sub) continue;
 487
 488     // Loop over the rings in the detector
 489     for (Int_t j = 0; j < 2; j++) {
 490       Float_t     rZ = 0;
 491       AliFMDRing* r  = 0;
 492       switch (j) {
 493       case 0: r  = sub->GetInner(); rZ = sub->GetInnerZ(); break;
 494       case 1: r  = sub->GetOuter(); rZ = sub->GetOuterZ(); break;
 495       }
 496       if (!r) continue;
 497
 498       // Calculate low/high theta and eta
 499       // FIXME: Is this right?
 500       Float_t realZ    = zVertex + rZ;
 501       Float_t thetaOut = TMath::ATan2(r->GetHighR(), realZ);
 502       Float_t thetaIn  = TMath::ATan2(r->GetLowR(), realZ);
 503       Float_t etaOut   = - TMath::Log(TMath::Tan(thetaOut / 2));
 504       Float_t etaIn    = - TMath::Log(TMath::Tan(thetaIn / 2));
 505       if (TMath::Abs(etaOut) > TMath::Abs(etaIn)) {
 506         Float_t tmp = etaIn;
 507         etaIn       = etaOut;
 508         etaOut      = tmp;
 509       }
 510
 511       //-------------------------------------------------------------
 512       //
 513       // Here starts poisson method
 514       //
 515       // Calculate eta step per strip, number of eta steps, number of
 516       // phi steps, and check the sign of the eta increment
 517       Float_t stripEta = (Float_t(r->GetNStrips()) / (etaIn - etaOut));
 518       Int_t   nEta     = Int_t(TMath::Abs(etaIn - etaOut) / fDeltaEta);
 519       Int_t   nPhi     = Int_t(360. / fDeltaPhi);
 520       Float_t sign     = TMath::Sign(Float_t(1.), etaIn);
 521
 522       AliDebug(10, Form("FMD%d%c Eta range: %f, %f %d Phi steps",
 523                         sub->GetId(), r->GetId(), etaOut, etaIn, nPhi));
 524
 525       // Loop over relevant phi values
 526       for (Int_t p = 0; p < nPhi; p++) {
 527         Float_t  minPhi    = p * fDeltaPhi;
 528         Float_t  maxPhi    = minPhi + fDeltaPhi;
 529         UShort_t minSector = UShort_t(minPhi / 360) * r->GetNSectors();
 530         UShort_t maxSector = UShort_t(maxPhi / 360) * r->GetNSectors();
 531
 532         AliDebug(10, Form(" Now in phi range %f, %f (sectors %d,%d)",
 533                           minPhi, maxPhi, minSector, maxSector));
 534         // Loop over relevant eta values
 535         for (Int_t e = nEta; e >= 0; --e) {
 536           Float_t  maxEta   = etaIn  - sign * e * fDeltaEta;
 537           Float_t  minEta   = maxEta - sign * fDeltaEta;
 538           if (sign > 0)  minEta = TMath::Max(minEta, etaOut);
 539           else           minEta = TMath::Min(minEta, etaOut);
 540           Float_t  theta1   = 2 * TMath::ATan(TMath::Exp(-minEta));
 541           Float_t  theta2   = 2 * TMath::ATan(TMath::Exp(-maxEta));
 542           Float_t  minR     = TMath::Abs(realZ * TMath::Tan(theta2));
 543           Float_t  maxR     = TMath::Abs(realZ * TMath::Tan(theta1));
 544           UShort_t minStrip = UShort_t((etaIn - maxEta) * stripEta + 0.5);
 545           UShort_t maxStrip = UShort_t((etaIn - minEta) * stripEta + 0.5);
 546
 547           AliDebug(10, Form("  Now in eta range %f, %f (strips %d, %d)\n"
 548                             "    [radii %f, %f, thetas %f, %f, sign %d]",
 549                             minEta, maxEta, minStrip, maxStrip,
 550                             minR, maxR, theta1, theta2, sign));
 551
 552           // Count number of empty strips
 553           Int_t   emptyStrips = 0;
 554           for (Int_t sector = minSector; sector < maxSector; sector++)
 555             for (Int_t strip = minStrip; strip < maxStrip; strip++)
 556               if (fAdcs(sub->GetId() - 1, r->GetId(), sector, strip)
 557                   < fThreshold) emptyStrips++;
 558
 559           // The total number of strips
 560           Float_t nTotal = (maxSector - minSector) * (maxStrip - minStrip);
 561
 562           // Log ratio of empty to total number of strips
 563           AliDebug(10, Form("Lambda= %d / %d = %f",
 564                             emptyStrips, nTotal,
 565                             Float_t(emptyStrips) / nTotal));
 566
 567           Double_t lambda = (emptyStrips > 0 ?
 568                              - TMath::Log(Double_t(emptyStrips) / nTotal) :
 569                              1);
 570
 571           // The reconstructed number of particles is then given by
 572           Int_t reconstructed = Int_t(lambda * nTotal + 0.5);
 573
 574           // Add a AliFMDParticles to the reconstruction tree.
 575           new((*fParticles)[nRecon])
 576             AliFMDParticles(sub->GetId(), r->GetId(),
 577                             minSector, maxSector, minStrip, maxStrip,
 578                             minEta, maxEta, minPhi, maxPhi,
 579                             reconstructed, AliFMDParticles::kPoission);
 580           nRecon++;
 581         } // phi
 582       } // eta
 583     } // ring
 584   } // detector
 585 }
 586
 587
 588 //____________________________________________________________________
 589 void
 590 AliFMDReconstructor::FillESD(AliRunLoader* /*fRunLoader*/,
 591                              AliESD* /*esd*/) const
 592 {
 593   // nothing to be done
 594
 595 }
 596
 597 //____________________________________________________________________
 598 //
 599 // EOF
 600 //