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           if (fBad(d, r, s, t)) {
 316             AliDebug(5, Form("Marking FMD%d%c[%2d,%3d] as bad", d, r, s, t));
 317             esd->SetMultiplicity(d, r, s, t, AliESDFMD::kInvalidMult);
 318           }
 319           if (param->IsDead(d, r, s, t)) {
 320             AliDebug(5, Form("Marking FMD%d%c[%2d,%3d] as dead", d, r, s, t));
 321             esd->SetMultiplicity(d, r, s, t, AliESDFMD::kInvalidMult);
 322             // esd->SetEta(d, r, s, t, AliESDFMD::kInvalidEta);
 323           }
 324           else if (esd->Multiplicity(d, r, s, t) == AliESDFMD::kInvalidMult) {
 325             AliDebug(10, Form("Setting null signal in FMD%d%c[%2d,%3d]",
 326                               d, r, s, t));
 327             esd->SetMultiplicity(d, r, s, t, 0);
 328           }
 329         }
 330       }
 331     }
 332   }
 333 }
 334
 335 //____________________________________________________________________
 336 void
 337 AliFMDReconstructor::Reconstruct(AliFMDRawReader& rawReader) const
 338 {
 339   AliFMDDebug(1, ("Reconstructing from FMD raw reader"));
 340   if (fZombie {
 341     AliWarning("I'm a zombie - cannot do anything");
 342     return;
 343   }
 344   fBad.Reset(false);
 345   UShort_t det, sec, str, fac;
 346   Short_t  adc, oldDet = -1;
 347   Bool_t   zs;
 348   Char_t   rng;
 349
 350   UseRecoParam(kTRUE);
 351   while (rawReader.NextSignal(det, rng, sec, str, adc, zs, fac)) {
 352     if (det != oldDet) {
 353       fZS[det-1]       = zs;
 354       fZSFactor[det-1] = fac;
 355       oldDet           = det;
 356     }
 357     ProcessSignal(det, rng, sec, str, adc);
 358   }
 359   UseRecoParam(kFALSE);
 360
 361 }
 362
 363 //____________________________________________________________________
 364 void
 365 AliFMDReconstructor::Reconstruct(AliRawReader* reader, TTree*) const
 366 {
 367   // Reconstruct directly from raw data (no intermediate output on
 368   // digit tree or rec point tree).
 369   //
 370   // Parameters:
 371   //   reader   Raw event reader
 372   //   ctree    Not used - 'cluster tree' to store rec-points in.
 373   AliFMDDebug(1, ("Reconstructing from raw reader"));
 374   if (fZombie {
 375     AliWarning("I'm a zombie - cannot do anything");
 376     return;
 377   }
 378   AliFMDRawReader rawReader(reader, 0);
 379   Reconstruct(rawReader);
 380 }
 381
 382 //____________________________________________________________________
 383 void
 384 AliFMDReconstructor::Digitize(AliRawReader* reader, TClonesArray* sdigits) const
 385 {
 386   // Reconstruct directly from raw data (no intermediate output on
 387   // digit tree or rec point tree).
 388   //
 389   // Parameters:
 390   //   reader   Raw event reader
 391   //   ctree    Not used.
 392   if (fZombie {
 393     AliWarning("I'm a zombie - cannot do anything");
 394     return;
 395   }
 396   AliFMDRawReader rawReader(reader, 0);
 397
 398   UShort_t det, sec, str, sam, rat, fac;
 399   Short_t  adc, oldDet = -1;
 400   Bool_t   zs;
 401   Char_t   rng;
 402
 403   UseRecoParam(kTRUE);
 404   while (rawReader.NextSample(det, rng, sec, str, sam, rat, adc, zs, fac)) {
 405     if (!rawReader.SelectSample(sam, rat)) continue;
 406     if (det != oldDet) {
 407       fZS[det-1]       = zs;
 408       fZSFactor[det-1] = fac;
 409       oldDet           = det;
 410     }
 411     DigitizeSignal(sdigits, det, rng, sec, str, sam, adc);
 412   }
 413   UseRecoParam(kFALSE);
 414 }
 415
 416 //____________________________________________________________________
 417 void
 418 AliFMDReconstructor::Reconstruct(TTree* digitsTree,
 419                                  TTree* clusterTree) const
 420 {
 421   // Reconstruct event from digits in tree
 422   // Get the FMD branch holding the digits.
 423   // FIXME: The vertex may not be known yet, so we may have to move
 424   // some of this to FillESD.
 425   //
 426   // Parameters:
 427   //   digitsTree       Pointer to a tree containing digits
 428   //   clusterTree      Pointer to output tree
 429   //
 430   if (fZombie {
 431     AliWarning("I'm a zombie - cannot do anything");
 432     return;
 433   }
 434   if (!fMult) fMult = new TClonesArray("AliFMDRecPoint");
 435
 436   AliFMDDebug(1, ("Reconstructing from digits in a tree"));
 437   GetVertex(fESD);
 438
 439
 440
 441   static TClonesArray* digits = new TClonesArray("AliFMDDigit");
 442   TBranch*      digitBranch   = digitsTree->GetBranch("FMD");
 443   if (!digitBranch) {
 444     Error("Exec", "No digit branch for the FMD found");
 445     return;
 446   }
 447   digitBranch->SetAddress(&digits);
 448
 449   if (digits)  digits->Clear();
 450   if (fMult)   fMult->Clear();
 451   if (fESDObj) fESDObj->Clear();
 452
 453   fNMult = 0;
 454   fTreeR = clusterTree;
 455   fTreeR->Branch("FMD", &fMult);
 456
 457   AliFMDDebug(5, ("Getting entry 0 from digit branch"));
 458   digitBranch->GetEntry(0);
 459
 460   AliFMDDebug(5, ("Processing digits"));
 461   UseRecoParam(kTRUE);
 462   ProcessDigits(digits);
 463   UseRecoParam(kFALSE);
 464
 465   Int_t written = clusterTree->Fill();
 466   AliFMDDebug(10, ("Filled %d bytes into cluster tree", written));
 467   // digits->Delete();
 468   // delete digits;
 469 }
 470
 471
 472 //____________________________________________________________________
 473 void
 474 AliFMDReconstructor::ProcessDigits(TClonesArray* digits,
 475                                    const AliFMDRawReader& rawRead) const
 476 {
 477   // For each digit, find the pseudo rapdity, azimuthal angle, and
 478   // number of corrected ADC counts, and pass it on to the algorithms
 479   // used.
 480   //
 481   // Parameters:
 482   //    digits  Array of digits
 483   //
 484   if (fZombie {
 485     AliWarning("I'm a zombie - cannot do anything");
 486     return;
 487   }
 488   AliFMDAltroMapping* map = AliFMDParameters::Instance()->GetAltroMap();
 489   for (size_t i = 1; i <= 3; i++) {
 490     fZS[i-1]       = rawRead.IsZeroSuppressed(map->Detector2DDL(i));
 491     fZSFactor[i-1] = rawRead.NoiseFactor(map->Detector2DDL(i));
 492   }
 493   UseRecoParam(kTRUE);
 494   ProcessDigits(digits);
 495   UseRecoParam(kFALSE);
 496 }
 497
 498 //____________________________________________________________________
 499 void
 500 AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
 501 {
 502   // For each digit, find the pseudo rapdity, azimuthal angle, and
 503   // number of corrected ADC counts, and pass it on to the algorithms
 504   // used.
 505   //
 506   // Parameters:
 507   //    digits  Array of digits
 508   //
 509   Int_t nDigits = digits->GetEntries();
 510   AliFMDDebug(2, ("Got %d digits", nDigits));
 511   fESDObj->SetNoiseFactor(fNoiseFactor);
 512   fESDObj->SetAngleCorrected(fAngleCorrect);
 513   fBad.Reset(false);
 514   for (Int_t i = 0; i < nDigits; i++) {
 515     AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
 516     if (!digit) continue;
 517     ProcessDigit(digit);
 518   }
 519 }
 520
 521 //____________________________________________________________________
 522 void
 523 AliFMDReconstructor::ProcessDigit(AliFMDDigit* digit) const
 524 {
 525   //
 526   // Process a single digit
 527   //
 528   // Parameters:
 529   //    digit Digiti to process
 530   //
 531   UShort_t det = digit->Detector();
 532   Char_t   rng = digit->Ring();
 533   UShort_t sec = digit->Sector();
 534   UShort_t str = digit->Strip();
 535   Short_t  adc = digit->Counts();
 536
 537   ProcessSignal(det, rng, sec, str, adc);
 538 }
 539
 540 //____________________________________________________________________
 541 void
 542 AliFMDReconstructor::ProcessSignal(UShort_t det,
 543                                    Char_t   rng,
 544                                    UShort_t sec,
 545                                    UShort_t str,
 546                                    Short_t  adc) const
 547 {
 548   // Process the signal from a single strip
 549   //
 550   // Parameters:
 551   //    det     Detector ID
 552   //    rng     Ring ID
 553   //    sec     Sector ID
 554   //    rng     Strip ID
 555   //    adc     ADC counts
 556   //
 557   if (adc >= AliFMDRawReader::kBadSignal) {
 558     AliFMDDebug(1, ("FMD%d%c[%2d,%3d] is marked bad", det, rng, sec, str));
 559     fBad(det,rng,sec,str) = true;
 560     return;
 561   }
 562
 563   // Check that the strip is not marked as dead
 564   AliFMDParameters* param  = AliFMDParameters::Instance();
 565   if (param->IsDead(det, rng, sec, str)) {
 566     AliFMDDebug(1, ("FMD%d%c[%2d,%3d] is dead", det, rng, sec, str));
 567     fBad(det,rng,sec,str) = true;
 568     return;
 569   }
 570
 571   // digit->Print();
 572   // Get eta and phi
 573   Float_t eta, phi;
 574   PhysicalCoordinates(det, rng, sec, str, eta, phi);
 575
 576   // Substract pedestal.
 577   UShort_t counts   = SubtractPedestal(det, rng, sec, str, adc);
 578   if(counts == USHRT_MAX) return;
 579
 580     // Gain match digits.
 581   Double_t edep     = Adc2Energy(det, rng, sec, str, eta, counts);
 582   // Get rid of nonsense energy
 583   if(edep < 0)  return;
 584
 585   // Make rough multiplicity
 586   Double_t mult     = Energy2Multiplicity(det, rng, sec, str, edep);
 587   // Get rid of nonsense mult
 588   //if (mult > 20) {
 589   //  AliWarning(Form("The mutliplicity in FMD%d%c[%2d,%3d]=%f > 20 "
 590   //                "(ADC: %d, Energy: %f)", det, rng, sec, str, mult,
 591   //                counts, edep));
 592   // }
 593   if (mult < 0)  return;
 594   AliFMDDebug(10, ("FMD%d%c[%2d,%3d]: "
 595                     "ADC: %d, Counts: %d, Energy: %f, Mult: %f",
 596                   det, rng, sec, str, adc, counts, edep, mult));
 597
 598   // Create a `RecPoint' on the output branch.
 599   if (fMult) {
 600     AliFMDRecPoint* m =
 601       new ((*fMult)[fNMult]) AliFMDRecPoint(det, rng, sec, str,
 602                                             eta, phi, edep, mult);
 603     (void)m; // Suppress warnings about unused variables.
 604     fNMult++;
 605   }
 606
 607   fESDObj->SetMultiplicity(det, rng, sec, str, mult);
 608   fESDObj->SetEta(det, rng, sec, str, eta);
 609
 610   if (fDiagAll) fDiagAll->Fill(adc, mult);
 611
 612 }
 613
 614 //____________________________________________________________________
 615 void
 616 AliFMDReconstructor::DigitizeSignal(TClonesArray* sdigits,
 617                                     UShort_t det,
 618                                     Char_t   rng,
 619                                     UShort_t sec,
 620                                     UShort_t str,
 621                                     UShort_t /* sam */,
 622                                     Short_t  adc) const
 623 {
 624   // Process the signal from a single strip
 625   //
 626   // Parameters:
 627   //    det     Detector ID
 628   //    rng     Ring ID
 629   //    sec     Sector ID
 630   //    rng     Strip ID
 631   //    adc     ADC counts
 632   //
 633   AliFMDParameters* param  = AliFMDParameters::Instance();
 634   // Check that the strip is not marked as dead
 635   if (param->IsDead(det, rng, sec, str)) {
 636     AliFMDDebug(10, ("FMD%d%c[%2d,%3d] is dead", det, rng, sec, str));
 637     return;
 638   }
 639
 640   // Substract pedestal.
 641   UShort_t counts   = SubtractPedestal(det, rng, sec, str, adc);
 642   if(counts == USHRT_MAX || counts == 0) return;
 643
 644     // Gain match digits.
 645   Double_t edep     = Adc2Energy(det, rng, sec, str, counts);
 646   // Get rid of nonsense energy
 647   if(edep < 0)  return;
 648
 649   Int_t n = sdigits->GetEntriesFast();
 650   // AliFMDSDigit* sdigit =
 651   new ((*sdigits)[n])
 652     AliFMDSDigit(det, rng, sec, str, edep, counts, counts, counts, counts);
 653   // sdigit->SetCount(sam, counts);
 654 }
 655
 656 //____________________________________________________________________
 657 UShort_t
 658 AliFMDReconstructor::SubtractPedestal(UShort_t det,
 659                                       Char_t   rng,
 660                                       UShort_t sec,
 661                                       UShort_t str,
 662                                       UShort_t adc,
 663                                       Float_t  noiseFactor,
 664                                       Bool_t   zsEnabled,
 665                                       UShort_t zsNoiseFactor) const
 666 {
 667   //
 668   // Subtract the pedestal off the ADC counts.
 669   //
 670   // Parameters:
 671   //    det           Detector number
 672   //    rng           Ring identifier
 673   //    sec           Sector number
 674   //    str           Strip number
 675   //    adc           ADC counts
 676   //    noiseFactor   If pedestal substracted pedestal is less then
 677   //        this times the noise, then consider this to be 0.
 678   //    zsEnabled     Whether zero-suppression is on.
 679   //    zsNoiseFactor Noise factor used in on-line pedestal
 680   //        subtraction.
 681   //
 682   // Return:
 683   //    The pedestal subtracted ADC counts (possibly 0), or @c
 684   //         USHRT_MAX in case of problems.
 685   //
 686   AliFMDParameters* param  = AliFMDParameters::Instance();
 687   Float_t           ped    = (zsEnabled ? 0 :
 688                                 param->GetPedestal(det, rng, sec, str));
 689   Float_t           noise  = param->GetPedestalWidth(det, rng, sec, str);
 690   if(ped < 0 || noise < 0) {
 691     AliWarningClass(Form("Invalid pedestal (%f) or noise (%f) "
 692                          "for FMD%d%c[%02d,%03d]",
 693                     ped, noise, det, rng, sec, str));
 694     return USHRT_MAX;
 695   }
 696   AliDebugClass(15, Form("Subtracting pedestal for FMD%d%c[%2d,%3d]=%4d "
 697                          "(%s w/factor %d, noise factor %f, "
 698                          "pedestal %8.2f+/-%8.2f)",
 699                          det, rng, sec, str, adc,
 700                          (zsEnabled ? "zs'ed" : "straight"),
 701                          zsNoiseFactor, noiseFactor, ped, noise));
 702
 703   Int_t counts = adc + Int_t(zsEnabled ? zsNoiseFactor * noise : - ped);
 704   counts =  TMath::Max(Int_t(counts), 0);
 705   // Calculate the noise factor for suppressing remenants of the noise
 706   // peak.  If we have done on-line zero suppression, we only check
 707   // for noise signals that are larger than the suppressed noise.  If
 708   // the noise factor used on line is larger than the factor used
 709   // here, we do not do this check at all.
 710   //
 711   // For example:
 712   //    Online factor  |  Read factor |  Result
 713   //    ---------------+--------------+-------------------------------
 714   //           2       |      3       | Check if signal > 1 * noise
 715   //           3       |      3       | Check if signal > 0
 716   //           3       |      2       | Check if signal > 0
 717   //
 718   // In this way, we make sure that we do not suppress away too much
 719   // data, and that the read-factor is the most stringent cut.
 720   Float_t nf = TMath::Max(0.F, noiseFactor - (zsEnabled ? zsNoiseFactor : 0));
 721   if (counts < noise * nf) counts = 0;
 722   if (counts > 0) AliDebugClass(15, "Got a hit strip");
 723
 724   UShort_t ret = counts < 0 ? 0 : counts;
 725   return ret;
 726 }
 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                                       Short_t  adc) const
 736 {
 737   // Member function to subtract the pedestal from a digit
 738   //
 739   // Parameters:
 740   //    det     Detector ID
 741   //    rng     Ring ID
 742   //    sec     Sector ID
 743   //    rng     Strip ID
 744   //    adc     # of ADC counts
 745   // Return:
 746   //    Pedestal subtracted signal or USHRT_MAX in case of problems
 747   //
 748   UShort_t counts = SubtractPedestal(det, rng, sec, str, adc,
 749                                      fNoiseFactor, fZS[det-1],
 750                                      fZSFactor[det-1]);
 751   if (fDiagStep1) fDiagStep1->Fill(adc, counts);
 752
 753   return counts;
 754 }
 755
 756 //____________________________________________________________________
 757 Float_t
 758 AliFMDReconstructor::Adc2Energy(UShort_t det,
 759                                 Char_t   rng,
 760                                 UShort_t sec,
 761                                 UShort_t str,
 762                                 UShort_t count) const
 763 {
 764   // Converts number of ADC counts to energy deposited.
 765   // Note, that this member function can be overloaded by derived
 766   // classes to do strip-specific look-ups in databases or the like,
 767   // to find the proper gain for a strip.
 768   //
 769   // In the first simple version, we calculate the energy deposited as
 770   //
 771   //    EnergyDeposited = cos(theta) * gain * count
 772   //
 773   // where
 774   //
 775   //           Pre_amp_MIP_Range
 776   //    gain = ----------------- * Energy_deposited_per_MIP
 777   //           ADC_channel_size
 778   //
 779   // is constant and the same for all strips.
 780   //
 781   // For the production we use the conversion measured in the NBI lab.
 782   // The total conversion is then:
 783   //
 784   //    gain = ADC / DAC
 785   //
 786   //                  EdepMip * count
 787   //      => energy = ----------------
 788   //                  gain * DACPerADC
 789   //
 790   // Parameters:
 791   //    det     Detector ID
 792   //    rng     Ring ID
 793   //    sec     Sector ID
 794   //    rng     Strip ID
 795   //    counts  Number of ADC counts over pedestal
 796   // Return
 797   //    The energy deposited in a single strip, or -1 in case of problems
 798   //
 799   if (count <= 0) return 0;
 800   AliFMDParameters* param = AliFMDParameters::Instance();
 801   Float_t           gain  = param->GetPulseGain(det, rng, sec, str);
 802   // 'Tagging' bad gains as bad energy
 803   if (gain < 0) {
 804     AliWarning(Form("Invalid gain (%f) for FMD%d%c[%02d,%03d]",
 805                     gain, det, rng, sec, str));
 806     return -1;
 807   }
 808   AliFMDDebug(5, ("Converting counts %d to energy (factor=%f, DAC2MIP=%f)",
 809                   count, gain,param->GetDACPerMIP()));
 810
 811   Double_t edep  = ((count * param->GetEdepMip())
 812                     / (gain * param->GetDACPerMIP()));
 813   return edep;
 814 }
 815
 816 //____________________________________________________________________
 817 Float_t
 818 AliFMDReconstructor::Adc2Energy(UShort_t det,
 819                                 Char_t   rng,
 820                                 UShort_t sec,
 821                                 UShort_t str,
 822                                 Float_t  eta,
 823                                 UShort_t count) const
 824 {
 825   // Converts number of ADC counts to energy deposited.
 826   // Note, that this member function can be overloaded by derived
 827   // classes to do strip-specific look-ups in databases or the like,
 828   // to find the proper gain for a strip.
 829   //
 830   // In the first simple version, we calculate the energy deposited as
 831   //
 832   //    EnergyDeposited = cos(theta) * gain * count
 833   //
 834   // where
 835   //
 836   //           Pre_amp_MIP_Range
 837   //    gain = ----------------- * Energy_deposited_per_MIP
 838   //           ADC_channel_size
 839   //
 840   // is constant and the same for all strips.
 841   //
 842   // For the production we use the conversion measured in the NBI lab.
 843   // The total conversion is then:
 844   //
 845   //    gain = ADC / DAC
 846   //
 847   //                  EdepMip * count
 848   //      => energy = ----------------
 849   //                  gain * DACPerADC
 850   //
 851   // Parameters:
 852   //    det     Detector ID
 853   //    rng     Ring ID
 854   //    sec     Sector ID
 855   //    rng     Strip ID
 856   //    eta     Psuedo-rapidity
 857   //    counts  Number of ADC counts over pedestal
 858   // Return
 859   //    The energy deposited in a single strip, or -1 in case of problems
 860   //
 861   Double_t edep = Adc2Energy(det, rng, sec, str, count);
 862
 863   if (fDiagStep2) fDiagStep2->Fill(count, edep);
 864   if (fAngleCorrect) {
 865     Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
 866     Double_t corr  = TMath::Abs(TMath::Cos(theta));
 867     Double_t cedep = corr * edep;
 868     AliFMDDebug(10, ("correcting for path %f * %f = %f (eta=%f, theta=%f)",
 869                       edep, corr, cedep, eta, theta));
 870     if (fDiagStep3) fDiagStep3->Fill(edep, cedep);
 871     edep = cedep;
 872   }
 873   return edep;
 874 }
 875
 876 //____________________________________________________________________
 877 Float_t
 878 AliFMDReconstructor::Energy2Multiplicity(UShort_t /*det*/,
 879                                          Char_t   /*rng*/,
 880                                          UShort_t /*sec*/,
 881                                          UShort_t /*str*/,
 882                                          Float_t  edep) const
 883 {
 884   // Converts an energy signal to number of particles.
 885   // Note, that this member function can be overloaded by derived
 886   // classes to do strip-specific look-ups in databases or the like,
 887   // to find the proper gain for a strip.
 888   //
 889   // In this simple version, we calculate the multiplicity as
 890   //
 891   //   multiplicity = Energy_deposited / Energy_deposited_per_MIP
 892   //
 893   // where
 894   //
 895   //   Energy_deposited_per_MIP = 1.664 * SI_density * SI_thickness
 896   //
 897   // is constant and the same for all strips
 898   //
 899   // Parameters:
 900   //    det     Detector ID
 901   //    rng     Ring ID
 902   //    sec     Sector ID
 903   //    rng     Strip ID
 904   //    edep    Energy deposited in a single strip
 905   // Return
 906   //    The "bare" multiplicity corresponding to the energy deposited
 907   AliFMDParameters* param   = AliFMDParameters::Instance();
 908   Double_t          edepMIP = param->GetEdepMip();
 909   Float_t           mult    = edep / edepMIP;
 910 #if 0
 911   if (edep > 0)
 912     AliFMDDebug(15, ("Translating energy %f to multiplicity via "
 913                      "divider %f->%f", edep, edepMIP, mult));
 914 #endif
 915   if (fDiagStep4) fDiagStep4->Fill(edep, mult);
 916   return mult;
 917 }
 918
 919 //____________________________________________________________________
 920 void
 921 AliFMDReconstructor::PhysicalCoordinates(UShort_t det,
 922                                          Char_t   rng,
 923                                          UShort_t sec,
 924                                          UShort_t str,
 925                                          Float_t& eta,
 926                                          Float_t& phi) const
 927 {
 928   // Get the eta and phi of a digit
 929   //
 930   // Parameters:
 931   //    det     Detector ID
 932   //    rng     Ring ID
 933   //    sec     Sector ID
 934   //    rng     Strip ID
 935   //    eta     On return, contains the psuedo-rapidity of the strip
 936   //    phi     On return, contains the azimuthal angle of the strip
 937   //
 938   AliFMDGeometry* geom = AliFMDGeometry::Instance();
 939   Double_t x, y, z, r, theta, deta, dphi;
 940   geom->Detector2XYZ(det, rng, sec, str, x, y, z);
 941
 942   // Correct for vertex offset.
 943   z     -= fCurrentVertex;
 944   AliFMDGeometry::XYZ2REtaPhiTheta(x, y, z, r, deta, dphi, theta);
 945   eta = deta;
 946   phi = dphi;
 947 }
 948
 949 namespace {
 950   class ESDPrinter : public AliESDFMD::ForOne
 951   {
 952   public:
 953     ESDPrinter() {}
 954     Bool_t operator()(UShort_t d, Char_t r, UShort_t s, UShort_t t,
 955                       Float_t m, Float_t e)
 956     {
 957       if (m > 0 && m != AliESDFMD::kInvalidMult)
 958         printf("  FMD%d%c[%2d,%3d] = %6.3f / %6.3f\n", d, r, s, t, m, e);
 959       return kTRUE;
 960     }
 961   };
 962 }
 963
 964
 965 //____________________________________________________________________
 966 void
 967 AliFMDReconstructor::FillESD(TTree*  /* digitsTree */,
 968                              TTree*  /* clusterTree */,
 969                              AliESDEvent* esd) const
 970 {
 971   // nothing to be done
 972   // FIXME: The vertex may not be known when Reconstruct is executed,
 973   // so we may have to move some of that member function here.
 974   AliFMDDebug(2, ("Calling FillESD with two trees and one ESD"));
 975   if (fZombie {
 976     AliWarning("I'm a zombie - cannot do anything");
 977     return;
 978   }
 979   // fESDObj->Print();
 980
 981   // Fix up ESD so that only truely dead channels get the kInvalidMult flag.
 982   MarkDeadChannels(fESDObj);
 983
 984   Double_t oldVz = fCurrentVertex;
 985   GetVertex(esd);
 986   if (fVertexType != kNoVertex) {
 987     AliFMDDebug(2, ("Revertexing the ESD data to vz=%f (was %f)",
 988                     fCurrentVertex, oldVz));
 989     AliFMDESDRevertexer revertexer;
 990     revertexer.Revertex(fESDObj, fCurrentVertex);
 991   }
 992
 993   if (AliDebugLevel() > 10) {
 994     ESDPrinter p;
 995     fESDObj->ForEach(p);
 996   }
 997
 998   if (esd) {
 999     AliFMDDebug(2, ("Writing FMD data to ESD tree"));
1000     esd->SetFMDData(fESDObj);
1001   }
1002
1003   if (!fDiagnostics || !esd) return;
1004   static bool first = true;
1005   // This is most likely NOT the event number you'd like to use. It
1006   // has nothing to do with the 'real' event number.
1007   // - That's OK.  We just use it for the name of the directory -
1008   // nothing else.  Christian
1009   Int_t evno = esd->GetEventNumberInFile();
1010   AliFMDDebug(3, ("Writing diagnostics histograms to FMD.Diag.root/%03d",evno));
1011   TFile f("FMD.Diag.root", (first ? "RECREATE" : "UPDATE"));
1012   first = false;
1013   f.cd();
1014   TDirectory* d = f.mkdir(Form("%03d", evno),
1015                           Form("Diagnostics histograms for event # %d", evno));
1016   d->cd();
1017   if (fDiagStep1) fDiagStep1->Write();
1018   if (fDiagStep2) fDiagStep2->Write();
1019   if (fDiagStep3) fDiagStep3->Write();
1020   if (fDiagStep4) fDiagStep4->Write();
1021   if (fDiagAll)   fDiagAll->Write();
1022   d->Write();
1023   f.Write();
1024   f.Close();
1025
1026   if (fDiagStep1) fDiagStep1->Reset();
1027   if (fDiagStep2) fDiagStep2->Reset();
1028   if (fDiagStep3) fDiagStep3->Reset();
1029   if (fDiagStep4) fDiagStep4->Reset();
1030   if (fDiagAll)   fDiagAll->Reset();
1031 }
1032
1033 //____________________________________________________________________
1034 void
1035 AliFMDReconstructor::FillESD(AliRawReader*, TTree* clusterTree,
1036                              AliESDEvent* esd) const
1037 {
1038   //
1039   // Forwards to above member function
1040   //
1041   if (fZombie {
1042     AliWarning("I'm a zombie - cannot do anything");
1043     return;
1044   }
1045   TTree* dummy = 0;
1046   FillESD(dummy, clusterTree, esd);
1047 }
1048 //____________________________________________________________________
1049 //
1050 // EOF
1051 //