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