FMD/AliFMDReconstructor.cxx

   1 //____________________________________________________________________
   2 //
   3 // This is a class that constructs AliFMDRecPoint objects from of Digits
   4 // This class reads either digits from a TClonesArray or raw data from
   5 // a DDL file (or similar), and stores the read ADC counts in an
   6 // internal cache (fAdcs).   The rec-points are made via the naiive
   7 // method.
   8 //
   9 //-- Authors: Evgeny Karpechev(INR) and Alla Maevsksia
  10 //  Latest changes by Christian Holm Christensen <cholm@nbi.dk>
  11 //
  12 //
  13 //____________________________________________________________________
  14 /**************************************************************************
  15  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  16  *                                                                        *
  17  * Author: The ALICE Off-line Project.                                    *
  18  * Contributors are mentioned in the code where appropriate.              *
  19  *                                                                        *
  20  * Permission to use, copy, modify and distribute this software and its   *
  21  * documentation strictly for non-commercial purposes is hereby granted   *
  22  * without fee, provided that the above copyright notice appears in all   *
  23  * copies and that both the copyright notice and this permission notice   *
  24  * appear in the supporting documentation. The authors make no claims     *
  25  * about the suitability of this software for any purpose. It is          *
  26  * provided "as is" without express or implied warranty.                  *
  27  **************************************************************************/
  28 /* $Id$ */
  29 /**
  30  * @file    AliFMDReconstructor.cxx
  31  * @author  Christian Holm Christensen <cholm@nbi.dk>
  32  * @date    Mon Mar 27 12:47:09 2006
  33  * @brief   FMD reconstruction
  34 */
  35
  36 // #include <AliLog.h>                        // ALILOG_H
  37 // #include <AliRun.h>                        // ALIRUN_H
  38 #include "AliFMDDebug.h"
  39 #include "AliFMDGeometry.h"                // ALIFMDGEOMETRY_H
  40 #include "AliFMDParameters.h"              // ALIFMDPARAMETERS_H
  41 #include "AliFMDAltroMapping.h"            // ALIFMDALTROMAPPING_H
  42 #include "AliFMDDigit.h"                   // ALIFMDDIGIT_H
  43 #include "AliFMDSDigit.h"                  // ALIFMDDIGIT_H
  44 #include "AliFMDReconstructor.h"           // ALIFMDRECONSTRUCTOR_H
  45 #include "AliFMDRecoParam.h"               // ALIFMDRECOPARAM_H
  46 #include "AliFMDRawReader.h"               // ALIFMDRAWREADER_H
  47 #include "AliFMDRecPoint.h"                // ALIFMDMULTNAIIVE_H
  48 #include "AliESDEvent.h"                   // ALIESDEVENT_H
  49 #include "AliESDVertex.h"                  // ALIESDVERTEX_H
  50 #include "AliESDTZERO.h"                   // ALIESDVERTEX_H
  51 #include <AliESDFMD.h>                     // ALIESDFMD_H
  52 #include <TMath.h>
  53 #include <TH1.h>
  54 #include <TH2.h>
  55 #include <TFile.h>
  56 #include <climits>
  57 #include "AliFMDESDRevertexer.h"
  58
  59
  60 class AliRawReader;
  61
  62 //____________________________________________________________________
  63 ClassImp(AliFMDReconstructor)
  64 #if 0
  65   ; // This is here to keep Emacs for indenting the next line
  66 #endif
  67
  68 //____________________________________________________________________
  69 AliFMDReconstructor::AliFMDReconstructor()
  70   : AliReconstructor(),
  71     fMult(0x0),
  72     fNMult(0),
  73     fTreeR(0x0),
  74     fCurrentVertex(0),
  75     fESDObj(0x0),
  76     fNoiseFactor(0),
  77     fAngleCorrect(kTRUE),
  78     fVertexType(kNoVertex),
  79     fESD(0x0),
  80     fDiagnostics(kFALSE),
  81     fDiagStep1(0),
  82     fDiagStep2(0),
  83     fDiagStep3(0),
  84     fDiagStep4(0),
  85     fDiagAll(0),
  86     fBad(0),
  87     fZombie(false)
  88 {
  89   // Make a new FMD reconstructor object - default CTOR.
  90   SetNoiseFactor();
  91   SetAngleCorrect();
  92   if (AliDebugLevel() > 0) fDiagnostics = kTRUE;
  93   for(Int_t det = 1; det<=3; det++) {
  94     fZS[det-1]       = kFALSE;
  95     fZSFactor[det-1] = 0;
  96   }
  97 }
  98
  99 //____________________________________________________________________
 100 AliFMDReconstructor::~AliFMDReconstructor()
 101 {
 102   // Destructor
 103   if (fMult)   fMult->Delete();
 104   if (fMult)   delete fMult;
 105   if (fESDObj) delete fESDObj;
 106 }
 107
 108 //____________________________________________________________________
 109 void
 110 AliFMDReconstructor::Init()
 111 {
 112   // Initialize the reconstructor
 113
 114   // Initialize the geometry
 115   AliFMDGeometry* geom = AliFMDGeometry::Instance();
 116   geom->Init();
 117   geom->InitTransformations();
 118
 119   // Initialize the parameters
 120   AliFMDParameters* param = AliFMDParameters::Instance();
 121   if (param->Init() != 0) {
 122     AliError("Failed to initialize parameters, making zombie");
 123     fZombie = true;
 124   }
 125   else
 126     fZombie = false;
 127
 128   // Current vertex position
 129   fCurrentVertex = 0;
 130   // Create array of reconstructed strip multiplicities
 131   // fMult = new TClonesArray("AliFMDRecPoint", 51200);
 132   // Create ESD output object
 133   fESDObj = new AliESDFMD;
 134
 135   // Check if we need diagnostics histograms
 136   if (!fDiagnostics) return;
 137   AliInfo("Making diagnostics histograms");
 138   if (!fDiagStep1) {
 139     fDiagStep1   = new TH2I("diagStep1", "Read ADC vs. Noise surpressed ADC",
 140                             1024, -.5, 1023.5, 1024, -.5, 1023.5);
 141     fDiagStep1->SetDirectory(0);
 142     fDiagStep1->GetXaxis()->SetTitle("ADC (read)");
 143     fDiagStep1->GetYaxis()->SetTitle(Form("ADC (noise surpressed %4.f)",
 144                                           fNoiseFactor));
 145   }
 146   if (!fDiagStep2) {
 147     fDiagStep2  = new TH2F("diagStep2",  "ADC vs Edep deduced",
 148                            1024, -.5, 1023.5, 100, 0, 2);
 149     fDiagStep2->SetDirectory(0);
 150     fDiagStep2->GetXaxis()->SetTitle("ADC (noise surpressed)");
 151     fDiagStep2->GetYaxis()->SetTitle("#Delta E [GeV]");
 152   }
 153   if (!fDiagStep3) {
 154     fDiagStep3  = new TH2F("diagStep3",  "Edep vs Edep path corrected",
 155                            100, 0., 2., 100, 0., 2.);
 156     fDiagStep3->SetDirectory(0);
 157     fDiagStep3->GetXaxis()->SetTitle("#Delta E [GeV]");
 158     fDiagStep3->GetYaxis()->SetTitle("#Delta E/#Delta x #times #delta x [GeV]");
 159   }
 160   if (!fDiagStep4) {
 161     fDiagStep4  = new TH2F("diagStep4",  "Edep vs Multiplicity deduced",
 162                            100, 0., 2., 100, -.1, 19.9);
 163     fDiagStep4->SetDirectory(0);
 164     fDiagStep4->GetXaxis()->SetTitle("#Delta E/#Delta x #times #delta x [GeV]");
 165     fDiagStep4->GetYaxis()->SetTitle("Multiplicity");
 166     fDiagAll    = new TH2F("diagAll",    "Read ADC vs Multiplicity deduced",
 167                            1024, -.5, 1023.5, 100, -.1, 19.9);
 168   }
 169   if (fDiagAll) {
 170     fDiagAll->SetDirectory(0);
 171     fDiagAll->GetXaxis()->SetTitle("ADC (read)");
 172     fDiagAll->GetYaxis()->SetTitle("Multiplicity");
 173   }
 174 }
 175
 176 //____________________________________________________________________
 177 void
 178 AliFMDReconstructor::ConvertDigits(AliRawReader* reader,
 179                                    TTree* digitsTree) const
 180 {
 181   // Convert Raw digits to AliFMDDigit's in a tree
 182   if (fZombie) {
 183     AliWarning("I'm a zombie - cannot do anything");
 184     return;
 185   }
 186   AliFMDDebug(1, ("Reading raw data into digits tree"));
 187   if (!digitsTree) {
 188     AliError("No digits tree passed");
 189     return;
 190   }
 191   static TClonesArray* array = new TClonesArray("AliFMDDigit");
 192   digitsTree->Branch("FMD", &array);
 193   array->Clear();
 194
 195   AliFMDRawReader rawRead(reader, digitsTree);
 196   // rawRead.SetSampleRate(fFMD->GetSampleRate());
 197   // rawRead.Exec();
 198   rawRead.ReadAdcs(array);
 199
 200   Int_t nWrite = digitsTree->Fill();
 201   AliFMDDebug(1, ("Got a grand total of %d digits, wrote %d bytes to tree",
 202                    array->GetEntriesFast(), nWrite));
 203
 204
 205   AliFMDAltroMapping* map = AliFMDParameters::Instance()->GetAltroMap();
 206   for (size_t i = 1; i <= 3; i++) {
 207     fZS[i-1]       = rawRead.IsZeroSuppressed(map->Detector2DDL(i));
 208     fZSFactor[i-1] = rawRead.NoiseFactor(map->Detector2DDL(i));
 209   }
 210 }
 211
 212 //____________________________________________________________________
 213 void
 214 AliFMDReconstructor::GetVertex(AliESDEvent* esd) const
 215 {
 216   // Return the vertex to use.
 217   // This is obtained from the ESD object.
 218   // If not found, a warning is issued.
 219   fVertexType    = kNoVertex;
 220   fCurrentVertex = 0;
 221   if (!esd) return;
 222
 223   const AliESDVertex* vertex = esd->GetPrimaryVertex();
 224   if (!vertex)        vertex = esd->GetPrimaryVertexSPD();
 225   if (!vertex)        vertex = esd->GetPrimaryVertexTPC();
 226   if (!vertex)        vertex = esd->GetVertex();
 227
 228   if (vertex) {
 229     AliFMDDebug(2, ("Got %s (%s) from ESD: %f",
 230                     vertex->GetName(), vertex->GetTitle(), vertex->GetZv()));
 231     fCurrentVertex = vertex->GetZv();
 232     fVertexType    = kESDVertex;
 233     return;
 234   }
 235   else if (esd->GetESDTZERO()) {
 236     AliFMDDebug(2, ("Got primary vertex from T0: %f", esd->GetT0zVertex()));
 237     fCurrentVertex = esd->GetT0zVertex();
 238     fVertexType    = kESDVertex;
 239     return;
 240   }
 241   AliWarning("Didn't get any vertex from ESD or generator");
 242 }
 243
 244 //____________________________________________________________________
 245 Int_t
 246 AliFMDReconstructor::GetIdentifier() const
 247 {
 248   // Get the detector identifier.
 249   // Note the actual value is cached so that we do not
 250   // need to do many expensive string comparisons.
 251   static Int_t idx = AliReconstruction::GetDetIndex(GetDetectorName());
 252   return idx;
 253 }
 254
 255 //____________________________________________________________________
 256 const AliFMDRecoParam*
 257 AliFMDReconstructor::GetParameters() const
 258 {
 259   // Get the reconstruction parameters.
 260   //
 261   // Return:
 262   //   Pointer to reconstruction parameters or null if not found or wrong type
 263   Int_t iDet = GetIdentifier(); // Was 12 - but changed on Cvetans request
 264   const AliDetectorRecoParam* params = AliReconstructor::GetRecoParam(iDet);
 265   if (!params || params->IsA() != AliFMDRecoParam::Class()) return 0;
 266   return static_cast<const AliFMDRecoParam*>(params);
 267 }
 268
 269 //____________________________________________________________________
 270 void
 271 AliFMDReconstructor::UseRecoParam(Bool_t set) const
 272 {
 273   //
 274   // Set-up reconstructor to use values from reconstruction
 275   // parameters, if present, for this event.   If the argument @a set
 276   // is @c false, then restore preset values.
 277   //
 278   // Parameters:
 279   //    set
 280   //
 281   static Float_t savedNoiseFactor  = fNoiseFactor;
 282   static Bool_t  savedAngleCorrect = fAngleCorrect;
 283   if (set) {
 284     const AliFMDRecoParam* params  = GetParameters();
 285     if (params) {
 286       fNoiseFactor  = params->NoiseFactor();
 287       fAngleCorrect = params->AngleCorrect();
 288     }
 289     return;
 290   }
 291   fNoiseFactor  = savedNoiseFactor;
 292   fAngleCorrect = savedAngleCorrect;
 293 }
 294
 295
 296 //____________________________________________________________________
 297 void
 298 AliFMDReconstructor::MarkDeadChannels(AliESDFMD* esd) const
 299 {
 300   // Loop over all entries of the ESD and mark
 301   // those that are dead as such
 302   // - otherwise put in the zero signal.
 303   AliFMDParameters* param = AliFMDParameters::Instance();
 304
 305   for (UShort_t d = 1; d <= 3; d++) {
 306     UShort_t nR = (d == 1 ? 1 : 2);
 307
 308     for (UShort_t q = 0; q < nR; q++) {
 309       Char_t   r  = (q == 0 ? 'I' : 'O');
 310       UShort_t nS = (q == 0 ?  20 :  40);
 311       UShort_t nT = (q == 0 ? 512 : 256);
 312
 313       for (UShort_t s = 0; s < nS; s++) {
 314         for (UShort_t t = 0; t < nT; t++) {
 315           // A strip can be marked `bad' for two reasons:
 316           //
 317           // - The raw reader fails to read the value
 318           // - The strip is marked in OCDB as bad (IsDead)
 319           //
 320           // Hence, a dead strip will always be marked kInvalid here,
 321           // while a non-dead bad-read strip will be filled with 0.
 322           if (fBad(d, r, s, t)) {
 323             AliDebug(5, Form("Marking FMD%d%c[%2d,%3d] as bad", d, r, s, t));
 324             esd->SetMultiplicity(d, r, s, t, AliESDFMD::kInvalidMult);
 325           }
 326           if (param->IsDead(d, r, s, t)) {
 327             AliDebug(5, Form("Marking FMD%d%c[%2d,%3d] as dead", d, r, s, t));
 328             esd->SetMultiplicity(d, r, s, t, AliESDFMD::kInvalidMult);
 329             // esd->SetEta(d, r, s, t, AliESDFMD::kInvalidEta);
 330           }
 331           else if (esd->Multiplicity(d, r, s, t) == AliESDFMD::kInvalidMult) {
 332             AliDebug(10, Form("Setting null signal in FMD%d%c[%2d,%3d]",
 333                               d, r, s, t));
 334             esd->SetMultiplicity(d, r, s, t, 0);
 335           }
 336         }
 337       }
 338     }
 339   }
 340 }
 341
 342 //____________________________________________________________________
 343 Bool_t
 344 AliFMDReconstructor::PreReconstruct() const
 345 {
 346   AliFMDDebug(1, ("Before reoconstructing"));
 347   if (fZombie) {
 348     AliWarning("I'm a zombie - cannot do anything");
 349     return false;
 350   }
 351
 352   // Get our vertex
 353   GetVertex(fESD);
 354
 355   // Reset bad flags
 356   fBad.Reset(false);
 357
 358   // Reset the output ESD
 359   if (fESDObj) {
 360     fESDObj->Clear();
 361
 362     // Pre-set eta values
 363     for (UShort_t d=1; d<=3; d++) {
 364       UShort_t nQ = (d == 1 ? 1 : 2);
 365       for (UShort_t q=0; q<nQ; q++) {
 366         UShort_t nStr = (q == 0 ? 512 : 256);
 367         Char_t   r    = (q == 0 ? 'I' : 'O');
 368
 369         for (UShort_t t = 0; t < nStr; t++) {
 370           Float_t eta, phi;
 371           // Always use sector 0
 372           PhysicalCoordinates(d, r, 0, t, eta, phi);
 373           fESDObj->SetEta(d, r, 0, t, eta);
 374         }
 375       }
 376     }
 377   }
 378
 379
 380   return true;
 381 }
 382
 383 //____________________________________________________________________
 384 void
 385 AliFMDReconstructor::Reconstruct(AliFMDRawReader& rawReader) const
 386 {
 387   // Reconstruct directly from raw data (no intermediate output on
 388   // digit tree or rec point tree).
 389   //
 390   // Parameters:
 391   //   rawReader        FMD Raw event reader
 392   AliFMDDebug(1, ("Reconstructing from FMD raw reader"));
 393   if (!PreReconstruct()) return;
 394
 395   UShort_t det, sec, str, fac;
 396   Short_t  adc, oldDet = -1;
 397   Bool_t   zs;
 398   Char_t   rng;
 399
 400   UseRecoParam(kTRUE);
 401   while (rawReader.NextSignal(det, rng, sec, str, adc, zs, fac)) {
 402     if (det != oldDet) {
 403       fZS[det-1]       = zs;
 404       fZSFactor[det-1] = fac;
 405       oldDet           = det;
 406     }
 407     ProcessSignal(det, rng, sec, str, adc);
 408   }
 409   UseRecoParam(kFALSE);
 410
 411 }
 412
 413 //____________________________________________________________________
 414 void
 415 AliFMDReconstructor::Reconstruct(AliRawReader* reader, TTree*) const
 416 {
 417   // Reconstruct directly from raw data (no intermediate output on
 418   // digit tree or rec point tree).
 419   //
 420   // Parameters:
 421   //   reader   Raw event reader
 422   //   ctree    Not used - 'cluster tree' to store rec-points in.
 423   AliFMDDebug(1, ("Reconstructing from raw reader"));
 424   if (fZombie) {
 425     AliWarning("I'm a zombie - cannot do anything");
 426     return;
 427   }
 428   AliFMDRawReader rawReader(reader, 0);
 429   Reconstruct(rawReader);
 430 }
 431
 432 //____________________________________________________________________
 433 void
 434 AliFMDReconstructor::Digitize(AliRawReader* reader, TClonesArray* sdigits) const
 435 {
 436   // Reconstruct directly from raw data (no intermediate output on
 437   // digit tree or rec point tree).
 438   //
 439   // Parameters:
 440   //   reader   Raw event reader
 441   //   ctree    Not used.
 442   if (fZombie) {
 443     AliWarning("I'm a zombie - cannot do anything");
 444     return;
 445   }
 446   AliFMDRawReader rawReader(reader, 0);
 447
 448   UShort_t det, sec, str, sam, rat, fac;
 449   Short_t  adc, oldDet = -1;
 450   Bool_t   zs;
 451   Char_t   rng;
 452
 453   UseRecoParam(kTRUE);
 454   while (rawReader.NextSample(det, rng, sec, str, sam, rat, adc, zs, fac)) {
 455     if (!rawReader.SelectSample(sam, rat)) continue;
 456     if (det != oldDet) {
 457       fZS[det-1]       = zs;
 458       fZSFactor[det-1] = fac;
 459       oldDet           = det;
 460     }
 461     DigitizeSignal(sdigits, det, rng, sec, str, sam, adc);
 462   }
 463   UseRecoParam(kFALSE);
 464 }
 465
 466 //____________________________________________________________________
 467 void
 468 AliFMDReconstructor::Reconstruct(TTree* digitsTree,
 469                                  TTree* clusterTree) const
 470 {
 471   // Reconstruct event from digits in tree
 472   // Get the FMD branch holding the digits.
 473   // FIXME: The vertex may not be known yet, so we may have to move
 474   // some of this to FillESD.
 475   //
 476   // Parameters:
 477   //   digitsTree       Pointer to a tree containing digits
 478   //   clusterTree      Pointer to output tree
 479   //
 480   if (!PreReconstruct()) return;
 481
 482   if (!fMult) fMult = new TClonesArray("AliFMDRecPoint");
 483
 484   AliFMDDebug(1, ("Reconstructing from digits in a tree"));
 485
 486   // Get the digitis array
 487   static TClonesArray* digits      = new TClonesArray("AliFMDDigit");
 488   TBranch*             digitBranch = digitsTree->GetBranch("FMD");
 489   if (!digitBranch) {
 490     Error("Exec", "No digit branch for the FMD found");
 491     return;
 492   }
 493   digitBranch->SetAddress(&digits);
 494
 495   if (digits)  digits->Clear();
 496   if (fMult)   fMult->Clear();
 497
 498   // Create rec-point output branch
 499   fNMult = 0;
 500   fTreeR = clusterTree;
 501   fTreeR->Branch("FMD", &fMult);
 502
 503   AliFMDDebug(5, ("Getting entry 0 from digit branch"));
 504   digitBranch->GetEntry(0);
 505
 506   AliFMDDebug(5, ("Processing digits"));
 507   UseRecoParam(kTRUE);
 508   ProcessDigits(digits);
 509   UseRecoParam(kFALSE);
 510
 511   Int_t written = clusterTree->Fill();
 512   AliFMDDebug(10, ("Filled %d bytes into cluster tree", written));
 513   // digits->Delete();
 514   // delete digits;
 515 }
 516
 517
 518 //____________________________________________________________________
 519 void
 520 AliFMDReconstructor::ProcessDigits(TClonesArray* digits,
 521                                    const AliFMDRawReader& rawRead) const
 522 {
 523   // For each digit, find the pseudo rapdity, azimuthal angle, and
 524   // number of corrected ADC counts, and pass it on to the algorithms
 525   // used.
 526   //
 527   // Parameters:
 528   //    digits  Array of digits
 529   //
 530   if (fZombie) {
 531     AliWarning("I'm a zombie - cannot do anything");
 532     return;
 533   }
 534   AliFMDAltroMapping* map = AliFMDParameters::Instance()->GetAltroMap();
 535   for (size_t i = 1; i <= 3; i++) {
 536     fZS[i-1]       = rawRead.IsZeroSuppressed(map->Detector2DDL(i));
 537     fZSFactor[i-1] = rawRead.NoiseFactor(map->Detector2DDL(i));
 538   }
 539   UseRecoParam(kTRUE);
 540   ProcessDigits(digits);
 541   UseRecoParam(kFALSE);
 542 }
 543
 544 //____________________________________________________________________
 545 void
 546 AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
 547 {
 548   // For each digit, find the pseudo rapdity, azimuthal angle, and
 549   // number of corrected ADC counts, and pass it on to the algorithms
 550   // used.
 551   //
 552   // Parameters:
 553   //    digits  Array of digits
 554   //
 555   Int_t nDigits = digits->GetEntries();
 556   AliFMDDebug(2, ("Got %d digits", nDigits));
 557   fESDObj->SetNoiseFactor(fNoiseFactor);
 558   fESDObj->SetAngleCorrected(fAngleCorrect);
 559   fBad.Reset(false);
 560   for (Int_t i = 0; i < nDigits; i++) {
 561     AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
 562     if (!digit) continue;
 563     ProcessDigit(digit);
 564   }
 565 }
 566
 567 //____________________________________________________________________
 568 void
 569 AliFMDReconstructor::ProcessDigit(AliFMDDigit* digit) const
 570 {
 571   //
 572   // Process a single digit
 573   //
 574   // Parameters:
 575   //    digit Digiti to process
 576   //
 577   UShort_t det = digit->Detector();
 578   Char_t   rng = digit->Ring();
 579   UShort_t sec = digit->Sector();
 580   UShort_t str = digit->Strip();
 581   Short_t  adc = digit->Counts();
 582
 583   ProcessSignal(det, rng, sec, str, adc);
 584 }
 585
 586 //____________________________________________________________________
 587 void
 588 AliFMDReconstructor::ProcessSignal(UShort_t det,
 589                                    Char_t   rng,
 590                                    UShort_t sec,
 591                                    UShort_t str,
 592                                    Short_t  adc) const
 593 {
 594   // Process the signal from a single strip
 595   //
 596   // Parameters:
 597   //    det     Detector ID
 598   //    rng     Ring ID
 599   //    sec     Sector ID
 600   //    rng     Strip ID
 601   //    adc     ADC counts
 602   //
 603   if (adc >= AliFMDRawReader::kBadSignal) {
 604     AliFMDDebug(3, ("FMD%d%c[%2d,%3d] is marked bad", det, rng, sec, str));
 605     fBad(det,rng,sec,str) = true;
 606     return;
 607   }
 608
 609   // Check that the strip is not marked as dead
 610   AliFMDParameters* param  = AliFMDParameters::Instance();
 611   if (param->IsDead(det, rng, sec, str)) {
 612     AliFMDDebug(3, ("FMD%d%c[%2d,%3d] is dead", det, rng, sec, str));
 613     fBad(det,rng,sec,str) = true;
 614     return;
 615   }
 616
 617   // digit->Print();
 618   // Get eta and phi
 619   Float_t eta = fESDObj->Eta(det, rng, 0, str);
 620
 621   // Substract pedestal.
 622   UShort_t counts   = SubtractPedestal(det, rng, sec, str, adc);
 623   if(counts == USHRT_MAX) return;
 624
 625     // Gain match digits.
 626   Double_t edep     = Adc2Energy(det, rng, sec, str, eta, counts);
 627   // Get rid of nonsense energy
 628   if(edep < 0)  return;
 629
 630   // Make rough multiplicity
 631   Double_t mult     = Energy2Multiplicity(det, rng, sec, str, edep);
 632   // Get rid of nonsense mult
 633   //if (mult > 20) {
 634   //  AliWarning(Form("The mutliplicity in FMD%d%c[%2d,%3d]=%f > 20 "
 635   //                "(ADC: %d, Energy: %f)", det, rng, sec, str, mult,
 636   //                counts, edep));
 637   // }
 638   if (mult < 0)  return;
 639   AliFMDDebug(10, ("FMD%d%c[%2d,%3d]: "
 640                     "ADC: %d, Counts: %d, Energy: %f, Mult: %f",
 641                   det, rng, sec, str, adc, counts, edep, mult));
 642
 643   // Create a `RecPoint' on the output branch.
 644   if (fMult) {
 645     Float_t phi;
 646     PhysicalCoordinates(det, rng, sec, str, eta, phi);
 647
 648     AliFMDRecPoint* m =
 649       new ((*fMult)[fNMult]) AliFMDRecPoint(det, rng, sec, str,
 650                                             eta, phi, edep, mult);
 651     (void)m; // Suppress warnings about unused variables.
 652     fNMult++;
 653   }
 654
 655   fESDObj->SetMultiplicity(det, rng, sec, str, mult);
 656   // fESDObj->SetEta(det, rng, sec, str, eta);
 657
 658   if (fDiagAll) fDiagAll->Fill(adc, mult);
 659
 660 }
 661
 662 //____________________________________________________________________
 663 void
 664 AliFMDReconstructor::DigitizeSignal(TClonesArray* sdigits,
 665                                     UShort_t det,
 666                                     Char_t   rng,
 667                                     UShort_t sec,
 668                                     UShort_t str,
 669                                     UShort_t /* sam */,
 670                                     Short_t  adc) const
 671 {
 672   // Process the signal from a single strip
 673   //
 674   // Parameters:
 675   //    det     Detector ID
 676   //    rng     Ring ID
 677   //    sec     Sector ID
 678   //    rng     Strip ID
 679   //    adc     ADC counts
 680   //
 681   AliFMDParameters* param  = AliFMDParameters::Instance();
 682   // Check that the strip is not marked as dead
 683   if (param->IsDead(det, rng, sec, str)) {
 684     AliFMDDebug(10, ("FMD%d%c[%2d,%3d] is dead", det, rng, sec, str));
 685     return;
 686   }
 687
 688   // Substract pedestal.
 689   UShort_t counts   = SubtractPedestal(det, rng, sec, str, adc);
 690   if(counts == USHRT_MAX || counts == 0) return;
 691
 692     // Gain match digits.
 693   Double_t edep     = Adc2Energy(det, rng, sec, str, counts);
 694   // Get rid of nonsense energy
 695   if(edep < 0)  return;
 696
 697   Int_t n = sdigits->GetEntriesFast();
 698   // AliFMDSDigit* sdigit =
 699   new ((*sdigits)[n])
 700     AliFMDSDigit(det, rng, sec, str, edep, counts, counts, counts, counts);
 701   // sdigit->SetCount(sam, counts);
 702 }
 703
 704 //____________________________________________________________________
 705 UShort_t
 706 AliFMDReconstructor::SubtractPedestal(UShort_t det,
 707                                       Char_t   rng,
 708                                       UShort_t sec,
 709                                       UShort_t str,
 710                                       UShort_t adc,
 711                                       Float_t  noiseFactor,
 712                                       Bool_t   zsEnabled,
 713                                       UShort_t zsNoiseFactor) const
 714 {
 715   //
 716   // Subtract the pedestal off the ADC counts.
 717   //
 718   // Parameters:
 719   //    det           Detector number
 720   //    rng           Ring identifier
 721   //    sec           Sector number
 722   //    str           Strip number
 723   //    adc           ADC counts
 724   //    noiseFactor   If pedestal substracted pedestal is less then
 725   //        this times the noise, then consider this to be 0.
 726   //    zsEnabled     Whether zero-suppression is on.
 727   //    zsNoiseFactor Noise factor used in on-line pedestal
 728   //        subtraction.
 729   //
 730   // Return:
 731   //    The pedestal subtracted ADC counts (possibly 0), or @c
 732   //         USHRT_MAX in case of problems.
 733   //
 734   AliFMDParameters* param  = AliFMDParameters::Instance();
 735   Float_t           ped    = (zsEnabled ? 0 :
 736                               param->GetPedestal(det, rng, sec, str));
 737   Float_t           noise  = param->GetPedestalWidth(det, rng, sec, str);
 738   if(ped < 0 || noise < 0) {
 739     AliWarningClass(Form("Invalid pedestal (%f) or noise (%f) "
 740                          "for FMD%d%c[%02d,%03d]",
 741                     ped, noise, det, rng, sec, str));
 742     return USHRT_MAX;
 743   }
 744   AliDebugClass(15, Form("Subtracting pedestal for FMD%d%c[%2d,%3d]=%4d "
 745                          "(%s w/factor %d, noise factor %f, "
 746                          "pedestal %8.2f+/-%8.2f)",
 747                          det, rng, sec, str, adc,
 748                          (zsEnabled ? "zs'ed" : "straight"),
 749                          zsNoiseFactor, noiseFactor, ped, noise));
 750
 751   Int_t counts = adc + Int_t(zsEnabled ? zsNoiseFactor * noise : - ped);
 752   counts =  TMath::Max(Int_t(counts), 0);
 753   // Calculate the noise factor for suppressing remenants of the noise
 754   // peak.  If we have done on-line zero suppression, we only check
 755   // for noise signals that are larger than the suppressed noise.  If
 756   // the noise factor used on line is larger than the factor used
 757   // here, we do not do this check at all.
 758   //
 759   // For example:
 760   //    Online factor  |  Read factor |  Result
 761   //    ---------------+--------------+-------------------------------
 762   //           2       |      3       | Check if signal > 1 * noise
 763   //           3       |      3       | Check if signal > 0
 764   //           3       |      2       | Check if signal > 0
 765   //
 766   // In this way, we make sure that we do not suppress away too much
 767   // data, and that the read-factor is the most stringent cut.
 768   Float_t nf = TMath::Max(0.F, noiseFactor - (zsEnabled ? zsNoiseFactor : 0));
 769   if (counts < noise * nf) counts = 0;
 770   if (counts > 0) AliDebugClass(15, "Got a hit strip");
 771
 772   UShort_t ret = counts < 0 ? 0 : counts;
 773   return ret;
 774 }
 775
 776
 777 //____________________________________________________________________
 778 UShort_t
 779 AliFMDReconstructor::SubtractPedestal(UShort_t det,
 780                                       Char_t   rng,
 781                                       UShort_t sec,
 782                                       UShort_t str,
 783                                       Short_t  adc) const
 784 {
 785   // Member function to subtract the pedestal from a digit
 786   //
 787   // Parameters:
 788   //    det     Detector ID
 789   //    rng     Ring ID
 790   //    sec     Sector ID
 791   //    rng     Strip ID
 792   //    adc     # of ADC counts
 793   // Return:
 794   //    Pedestal subtracted signal or USHRT_MAX in case of problems
 795   //
 796   UShort_t counts = SubtractPedestal(det, rng, sec, str, adc,
 797                                      fNoiseFactor, fZS[det-1],
 798                                      fZSFactor[det-1]);
 799   if (fDiagStep1) fDiagStep1->Fill(adc, counts);
 800
 801   return counts;
 802 }
 803
 804 //____________________________________________________________________
 805 Float_t
 806 AliFMDReconstructor::Adc2Energy(UShort_t det,
 807                                 Char_t   rng,
 808                                 UShort_t sec,
 809                                 UShort_t str,
 810                                 UShort_t count) const
 811 {
 812   // Converts number of ADC counts to energy deposited.
 813   // Note, that this member function can be overloaded by derived
 814   // classes to do strip-specific look-ups in databases or the like,
 815   // to find the proper gain for a strip.
 816   //
 817   // In the first simple version, we calculate the energy deposited as
 818   //
 819   //    EnergyDeposited = cos(theta) * gain * count
 820   //
 821   // where
 822   //
 823   //           Pre_amp_MIP_Range
 824   //    gain = ----------------- * Energy_deposited_per_MIP
 825   //           ADC_channel_size
 826   //
 827   // is constant and the same for all strips.
 828   //
 829   // For the production we use the conversion measured in the NBI lab.
 830   // The total conversion is then:
 831   //
 832   //    gain = ADC / DAC
 833   //
 834   //                  EdepMip * count
 835   //      => energy = ----------------
 836   //                  gain * DACPerADC
 837   //
 838   // Parameters:
 839   //    det     Detector ID
 840   //    rng     Ring ID
 841   //    sec     Sector ID
 842   //    rng     Strip ID
 843   //    counts  Number of ADC counts over pedestal
 844   // Return
 845   //    The energy deposited in a single strip, or -1 in case of problems
 846   //
 847   if (count <= 0) return 0;
 848   AliFMDParameters* param = AliFMDParameters::Instance();
 849   Float_t           gain  = param->GetPulseGain(det, rng, sec, str);
 850   // 'Tagging' bad gains as bad energy
 851   if (gain < 0) {
 852     AliWarning(Form("Invalid gain (%f) for FMD%d%c[%02d,%03d]",
 853                     gain, det, rng, sec, str));
 854     return -1;
 855   }
 856   AliFMDDebug(5, ("Converting counts %d to energy (factor=%f, DAC2MIP=%f)",
 857                   count, gain,param->GetDACPerMIP()));
 858
 859   Double_t edep  = ((count * param->GetEdepMip())
 860                     / (gain * param->GetDACPerMIP()));
 861   return edep;
 862 }
 863
 864 //____________________________________________________________________
 865 Float_t
 866 AliFMDReconstructor::Adc2Energy(UShort_t det,
 867                                 Char_t   rng,
 868                                 UShort_t sec,
 869                                 UShort_t str,
 870                                 Float_t  eta,
 871                                 UShort_t count) const
 872 {
 873   // Converts number of ADC counts to energy deposited.
 874   // Note, that this member function can be overloaded by derived
 875   // classes to do strip-specific look-ups in databases or the like,
 876   // to find the proper gain for a strip.
 877   //
 878   // In the first simple version, we calculate the energy deposited as
 879   //
 880   //    EnergyDeposited = cos(theta) * gain * count
 881   //
 882   // where
 883   //
 884   //           Pre_amp_MIP_Range
 885   //    gain = ----------------- * Energy_deposited_per_MIP
 886   //           ADC_channel_size
 887   //
 888   // is constant and the same for all strips.
 889   //
 890   // For the production we use the conversion measured in the NBI lab.
 891   // The total conversion is then:
 892   //
 893   //    gain = ADC / DAC
 894   //
 895   //                  EdepMip * count
 896   //      => energy = ----------------
 897   //                  gain * DACPerADC
 898   //
 899   // Parameters:
 900   //    det     Detector ID
 901   //    rng     Ring ID
 902   //    sec     Sector ID
 903   //    rng     Strip ID
 904   //    eta     Psuedo-rapidity
 905   //    counts  Number of ADC counts over pedestal
 906   // Return
 907   //    The energy deposited in a single strip, or -1 in case of problems
 908   //
 909   Double_t edep = Adc2Energy(det, rng, sec, str, count);
 910
 911   if (fDiagStep2) fDiagStep2->Fill(count, edep);
 912   if (fAngleCorrect) {
 913     Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
 914     Double_t corr  = TMath::Abs(TMath::Cos(theta));
 915     Double_t cedep = corr * edep;
 916     AliFMDDebug(10, ("correcting for path %f * %f = %f (eta=%f, theta=%f)",
 917                       edep, corr, cedep, eta, theta));
 918     if (fDiagStep3) fDiagStep3->Fill(edep, cedep);
 919     edep = cedep;
 920   }
 921   return edep;
 922 }
 923
 924 //____________________________________________________________________
 925 Float_t
 926 AliFMDReconstructor::Energy2Multiplicity(UShort_t /*det*/,
 927                                          Char_t   /*rng*/,
 928                                          UShort_t /*sec*/,
 929                                          UShort_t /*str*/,
 930                                          Float_t  edep) const
 931 {
 932   // Converts an energy signal to number of particles.
 933   // Note, that this member function can be overloaded by derived
 934   // classes to do strip-specific look-ups in databases or the like,
 935   // to find the proper gain for a strip.
 936   //
 937   // In this simple version, we calculate the multiplicity as
 938   //
 939   //   multiplicity = Energy_deposited / Energy_deposited_per_MIP
 940   //
 941   // where
 942   //
 943   //   Energy_deposited_per_MIP = 1.664 * SI_density * SI_thickness
 944   //
 945   // is constant and the same for all strips
 946   //
 947   // Parameters:
 948   //    det     Detector ID
 949   //    rng     Ring ID
 950   //    sec     Sector ID
 951   //    rng     Strip ID
 952   //    edep    Energy deposited in a single strip
 953   // Return
 954   //    The "bare" multiplicity corresponding to the energy deposited
 955   AliFMDParameters* param   = AliFMDParameters::Instance();
 956   Double_t          edepMIP = param->GetEdepMip();
 957   Float_t           mult    = edep / edepMIP;
 958 #if 0
 959   if (edep > 0)
 960     AliFMDDebug(15, ("Translating energy %f to multiplicity via "
 961                      "divider %f->%f", edep, edepMIP, mult));
 962 #endif
 963   if (fDiagStep4) fDiagStep4->Fill(edep, mult);
 964   return mult;
 965 }
 966
 967 //____________________________________________________________________
 968 void
 969 AliFMDReconstructor::PhysicalCoordinates(UShort_t det,
 970                                          Char_t   rng,
 971                                          UShort_t sec,
 972                                          UShort_t str,
 973                                          Float_t& eta,
 974                                          Float_t& phi) const
 975 {
 976   // Get the eta and phi of a digit
 977   //
 978   // Parameters:
 979   //    det     Detector ID
 980   //    rng     Ring ID
 981   //    sec     Sector ID
 982   //    rng     Strip ID
 983   //    eta     On return, contains the psuedo-rapidity of the strip
 984   //    phi     On return, contains the azimuthal angle of the strip
 985   //
 986   AliFMDGeometry* geom = AliFMDGeometry::Instance();
 987   Double_t x, y, z, r, theta, deta, dphi;
 988   geom->Detector2XYZ(det, rng, sec, str, x, y, z);
 989
 990   // Correct for vertex offset.
 991   z     -= fCurrentVertex;
 992   AliFMDGeometry::XYZ2REtaPhiTheta(x, y, z, r, deta, dphi, theta);
 993   eta = deta;
 994   phi = dphi;
 995 }
 996
 997 namespace {
 998   class ESDPrinter : public AliESDFMD::ForOne
 999   {
1000   public:
1001     ESDPrinter() {}
1002     Bool_t operator()(UShort_t d, Char_t r, UShort_t s, UShort_t t,
1003                       Float_t m, Float_t e)
1004     {
1005       if (m > 0 && m != AliESDFMD::kInvalidMult)
1006         printf("  FMD%d%c[%2d,%3d] = %6.3f / %6.3f\n", d, r, s, t, m, e);
1007       return kTRUE;
1008     }
1009   };
1010 }
1011
1012
1013 //____________________________________________________________________
1014 void
1015 AliFMDReconstructor::FillESD(TTree*  /* digitsTree */,
1016                              TTree*  /* clusterTree */,
1017                              AliESDEvent* esd) const
1018 {
1019   // nothing to be done
1020   // FIXME: The vertex may not be known when Reconstruct is executed,
1021   // so we may have to move some of that member function here.
1022   AliFMDDebug(2, ("Calling FillESD with two trees and one ESD"));
1023   if (fZombie) {
1024     AliWarning("I'm a zombie - cannot do anything");
1025     return;
1026   }
1027   // fESDObj->Print();
1028
1029   // Fix up ESD so that only truely dead channels get the kInvalidMult flag.
1030   MarkDeadChannels(fESDObj);
1031
1032   Double_t oldVz = fCurrentVertex;
1033   GetVertex(esd);
1034   if (fVertexType != kNoVertex) {
1035     AliFMDDebug(2, ("Revertexing the ESD data to vz=%f (was %f)",
1036                     fCurrentVertex, oldVz));
1037     AliFMDESDRevertexer revertexer;
1038     revertexer.Revertex(fESDObj, fCurrentVertex);
1039   }
1040
1041   if (AliDebugLevel() > 10) {
1042     ESDPrinter p;
1043     fESDObj->ForEach(p);
1044   }
1045
1046   if (esd) {
1047     AliFMDDebug(2, ("Writing FMD data to ESD tree"));
1048     esd->SetFMDData(fESDObj);
1049   }
1050
1051   if (!fDiagnostics || !esd) return;
1052   static bool first = true;
1053   // This is most likely NOT the event number you'd like to use. It
1054   // has nothing to do with the 'real' event number.
1055   // - That's OK.  We just use it for the name of the directory -
1056   // nothing else.  Christian
1057   Int_t evno = esd->GetEventNumberInFile();
1058   AliFMDDebug(3, ("Writing diagnostics histograms to FMD.Diag.root/%03d",evno));
1059   TFile f("FMD.Diag.root", (first ? "RECREATE" : "UPDATE"));
1060   first = false;
1061   f.cd();
1062   TDirectory* d = f.mkdir(Form("%03d", evno),
1063                           Form("Diagnostics histograms for event # %d", evno));
1064   d->cd();
1065   if (fDiagStep1) fDiagStep1->Write();
1066   if (fDiagStep2) fDiagStep2->Write();
1067   if (fDiagStep3) fDiagStep3->Write();
1068   if (fDiagStep4) fDiagStep4->Write();
1069   if (fDiagAll)   fDiagAll->Write();
1070   d->Write();
1071   f.Write();
1072   f.Close();
1073
1074   if (fDiagStep1) fDiagStep1->Reset();
1075   if (fDiagStep2) fDiagStep2->Reset();
1076   if (fDiagStep3) fDiagStep3->Reset();
1077   if (fDiagStep4) fDiagStep4->Reset();
1078   if (fDiagAll)   fDiagAll->Reset();
1079 }
1080
1081 //____________________________________________________________________
1082 void
1083 AliFMDReconstructor::FillESD(AliRawReader*, TTree* clusterTree,
1084                              AliESDEvent* esd) const
1085 {
1086   //
1087   // Forwards to above member function
1088   //
1089   if (fZombie) {
1090     AliWarning("I'm a zombie - cannot do anything");
1091     return;
1092   }
1093   TTree* dummy = 0;
1094   FillESD(dummy, clusterTree, esd);
1095 }
1096 //____________________________________________________________________
1097 //
1098 // EOF
1099 //