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