FMD/AliFMDDigitizer.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2004, 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 class contains the procedures simulation ADC  signal for the
  21 //  Forward Multiplicity detector  : Hits->Digits and Hits->SDigits
  22 //
  23 //  Digits consists of
  24 //   - Detector #
  25 //   - Ring ID
  26 //   - Sector #
  27 //   - Strip #
  28 //   - ADC count in this channel
  29 //
  30 //  Digits consists of
  31 //   - Detector #
  32 //   - Ring ID
  33 //   - Sector #
  34 //   - Strip #
  35 //   - Total energy deposited in the strip
  36 //   - ADC count in this channel
  37 //
  38 // As the Digits and SDigits have so much in common, the classes
  39 // AliFMDDigitizer and AliFMDSDigitizer are implemented via a base
  40 // class AliFMDBaseDigitizer.
  41 //
  42 //                 +---------------------+
  43 //                 | AliFMDBaseDigitizer |
  44 //                 +---------------------+
  45 //                           ^
  46 //                           |
  47 //                +----------+---------+
  48 //                |                    |
  49 //      +-----------------+     +------------------+
  50 //      | AliFMDDigitizer |     | AliFMDSDigitizer |
  51 //      +-----------------+     +------------------+
  52 //
  53 // These classes has several paramters:
  54 //
  55 //     fPedestal
  56 //     fPedestalWidth
  57 //         (Only AliFMDDigitizer)
  58 //         Mean and width of the pedestal.  The pedestal is simulated
  59 //         by a Guassian, but derived classes my override MakePedestal
  60 //         to simulate it differently (or pick it up from a database).
  61 //
  62 //     fVA1MipRange
  63 //         The dymamic MIP range of the VA1_ALICE pre-amplifier chip
  64 //
  65 //     fAltroChannelSize
  66 //         The largest number plus one that can be stored in one
  67 //         channel in one time step in the ALTRO ADC chip.
  68 //
  69 //     fSampleRate
  70 //         How many times the ALTRO ADC chip samples the VA1_ALICE
  71 //         pre-amplifier signal.   The VA1_ALICE chip is read-out at
  72 //         10MHz, while it's possible to drive the ALTRO chip at
  73 //         25MHz.  That means, that the ALTRO chip can have time to
  74 //         sample each VA1_ALICE signal up to 2 times.  Although it's
  75 //         not certain this feature will be used in the production,
  76 //         we'd like have the option, and so it should be reflected in
  77 //         the code.
  78 //
  79 //
  80 // The shaping function of the VA1_ALICE is generally given by
  81 //
  82 //      f(x) = A(1 - exp(-Bx))
  83 //
  84 // where A is the total charge collected in the pre-amp., and B is a
  85 // paramter that depends on the shaping time of the VA1_ALICE circut.
  86 //
  87 // When simulating the shaping function of the VA1_ALICe
  88 // pre-amp. chip, we have to take into account, that the shaping
  89 // function depends on the previous value of read from the pre-amp.
  90 //
  91 // That results in the following algorithm:
  92 //
  93 //    last = 0;
  94 //    FOR charge IN pre-amp. charge train DO
  95 //      IF last < charge THEN
  96 //        f(t) = (charge - last) * (1 - exp(-B * t)) + last
  97 //      ELSE
  98 //        f(t) = (last - charge) * exp(-B * t) + charge)
  99 //      ENDIF
 100 //      FOR i IN # samples DO
 101 //        adc_i = f(i / (# samples))
 102 //      DONE
 103 //      last = charge
 104 //   DONE
 105 //
 106 // Here,
 107 //
 108 //   pre-amp. charge train
 109 //       is a series of 128 charges read from the VA1_ALICE chip
 110 //
 111 //   # samples
 112 //       is the number of times the ALTRO ADC samples each of the 128
 113 //       charges from the pre-amp.
 114 //
 115 // Where Q is the total charge collected by the VA1_ALICE
 116 // pre-amplifier.   Q is then given by
 117 //
 118 //           E S
 119 //      Q =  - -
 120 //           e R
 121 //
 122 // where E is the total energy deposited in a silicon strip, R is the
 123 // dynamic range of the VA1_ALICE pre-amp (fVA1MipRange), e is the
 124 // energy deposited by a single MIP, and S ALTRO channel size in each
 125 // time step (fAltroChannelSize).
 126 //
 127 // The energy deposited per MIP is given by
 128 //
 129 //      e = M * rho * w
 130 //
 131 // where M is the universal number 1.664, rho is the density of
 132 // silicon, and w is the depth of the silicon sensor.
 133 //
 134 // The final ADC count is given by
 135 //
 136 //      C' = C + P
 137 //
 138 // where P is the (randomized) pedestal (see MakePedestal)
 139 //
 140 // This class uses the class template AliFMDMap<Type> to make an
 141 // internal cache of the energy deposted of the hits.  The class
 142 // template is instantasized as
 143 //
 144 //  typedef AliFMDMap<std::pair<Float_t, UShort_t> > AliFMDEdepMap;
 145 //
 146 // The first member of the values is the summed energy deposition in a
 147 // given strip, while the second member of the values is the number of
 148 // hits in a given strip.  Using the second member, it's possible to
 149 // do some checks on just how many times a strip got hit, and what
 150 // kind of error we get in our reconstructed hits.  Note, that this
 151 // information is currently not written to the digits tree.  I think a
 152 // QA (Quality Assurance) digit tree is better suited for that task.
 153 // However, the information is there to be used in the future.
 154 //
 155 //
 156 // Latest changes by Christian Holm Christensen
 157 //
 158 //////////////////////////////////////////////////////////////////////////////
 159
 160 //      /1
 161 //      |           A(-1 + B + exp(-B))
 162 //      | f(x) dx = ------------------- = 1
 163 //      |                    B
 164 //      / 0
 165 //
 166 // and B is the a parameter defined by the shaping time (fShapingTime).
 167 //
 168 // Solving the above equation, for A gives
 169 //
 170 //                 B
 171 //      A = ----------------
 172 //          -1 + B + exp(-B)
 173 //
 174 // So, if we define the function g: [0,1] -> [0:1] by
 175 //
 176 //               / v
 177 //               |              Bu + exp(-Bu) - Bv - exp(-Bv)
 178 //      g(u,v) = | f(x) dx = -A -----------------------------
 179 //               |                            B
 180 //               / u
 181 //
 182 // we can evaluate the ALTRO sample of the VA1_ALICE pre-amp between
 183 // any two times (u, v), by
 184 //
 185 //
 186 //                                B         Bu + exp(-Bu) - Bv - exp(-Bv)
 187 //      C = Q g(u,v) = - Q ---------------- -----------------------------
 188 //                         -1 + B + exp(-B)              B
 189 //
 190 //               Bu + exp(-Bu) - Bv - exp(-Bv)
 191 //        = -  Q -----------------------------
 192 //                    -1 + B + exp(-B)
 193 //
 194
 195 #include <TTree.h>              // ROOT_TTree
 196 #include <TRandom.h>            // ROOT_TRandom
 197 #include <AliLog.h>             // ALILOG_H
 198 #include "AliFMDDigitizer.h"    // ALIFMDDIGITIZER_H
 199 #include "AliFMD.h"             // ALIFMD_H
 200 #include "AliFMDGeometry.h"     // ALIFMDGEOMETRY_H
 201 #include "AliFMDDetector.h"     // ALIFMDDETECTOR_H
 202 #include "AliFMDRing.h"         // ALIFMDRING_H
 203 #include "AliFMDHit.h"          // ALIFMDHIT_H
 204 #include "AliFMDDigit.h"        // ALIFMDDIGIT_H
 205 #include "AliFMDParameters.h"   // ALIFMDPARAMETERS_H
 206 #include <AliRunDigitizer.h>    // ALIRUNDIGITIZER_H
 207 #include <AliRun.h>             // ALIRUN_H
 208 #include <AliLoader.h>          // ALILOADER_H
 209 #include <AliRunLoader.h>       // ALIRUNLOADER_H
 210
 211 //====================================================================
 212 ClassImp(AliFMDBaseDigitizer)
 213 #if 0
 214   ; // This is here to keep Emacs for indenting the next line
 215 #endif
 216
 217 //____________________________________________________________________
 218 AliFMDBaseDigitizer::AliFMDBaseDigitizer()
 219   : fRunLoader(0)
 220 {
 221   // Default ctor - don't use it
 222 }
 223
 224 //____________________________________________________________________
 225 AliFMDBaseDigitizer::AliFMDBaseDigitizer(AliRunDigitizer* manager)
 226   : AliDigitizer(manager, "AliFMDBaseDigitizer", "FMD Digitizer base class"),
 227     fRunLoader(0),
 228     fEdep(AliFMDMap::kMaxDetectors,
 229           AliFMDMap::kMaxRings,
 230           AliFMDMap::kMaxSectors,
 231           AliFMDMap::kMaxStrips)
 232 {
 233   // Normal CTOR
 234   AliDebug(1," processed");
 235   SetShapingTime();
 236 }
 237
 238 //____________________________________________________________________
 239 AliFMDBaseDigitizer::AliFMDBaseDigitizer(const Char_t* name,
 240                                          const Char_t* title)
 241   : AliDigitizer(name, title),
 242     fRunLoader(0),
 243     fEdep(AliFMDMap::kMaxDetectors,
 244           AliFMDMap::kMaxRings,
 245           AliFMDMap::kMaxSectors,
 246           AliFMDMap::kMaxStrips)
 247 {
 248   // Normal CTOR
 249   AliDebug(1," processed");
 250   SetShapingTime();
 251 }
 252
 253 //____________________________________________________________________
 254 AliFMDBaseDigitizer::~AliFMDBaseDigitizer()
 255 {
 256   // Destructor
 257 }
 258
 259 //____________________________________________________________________
 260 Bool_t
 261 AliFMDBaseDigitizer::Init()
 262 {
 263   // Initialization
 264   return kTRUE;
 265 }
 266
 267
 268 //____________________________________________________________________
 269 UShort_t
 270 AliFMDBaseDigitizer::MakePedestal(UShort_t,
 271                                   Char_t,
 272                                   UShort_t,
 273                                   UShort_t) const
 274 {
 275   return 0;
 276 }
 277
 278 //____________________________________________________________________
 279 void
 280 AliFMDBaseDigitizer::SumContributions(AliFMD* fmd)
 281 {
 282   // Sum energy deposited contributions from each hit in a cache
 283   // (fEdep).
 284   if (!fRunLoader)
 285     Fatal("SumContributions", "no run loader");
 286
 287   // Clear array of deposited energies
 288   fEdep.Reset();
 289
 290   // Get the FMD loader
 291   AliLoader* inFMD = fRunLoader->GetLoader("FMDLoader");
 292   // And load the hits
 293   inFMD->LoadHits("READ");
 294
 295   // Get the tree of hits
 296   TTree* hitsTree = inFMD->TreeH();
 297   if (!hitsTree)  {
 298     // Try again
 299     inFMD->LoadHits("READ");
 300     hitsTree = inFMD->TreeH();
 301   }
 302
 303   // Get the FMD branch
 304   TBranch* hitsBranch = hitsTree->GetBranch("FMD");
 305   if (hitsBranch) fmd->SetHitsAddressBranch(hitsBranch);
 306   else            AliFatal("Branch FMD hit not found");
 307
 308   // Get a list of hits from the FMD manager
 309   TClonesArray *fmdHits = fmd->Hits();
 310
 311   // Get number of entries in the tree
 312   Int_t ntracks  = Int_t(hitsTree->GetEntries());
 313
 314   AliFMDParameters* param = AliFMDParameters::Instance();
 315   Int_t read = 0;
 316   // Loop over the tracks in the
 317   for (Int_t track = 0; track < ntracks; track++)  {
 318     // Read in entry number `track'
 319     read += hitsBranch->GetEntry(track);
 320
 321     // Get the number of hits
 322     Int_t nhits = fmdHits->GetEntries ();
 323     for (Int_t hit = 0; hit < nhits; hit++) {
 324       // Get the hit number `hit'
 325       AliFMDHit* fmdHit =
 326         static_cast<AliFMDHit*>(fmdHits->UncheckedAt(hit));
 327
 328       // Extract parameters
 329       UShort_t detector = fmdHit->Detector();
 330       Char_t   ring     = fmdHit->Ring();
 331       UShort_t sector   = fmdHit->Sector();
 332       UShort_t strip    = fmdHit->Strip();
 333       Float_t  edep     = fmdHit->Edep();
 334
 335       // Check if strip is `dead'
 336       if (param->IsDead(detector, ring, sector, strip)) continue;
 337
 338       // Give warning in case of double hit
 339       if (fEdep(detector, ring, sector, strip).fEdep != 0)
 340         AliDebug(5, Form("Double hit in %d%c(%d,%d)",
 341                          detector, ring, sector, strip));
 342
 343       // Sum energy deposition
 344       fEdep(detector, ring, sector, strip).fEdep  += edep;
 345       fEdep(detector, ring, sector, strip).fN     += 1;
 346       // Add this to the energy deposited for this strip
 347     }  // hit loop
 348   } // track loop
 349   AliDebug(1, Form("Size of cache: %d bytes, read %d bytes",
 350                    sizeof(fEdep), read));
 351 }
 352
 353 //____________________________________________________________________
 354 void
 355 AliFMDBaseDigitizer::DigitizeHits(AliFMD* fmd) const
 356 {
 357   // For the stored energy contributions in the cache (fEdep), convert
 358   // the energy signal to ADC counts, and store the created digit in
 359   // the digits array (AliFMD::fDigits)
 360   //
 361   AliFMDGeometry* geometry = AliFMDGeometry::Instance();
 362
 363   TArrayI counts(3);
 364   for (UShort_t detector=1; detector <= 3; detector++) {
 365     // Get pointer to subdetector
 366     AliFMDDetector* det = geometry->GetDetector(detector);
 367     if (!det) continue;
 368     for (UShort_t ringi = 0; ringi <= 1; ringi++) {
 369       Char_t ring = ringi == 0 ? 'I' : 'O';
 370       // Get pointer to Ring
 371       AliFMDRing* r = det->GetRing(ring);
 372       if (!r) continue;
 373
 374       // Get number of sectors
 375       UShort_t nSectors = UShort_t(360. / r->GetTheta());
 376       // Loop over the number of sectors
 377       for (UShort_t sector = 0; sector < nSectors; sector++) {
 378         // Get number of strips
 379         UShort_t nStrips = r->GetNStrips();
 380         // Loop over the stips
 381         Float_t last = 0;
 382         for (UShort_t strip = 0; strip < nStrips; strip++) {
 383           // Reset the counter array to the invalid value -1
 384           counts.Reset(-1);
 385           // Reset the last `ADC' value when we've get to the end of a
 386           // VA1_ALICE channel.
 387           if (strip % 128 == 0) last = 0;
 388
 389           Float_t edep = fEdep(detector, ring, sector, strip).fEdep;
 390           ConvertToCount(edep, last, detector, ring, sector, strip, counts);
 391           last = edep;
 392           AddDigit(fmd, detector, ring, sector, strip, edep,
 393                    UShort_t(counts[0]), Short_t(counts[1]),
 394                    Short_t(counts[2]));
 395 #if 0
 396           // This checks if the digit created will give the `right'
 397           // number of particles when reconstructed, using a naiive
 398           // approach.  It's here only as a quality check - nothing
 399           // else.
 400           CheckDigit(digit, fEdep(detector, ring, sector, strip).fN,
 401                      counts);
 402 #endif
 403         } // Strip
 404       } // Sector
 405     } // Ring
 406   } // Detector
 407 }
 408
 409 //____________________________________________________________________
 410 void
 411 AliFMDBaseDigitizer::ConvertToCount(Float_t   edep,
 412                                     Float_t   last,
 413                                     UShort_t  detector,
 414                                     Char_t    ring,
 415                                     UShort_t  sector,
 416                                     UShort_t  strip,
 417                                     TArrayI&  counts) const
 418 {
 419   // Convert the total energy deposited to a (set of) ADC count(s).
 420   //
 421   // This is done by
 422   //
 423   //               Energy_Deposited      ALTRO_Channel_Size
 424   //    ADC = -------------------------- ------------------- + pedestal
 425   //          Energy_Deposition_Of_1_MIP VA1_ALICE_MIP_Range
 426   //
 427   //               Energy_Deposited             fAltroChannelSize
 428   //        = --------------------------------- ----------------- + pedestal
 429   //          1.664 * Si_Thickness * Si_Density   fVA1MipRange
 430   //
 431   //
 432   //        = Energy_Deposited * ConversionFactor + pedestal
 433   //
 434   // However, this is modified by the response function of the
 435   // VA1_ALICE pre-amp. chip in case we are doing oversampling of the
 436   // VA1_ALICE output.
 437   //
 438   // In that case, we get N=fSampleRate values of the ADC, and the
 439   // `EnergyDeposited' is a function of which sample where are
 440   // calculating the ADC for
 441   //
 442   //     ADC_i = f(EnergyDeposited, i/N, Last) * ConversionFactor + pedestal
 443   //
 444   // where Last is the Energy deposited in the previous strip.
 445   //
 446   // Here, f is the shaping function of the VA1_ALICE.   This is given
 447   // by
 448   //
 449   //                    |   (E - l) * (1 - exp(-B * t) + l   if E > l
 450   //       f(E, t, l) = <
 451   //                    |   (l - E) * exp(-B * t) + E        otherwise
 452   //
 453   //
 454   //                  = E + (l - E) * ext(-B * t)
 455   //
 456   AliFMDParameters* param = AliFMDParameters::Instance();
 457   Float_t  convF          = 1/param->GetPulseGain(detector,ring,sector,strip);
 458   UShort_t ped            = MakePedestal(detector,ring,sector,strip);
 459   UInt_t   maxAdc         = param->GetAltroChannelSize();
 460   UShort_t rate           = param->GetSampleRate(AliFMDParameters::kBaseDDL);
 461   UShort_t size           = param->GetAltroChannelSize();
 462
 463   // In case we don't oversample, just return the end value.
 464   if (rate == 1) {
 465     counts[0] = UShort_t(TMath::Min(edep * convF + ped, Float_t(size)));
 466     return;
 467   }
 468
 469   // Create a pedestal
 470   Float_t b = fShapingTime;
 471   for (Ssiz_t i = 0; i < rate;  i++) {
 472     Float_t t = Float_t(i) / rate;
 473     Float_t s = edep + (last - edep) * TMath::Exp(-b * t);
 474     counts[i] = UShort_t(TMath::Min(s * convF + ped, Float_t(maxAdc)));
 475   }
 476 }
 477
 478
 479 //====================================================================
 480 ClassImp(AliFMDDigitizer)
 481
 482 //____________________________________________________________________
 483 AliFMDDigitizer::AliFMDDigitizer()
 484   : AliFMDBaseDigitizer()
 485 {
 486   // Default ctor - don't use it
 487 }
 488
 489 //____________________________________________________________________
 490 AliFMDDigitizer::AliFMDDigitizer(AliRunDigitizer* manager)
 491   : AliFMDBaseDigitizer(manager)
 492 {
 493   // Normal CTOR
 494   AliDebug(1," processed");
 495 }
 496
 497 //____________________________________________________________________
 498 void
 499 AliFMDDigitizer::Exec(Option_t*)
 500 {
 501   // Get the output manager
 502   TString outFolder(fManager->GetOutputFolderName());
 503   AliRunLoader* out =
 504     AliRunLoader::GetRunLoader(outFolder.Data());
 505   // Get the FMD output manager
 506   AliLoader* outFMD = out->GetLoader("FMDLoader");
 507
 508   // Get the input loader
 509   TString inFolder(fManager->GetInputFolderName(0));
 510   fRunLoader =
 511     AliRunLoader::GetRunLoader(inFolder.Data());
 512   if (!fRunLoader) {
 513     AliError("Can not find Run Loader for input stream 0");
 514     return;
 515   }
 516   // Get the AliRun object
 517   if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
 518
 519   // Get the AliFMD object
 520   AliFMD* fmd =
 521     static_cast<AliFMD*>(fRunLoader->GetAliRun()->GetDetector("FMD"));
 522   if (!fmd) {
 523     AliError("Can not get FMD from gAlice");
 524     return;
 525   }
 526
 527   Int_t nFiles= fManager->GetNinputs();
 528   for (Int_t inputFile = 0; inputFile < nFiles; inputFile++) {
 529     AliDebug(1,Form(" Digitizing event number %d",
 530                     fManager->GetOutputEventNr()));
 531     // Get the current loader
 532     fRunLoader =
 533       AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inputFile));
 534     if (!fRunLoader) Fatal("Exec", "no run loader");
 535     // Cache contriutions
 536     SumContributions(fmd);
 537   }
 538   // Digitize the event
 539   DigitizeHits(fmd);
 540
 541   // Load digits from the tree
 542   outFMD->LoadDigits("update");
 543   // Get the tree of digits
 544   TTree* digitTree = outFMD->TreeD();
 545   if (!digitTree) {
 546     outFMD->MakeTree("D");
 547     digitTree = outFMD->TreeD();
 548   }
 549   digitTree->Reset();
 550   // Make a branch in the tree
 551   TClonesArray* digits = fmd->Digits();
 552   fmd->MakeBranchInTree(digitTree, fmd->GetName(), &(digits), 4000, 0);
 553   // TBranch* digitBranch = digitTree->GetBranch(fmd->GetName());
 554   // Fill the tree
 555   Int_t write = 0;
 556   write = digitTree->Fill();
 557   AliDebug(1, Form("Wrote %d bytes to digit tree", write));
 558
 559   // Write the digits to disk
 560   outFMD->WriteDigits("OVERWRITE");
 561   outFMD->UnloadHits();
 562   outFMD->UnloadDigits();
 563
 564   // Reset the digits in the AliFMD object
 565   fmd->ResetDigits();
 566 }
 567
 568
 569 //____________________________________________________________________
 570 UShort_t
 571 AliFMDDigitizer::MakePedestal(UShort_t  detector,
 572                               Char_t    ring,
 573                               UShort_t  sector,
 574                               UShort_t  strip) const
 575 {
 576   // Make a pedestal
 577   AliFMDParameters* param =AliFMDParameters::Instance();
 578   Float_t           mean  =param->GetPedestal(detector,ring,sector,strip);
 579   Float_t           width =param->GetPedestalWidth(detector,ring,sector,strip);
 580   return UShort_t(TMath::Max(gRandom->Gaus(mean, width), 0.));
 581 }
 582
 583 //____________________________________________________________________
 584 void
 585 AliFMDDigitizer::AddDigit(AliFMD*  fmd,
 586                           UShort_t detector,
 587                           Char_t   ring,
 588                           UShort_t sector,
 589                           UShort_t strip,
 590                           Float_t  /* edep */,
 591                           UShort_t count1,
 592                           Short_t  count2,
 593                           Short_t  count3) const
 594 {
 595   // Add a digit
 596   fmd->AddDigitByFields(detector, ring, sector, strip, count1, count2, count3);
 597 }
 598
 599 //____________________________________________________________________
 600 void
 601 AliFMDDigitizer::CheckDigit(AliFMDDigit*    digit,
 602                             UShort_t        nhits,
 603                             const TArrayI&  counts)
 604 {
 605   // Check that digit is consistent
 606   AliFMDParameters* param = AliFMDParameters::Instance();
 607   UShort_t          det   = digit->Detector();
 608   Char_t            ring  = digit->Ring();
 609   UShort_t          sec   = digit->Sector();
 610   UShort_t          str   = digit->Strip();
 611   Float_t           mean  = param->GetPedestal(det,ring,sec,str);
 612   Float_t           width = param->GetPedestalWidth(det,ring,sec,str);
 613   UShort_t          range = param->GetVA1MipRange();
 614   UShort_t          size  = param->GetAltroChannelSize();
 615   Int_t             integral = counts[0];
 616   if (counts[1] >= 0) integral += counts[1];
 617   if (counts[2] >= 0) integral += counts[2];
 618   integral -= Int_t(mean + 2 * width);
 619   if (integral < 0) integral = 0;
 620
 621   Float_t convF = Float_t(range) / size;
 622   Float_t mips  = integral * convF;
 623   if (mips > Float_t(nhits) + .5 || mips < Float_t(nhits) - .5)
 624     Warning("CheckDigit", "Digit -> %4.2f MIPS != %d +/- .5 hits",
 625             mips, nhits);
 626 }
 627
 628 //====================================================================
 629 ClassImp(AliFMDSDigitizer)
 630
 631 //____________________________________________________________________
 632 AliFMDSDigitizer::AliFMDSDigitizer()
 633 {
 634   // Default ctor - don't use it
 635 }
 636
 637 //____________________________________________________________________
 638 AliFMDSDigitizer::AliFMDSDigitizer(const Char_t* headerFile,
 639                                    const Char_t* /* sdigfile */)
 640   : AliFMDBaseDigitizer("FMDSDigitizer", "FMD SDigitizer")
 641 {
 642   // Normal CTOR
 643   AliDebug(1," processed");
 644
 645   fRunLoader = AliRunLoader::GetRunLoader(); // Open(headerFile);
 646   if (!fRunLoader)
 647     Fatal("AliFMDSDigitizer", "cannot open session, header file '%s'",
 648           headerFile);
 649   AliLoader* loader = fRunLoader->GetLoader("FMDLoader");
 650   if (!loader)
 651     Fatal("AliFMDSDigitizer", "cannot find FMD loader in specified event");
 652
 653   // Add task to tasks folder
 654   loader->PostSDigitizer(this);
 655
 656 }
 657
 658 //____________________________________________________________________
 659 AliFMDSDigitizer::~AliFMDSDigitizer()
 660 {
 661   // Destructor
 662   AliLoader* loader = fRunLoader->GetLoader("FMDLoader");
 663   loader->CleanSDigitizer();
 664 }
 665
 666 //____________________________________________________________________
 667 void
 668 AliFMDSDigitizer::Exec(Option_t*)
 669 {
 670   // Get the output manager
 671   if (!fRunLoader) {
 672     Error("Exec", "Run loader is not set");
 673     return;
 674   }
 675   if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
 676   if (!fRunLoader->TreeE())     fRunLoader->LoadHeader();
 677
 678   AliLoader* fmdLoader = fRunLoader->GetLoader("FMDLoader");
 679   if (!fmdLoader) Fatal("Exec", "no FMD loader");
 680
 681   // Get the AliFMD object
 682   AliFMD* fmd =
 683     static_cast<AliFMD*>(fRunLoader->GetAliRun()->GetDetector("FMD"));
 684   if (!fmd) {
 685     AliError("Can not get FMD from gAlice");
 686     return;
 687   }
 688
 689   Int_t nEvents = Int_t(fRunLoader->TreeE()->GetEntries());
 690   for (Int_t event = 0; event < nEvents; event++) {
 691     AliDebug(1,Form(" Digitizing event number %d", event));
 692     // Get the current loader
 693     fRunLoader->GetEvent(event);
 694
 695     if (!fmdLoader->TreeS()) fmdLoader->MakeTree("S");
 696     // Make a branch
 697     fmd->MakeBranch("S");
 698
 699     // Cache contriutions
 700     SumContributions(fmd);
 701
 702     // Digitize the event
 703     DigitizeHits(fmd);
 704
 705     fmdLoader->TreeS()->Reset();
 706     fmdLoader->TreeS()->Fill();
 707     fmdLoader->WriteSDigits("OVERWRITE");
 708   }
 709 }
 710
 711 //____________________________________________________________________
 712 void
 713 AliFMDSDigitizer::AddDigit(AliFMD*  fmd,
 714                            UShort_t detector,
 715                            Char_t   ring,
 716                            UShort_t sector,
 717                            UShort_t strip,
 718                            Float_t  edep,
 719                            UShort_t count1,
 720                            Short_t  count2,
 721                            Short_t  count3) const
 722 {
 723   // Add a summable digit
 724   fmd->AddSDigitByFields(detector, ring, sector, strip, edep,
 725                          count1, count2, count3);
 726 }
 727
 728
 729
 730 //____________________________________________________________________
 731 //
 732 // EOF
 733 //
 734
 735
 736
 737