FMD/AliFMDReconstructor.cxx

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