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