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 // Import revertexer into a private namespace (to prevent conflicts)
  57 namespace {
  58 # include "AliFMDESDRevertexer.h"
  59 }
  60
  61
  62 class AliRawReader;
  63
  64 //____________________________________________________________________
  65 ClassImp(AliFMDReconstructor)
  66 #if 0
  67   ; // This is here to keep Emacs for indenting the next line
  68 #endif
  69
  70 //____________________________________________________________________
  71 AliFMDReconstructor::AliFMDReconstructor()
  72   : AliReconstructor(),
  73     fMult(0x0),
  74     fNMult(0),
  75     fTreeR(0x0),
  76     fCurrentVertex(0),
  77     fESDObj(0x0),
  78     fNoiseFactor(0),
  79     fAngleCorrect(kTRUE),
  80     fVertexType(kNoVertex),
  81     fESD(0x0),
  82     fDiagnostics(kTRUE),
  83     fDiagStep1(0),
  84     fDiagStep2(0),
  85     fDiagStep3(0),
  86     fDiagStep4(0),
  87     fDiagAll(0)
  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   param->Init();
 122
 123   // Current vertex position
 124   fCurrentVertex = 0;
 125   // Create array of reconstructed strip multiplicities
 126   fMult = new TClonesArray("AliFMDRecPoint", 51200);
 127   // Create ESD output object
 128   fESDObj = new AliESDFMD;
 129
 130   // Check if we need diagnostics histograms
 131   if (!fDiagnostics) return;
 132   AliInfo("Making diagnostics histograms");
 133   fDiagStep1   = new TH2I("diagStep1", "Read ADC vs. Noise surpressed ADC",
 134                         1024, -.5, 1023.5, 1024, -.5, 1023.5);
 135   fDiagStep1->SetDirectory(0);
 136   fDiagStep1->GetXaxis()->SetTitle("ADC (read)");
 137   fDiagStep1->GetYaxis()->SetTitle(Form("ADC (noise surpressed %4.f)",
 138                                         fNoiseFactor));
 139   fDiagStep2  = new TH2F("diagStep2",  "ADC vs Edep deduced",
 140                         1024, -.5, 1023.5, 100, 0, 2);
 141   fDiagStep2->SetDirectory(0);
 142   fDiagStep2->GetXaxis()->SetTitle("ADC (noise surpressed)");
 143   fDiagStep2->GetYaxis()->SetTitle("#Delta E [GeV]");
 144   fDiagStep3  = new TH2F("diagStep3",  "Edep vs Edep path corrected",
 145                         100, 0., 2., 100, 0., 2.);
 146   fDiagStep3->SetDirectory(0);
 147   fDiagStep3->GetXaxis()->SetTitle("#Delta E [GeV]");
 148   fDiagStep3->GetYaxis()->SetTitle("#Delta E/#Delta x #times #delta x [GeV]");
 149   fDiagStep4  = new TH2F("diagStep4",  "Edep vs Multiplicity deduced",
 150                         100, 0., 2., 100, -.1, 19.9);
 151   fDiagStep4->SetDirectory(0);
 152   fDiagStep4->GetXaxis()->SetTitle("#Delta E/#Delta x #times #delta x [GeV]");
 153   fDiagStep4->GetYaxis()->SetTitle("Multiplicity");
 154   fDiagAll    = new TH2F("diagAll",    "Read ADC vs Multiplicity deduced",
 155                          1024, -.5, 1023.5, 100, -.1, 19.9);
 156   fDiagAll->SetDirectory(0);
 157   fDiagAll->GetXaxis()->SetTitle("ADC (read)");
 158   fDiagAll->GetYaxis()->SetTitle("Multiplicity");
 159 }
 160
 161 //____________________________________________________________________
 162 void
 163 AliFMDReconstructor::ConvertDigits(AliRawReader* reader,
 164                                    TTree* digitsTree) const
 165 {
 166   // Convert Raw digits to AliFMDDigit's in a tree
 167   AliFMDDebug(1, ("Reading raw data into digits tree"));
 168   AliFMDRawReader rawRead(reader, digitsTree);
 169   // rawRead.SetSampleRate(fFMD->GetSampleRate());
 170   rawRead.Exec();
 171   AliFMDAltroMapping* map = AliFMDParameters::Instance()->GetAltroMap();
 172   for (size_t i = 1; i <= 3; i++) {
 173     fZS[i]       = rawRead.IsZeroSuppressed(map->Detector2DDL(i));
 174     fZSFactor[i] = rawRead.NoiseFactor(map->Detector2DDL(i));
 175   }
 176 }
 177
 178 //____________________________________________________________________
 179 void
 180 AliFMDReconstructor::GetVertex(AliESDEvent* esd) const
 181 {
 182   // Return the vertex to use.
 183   // This is obtained from the ESD object.
 184   // If not found, a warning is issued.
 185   fVertexType    = kNoVertex;
 186   fCurrentVertex = 0;
 187   if (!esd) return;
 188
 189   const AliESDVertex* vertex = esd->GetPrimaryVertex();
 190   if (!vertex)        vertex = esd->GetPrimaryVertexSPD();
 191   if (!vertex)        vertex = esd->GetPrimaryVertexTPC();
 192   if (!vertex)        vertex = esd->GetVertex();
 193
 194   if (vertex) {
 195     AliFMDDebug(2, ("Got %s (%s) from ESD: %f",
 196                     vertex->GetName(), vertex->GetTitle(), vertex->GetZv()));
 197     fCurrentVertex = vertex->GetZv();
 198     fVertexType    = kESDVertex;
 199     return;
 200   }
 201   else if (esd->GetESDTZERO()) {
 202     AliFMDDebug(2, ("Got primary vertex from T0: %f", esd->GetT0zVertex()));
 203     fCurrentVertex = esd->GetT0zVertex();
 204     fVertexType    = kESDVertex;
 205     return;
 206   }
 207   AliWarning("Didn't get any vertex from ESD or generator");
 208 }
 209
 210 //____________________________________________________________________
 211 Int_t
 212 AliFMDReconstructor::GetIdentifier() const
 213 {
 214   return AliReconstruction::GetDetIndex(GetDetectorName());
 215 }
 216
 217 //____________________________________________________________________
 218 const AliFMDRecoParam*
 219 AliFMDReconstructor::GetParameters() const
 220 {
 221   Int_t iDet = 12; // GetIdentifier();
 222   const AliDetectorRecoParam* params = AliReconstructor::GetRecoParam(iDet);
 223   if (!params || params->IsA() != AliFMDRecoParam::Class()) return 0;
 224   return static_cast<const AliFMDRecoParam*>(params);
 225 }
 226
 227 //____________________________________________________________________
 228 void
 229 AliFMDReconstructor::UseRecoParam(Bool_t set) const
 230 {
 231   static Float_t savedNoiseFactor  = fNoiseFactor;
 232   static Bool_t  savedAngleCorrect = fAngleCorrect;
 233   if (set) {
 234     const AliFMDRecoParam* params  = GetParameters();
 235     if (params) {
 236       fNoiseFactor  = params->NoiseFactor();
 237       fAngleCorrect = params->AngleCorrect();
 238     }
 239     return;
 240   }
 241   fNoiseFactor  = savedNoiseFactor;
 242   fAngleCorrect = savedAngleCorrect;
 243 }
 244
 245
 246
 247 //____________________________________________________________________
 248 void
 249 AliFMDReconstructor::Reconstruct(AliRawReader* reader, TTree*) const
 250 {
 251   // Reconstruct directly from raw data (no intermediate output on
 252   // digit tree or rec point tree).
 253   //
 254   // Parameters:
 255   //   reader   Raw event reader
 256   //   ctree    Not used.
 257   AliFMDDebug(1, ("Reconstructing from raw reader"));
 258   AliFMDRawReader rawReader(reader, 0);
 259
 260   UShort_t det, sec, str, fac;
 261   Short_t  adc, oldDet = -1;
 262   Bool_t   zs;
 263   Char_t   rng;
 264
 265   UseRecoParam(kTRUE);
 266   while (rawReader.NextSignal(det, rng, sec, str, adc, zs, fac)) {
 267     if (det != oldDet) {
 268       fZS[det-1]       = zs;
 269       fZSFactor[det-1] = fac;
 270       oldDet           = det;
 271     }
 272     ProcessSignal(det, rng, sec, str, adc);
 273   }
 274   UseRecoParam(kFALSE);
 275
 276 }
 277
 278 //____________________________________________________________________
 279 void
 280 AliFMDReconstructor::Digitize(AliRawReader* reader, TClonesArray* sdigits) const
 281 {
 282   // Reconstruct directly from raw data (no intermediate output on
 283   // digit tree or rec point tree).
 284   //
 285   // Parameters:
 286   //   reader   Raw event reader
 287   //   ctree    Not used.
 288   AliFMDRawReader rawReader(reader, 0);
 289
 290   UShort_t det, sec, str, sam, rat, fac;
 291   Short_t  adc, oldDet = -1;
 292   Bool_t   zs;
 293   Char_t   rng;
 294
 295   UseRecoParam(kTRUE);
 296   while (rawReader.NextSample(det, rng, sec, str, sam, rat, adc, zs, fac)) {
 297     if (!rawReader.SelectSample(sam, rat)) continue;
 298     if (det != oldDet) {
 299       fZS[det-1]       = zs;
 300       fZSFactor[det-1] = fac;
 301       oldDet           = det;
 302     }
 303     DigitizeSignal(sdigits, det, rng, sec, str, sam, adc);
 304   }
 305   UseRecoParam(kFALSE);
 306 }
 307
 308 //____________________________________________________________________
 309 void
 310 AliFMDReconstructor::Reconstruct(TTree* digitsTree,
 311                                  TTree* clusterTree) const
 312 {
 313   // Reconstruct event from digits in tree
 314   // Get the FMD branch holding the digits.
 315   // FIXME: The vertex may not be known yet, so we may have to move
 316   // some of this to FillESD.
 317   //
 318   // Parameters:
 319   //   digitsTree       Pointer to a tree containing digits
 320   //   clusterTree      Pointer to output tree
 321   //
 322   AliFMDDebug(1, ("Reconstructing from digits in a tree"));
 323   GetVertex(fESD);
 324
 325   TBranch *digitBranch = digitsTree->GetBranch("FMD");
 326   if (!digitBranch) {
 327     Error("Exec", "No digit branch for the FMD found");
 328     return;
 329   }
 330   TClonesArray* digits = new TClonesArray("AliFMDDigit");
 331   digitBranch->SetAddress(&digits);
 332
 333   if (fMult)   fMult->Clear();
 334   if (fESDObj) fESDObj->Clear();
 335
 336   fNMult = 0;
 337   fTreeR = clusterTree;
 338   fTreeR->Branch("FMD", &fMult);
 339
 340   AliFMDDebug(5, ("Getting entry 0 from digit branch"));
 341   digitBranch->GetEntry(0);
 342
 343   AliFMDDebug(5, ("Processing digits"));
 344   UseRecoParam(kTRUE);
 345   ProcessDigits(digits);
 346   UseRecoParam(kFALSE);
 347
 348   Int_t written = clusterTree->Fill();
 349   AliFMDDebug(10, ("Filled %d bytes into cluster tree", written));
 350   digits->Delete();
 351   delete digits;
 352 }
 353
 354
 355 //____________________________________________________________________
 356 void
 357 AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
 358 {
 359   // For each digit, find the pseudo rapdity, azimuthal angle, and
 360   // number of corrected ADC counts, and pass it on to the algorithms
 361   // used.
 362   //
 363   // Parameters:
 364   //    digits  Array of digits
 365   //
 366   Int_t nDigits = digits->GetEntries();
 367   AliFMDDebug(2, ("Got %d digits", nDigits));
 368   fESDObj->SetNoiseFactor(fNoiseFactor);
 369   fESDObj->SetAngleCorrected(fAngleCorrect);
 370   for (Int_t i = 0; i < nDigits; i++) {
 371     AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
 372     if (!digit) continue;
 373     ProcessDigit(digit);
 374   }
 375 }
 376
 377 //____________________________________________________________________
 378 void
 379 AliFMDReconstructor::ProcessDigit(AliFMDDigit* digit) const
 380 {
 381   UShort_t det = digit->Detector();
 382   Char_t   rng = digit->Ring();
 383   UShort_t sec = digit->Sector();
 384   UShort_t str = digit->Strip();
 385   Short_t  adc = digit->Counts();
 386
 387   ProcessSignal(det, rng, sec, str, adc);
 388 }
 389
 390 //____________________________________________________________________
 391 void
 392 AliFMDReconstructor::ProcessSignal(UShort_t det,
 393                                    Char_t   rng,
 394                                    UShort_t sec,
 395                                    UShort_t str,
 396                                    Short_t  adc) const
 397 {
 398   // Process the signal from a single strip
 399   //
 400   // Parameters:
 401   //    det     Detector ID
 402   //    rng     Ring ID
 403   //    sec     Sector ID
 404   //    rng     Strip ID
 405   //    adc     ADC counts
 406   //
 407   AliFMDParameters* param  = AliFMDParameters::Instance();
 408   // Check that the strip is not marked as dead
 409   if (param->IsDead(det, rng, sec, str)) {
 410     AliFMDDebug(3, ("FMD%d%c[%2d,%3d] is dead", det, rng, sec, str));
 411     return;
 412   }
 413
 414   // digit->Print();
 415   // Get eta and phi
 416   Float_t eta, phi;
 417   PhysicalCoordinates(det, rng, sec, str, eta, phi);
 418
 419   // Substract pedestal.
 420   UShort_t counts   = SubtractPedestal(det, rng, sec, str, adc);
 421   if(counts == USHRT_MAX) return;
 422
 423     // Gain match digits.
 424   Double_t edep     = Adc2Energy(det, rng, sec, str, eta, counts);
 425   // Get rid of nonsense energy
 426   if(edep < 0)  return;
 427
 428   // Make rough multiplicity
 429   Double_t mult     = Energy2Multiplicity(det, rng, sec, str, edep);
 430   // Get rid of nonsense mult
 431   if (mult < 0)  return;
 432   AliFMDDebug(10, ("FMD%d%c[%2d,%3d]: "
 433                     "ADC: %d, Counts: %d, Energy: %f, Mult: %f",
 434                   det, rng, sec, str, adc, counts, edep, mult));
 435
 436   // Create a `RecPoint' on the output branch.
 437   if (fMult) {
 438     AliFMDRecPoint* m =
 439       new ((*fMult)[fNMult]) AliFMDRecPoint(det, rng, sec, str,
 440                                             eta, phi, edep, mult);
 441     (void)m; // Suppress warnings about unused variables.
 442     fNMult++;
 443   }
 444
 445   fESDObj->SetMultiplicity(det, rng, sec, str, mult);
 446   fESDObj->SetEta(det, rng, sec, str, eta);
 447
 448   if (fDiagAll) fDiagAll->Fill(adc, mult);
 449
 450 }
 451
 452 //____________________________________________________________________
 453 void
 454 AliFMDReconstructor::DigitizeSignal(TClonesArray* sdigits,
 455                                     UShort_t det,
 456                                     Char_t   rng,
 457                                     UShort_t sec,
 458                                     UShort_t str,
 459                                     UShort_t /* sam */,
 460                                     Short_t  adc) const
 461 {
 462   // Process the signal from a single strip
 463   //
 464   // Parameters:
 465   //    det     Detector ID
 466   //    rng     Ring ID
 467   //    sec     Sector ID
 468   //    rng     Strip ID
 469   //    adc     ADC counts
 470   //
 471   AliFMDParameters* param  = AliFMDParameters::Instance();
 472   // Check that the strip is not marked as dead
 473   if (param->IsDead(det, rng, sec, str)) {
 474     AliFMDDebug(10, ("FMD%d%c[%2d,%3d] is dead", det, rng, sec, str));
 475     return;
 476   }
 477
 478   // Substract pedestal.
 479   UShort_t counts   = SubtractPedestal(det, rng, sec, str, adc);
 480   if(counts == USHRT_MAX || counts == 0) return;
 481
 482     // Gain match digits.
 483   Double_t edep     = Adc2Energy(det, rng, sec, str, counts);
 484   // Get rid of nonsense energy
 485   if(edep < 0)  return;
 486
 487   Int_t n = sdigits->GetEntriesFast();
 488   // AliFMDSDigit* sdigit =
 489   new ((*sdigits)[n])
 490     AliFMDSDigit(det, rng, sec, str, edep, counts, counts, counts, counts);
 491   // sdigit->SetCount(sam, counts);
 492 }
 493
 494 //____________________________________________________________________
 495 UShort_t
 496 AliFMDReconstructor::SubtractPedestal(UShort_t det,
 497                                       Char_t   rng,
 498                                       UShort_t sec,
 499                                       UShort_t str,
 500                                       UShort_t adc,
 501                                       Float_t  noiseFactor,
 502                                       Bool_t   zsEnabled,
 503                                       UShort_t zsNoiseFactor) const
 504 {
 505   AliFMDParameters* param  = AliFMDParameters::Instance();
 506   Float_t           ped    = (zsEnabled ? 0 :
 507                                 param->GetPedestal(det, rng, sec, str));
 508   Float_t           noise  = param->GetPedestalWidth(det, rng, sec, str);
 509   if(ped < 0 || noise < 0) {
 510     AliWarningClass(Form("Invalid pedestal (%f) or noise (%f) "
 511                          "for FMD%d%c[%02d,%03d]",
 512                     ped, noise, det, rng, sec, str));
 513     return USHRT_MAX;
 514   }
 515   AliDebugClass(15, Form("Subtracting pedestal for FMD%d%c[%2d,%3d]=%4d "
 516                          "(%s w/factor %d, noise factor %f, "
 517                          "pedestal %8.2f+/-%8.2f)",
 518                          det, rng, sec, str, adc,
 519                          (zsEnabled ? "zs'ed" : "straight"),
 520                          zsNoiseFactor, noiseFactor, ped, noise));
 521
 522   Int_t counts = adc + Int_t(zsEnabled ? zsNoiseFactor * noise : - ped);
 523   counts =  TMath::Max(Int_t(counts), 0);
 524   // Calculate the noise factor for suppressing remenants of the noise
 525   // peak.  If we have done on-line zero suppression, we only check
 526   // for noise signals that are larger than the suppressed noise.  If
 527   // the noise factor used on line is larger than the factor used
 528   // here, we do not do this check at all.
 529   //
 530   // For example:
 531   //    Online factor  |  Read factor |  Result
 532   //    ---------------+--------------+-------------------------------
 533   //           2       |      3       | Check if signal > 1 * noise
 534   //           3       |      3       | Check if signal > 0
 535   //           3       |      2       | Check if signal > 0
 536   //
 537   // In this way, we make sure that we do not suppress away too much
 538   // data, and that the read-factor is the most stringent cut.
 539   Float_t nf = TMath::Max(0.F, noiseFactor - (zsEnabled ? zsNoiseFactor : 0));
 540   if (counts < noise * nf) counts = 0;
 541   if (counts > 0) AliDebugClass(15, "Got a hit strip");
 542
 543   return counts;
 544 }
 545
 546
 547 //____________________________________________________________________
 548 UShort_t
 549 AliFMDReconstructor::SubtractPedestal(UShort_t det,
 550                                       Char_t   rng,
 551                                       UShort_t sec,
 552                                       UShort_t str,
 553                                       Short_t  adc) const
 554 {
 555   // Member function to subtract the pedestal from a digit
 556   //
 557   // Parameters:
 558   //    det     Detector ID
 559   //    rng     Ring ID
 560   //    sec     Sector ID
 561   //    rng     Strip ID
 562   //    adc     # of ADC counts
 563   // Return:
 564   //    Pedestal subtracted signal or USHRT_MAX in case of problems
 565   //
 566   UShort_t counts = SubtractPedestal(det, rng, sec, str, adc,
 567                                      fNoiseFactor, fZS[det-1],
 568                                      fZSFactor[det-1]);
 569   if (fDiagStep1) fDiagStep1->Fill(adc, counts);
 570
 571   return counts;
 572 }
 573
 574 //____________________________________________________________________
 575 Float_t
 576 AliFMDReconstructor::Adc2Energy(UShort_t det,
 577                                 Char_t   rng,
 578                                 UShort_t sec,
 579                                 UShort_t str,
 580                                 UShort_t count) const
 581 {
 582   // Converts number of ADC counts to energy deposited.
 583   // Note, that this member function can be overloaded by derived
 584   // classes to do strip-specific look-ups in databases or the like,
 585   // to find the proper gain for a strip.
 586   //
 587   // In the first simple version, we calculate the energy deposited as
 588   //
 589   //    EnergyDeposited = cos(theta) * gain * count
 590   //
 591   // where
 592   //
 593   //           Pre_amp_MIP_Range
 594   //    gain = ----------------- * Energy_deposited_per_MIP
 595   //           ADC_channel_size
 596   //
 597   // is constant and the same for all strips.
 598   //
 599   // For the production we use the conversion measured in the NBI lab.
 600   // The total conversion is then:
 601   //
 602   //    gain = ADC / DAC
 603   //
 604   //                  EdepMip * count
 605   //      => energy = ----------------
 606   //                  gain * DACPerADC
 607   //
 608   // Parameters:
 609   //    det     Detector ID
 610   //    rng     Ring ID
 611   //    sec     Sector ID
 612   //    rng     Strip ID
 613   //    counts  Number of ADC counts over pedestal
 614   // Return
 615   //    The energy deposited in a single strip, or -1 in case of problems
 616   //
 617   if (count <= 0) return 0;
 618   AliFMDParameters* param = AliFMDParameters::Instance();
 619   Float_t           gain  = param->GetPulseGain(det, rng, sec, str);
 620   // 'Tagging' bad gains as bad energy
 621   if (gain < 0) {
 622     AliWarning(Form("Invalid gain (%f) for FMD%d%c[%02d,%03d]",
 623                     gain, det, rng, sec, str));
 624     return -1;
 625   }
 626   AliFMDDebug(5, ("Converting counts %d to energy (factor=%f, DAC2MIP=%f)",
 627                   count, gain,param->GetDACPerMIP()));
 628
 629   Double_t edep  = ((count * param->GetEdepMip())
 630                     / (gain * param->GetDACPerMIP()));
 631   return edep;
 632 }
 633
 634 //____________________________________________________________________
 635 Float_t
 636 AliFMDReconstructor::Adc2Energy(UShort_t det,
 637                                 Char_t   rng,
 638                                 UShort_t sec,
 639                                 UShort_t str,
 640                                 Float_t  eta,
 641                                 UShort_t count) const
 642 {
 643   // Converts number of ADC counts to energy deposited.
 644   // Note, that this member function can be overloaded by derived
 645   // classes to do strip-specific look-ups in databases or the like,
 646   // to find the proper gain for a strip.
 647   //
 648   // In the first simple version, we calculate the energy deposited as
 649   //
 650   //    EnergyDeposited = cos(theta) * gain * count
 651   //
 652   // where
 653   //
 654   //           Pre_amp_MIP_Range
 655   //    gain = ----------------- * Energy_deposited_per_MIP
 656   //           ADC_channel_size
 657   //
 658   // is constant and the same for all strips.
 659   //
 660   // For the production we use the conversion measured in the NBI lab.
 661   // The total conversion is then:
 662   //
 663   //    gain = ADC / DAC
 664   //
 665   //                  EdepMip * count
 666   //      => energy = ----------------
 667   //                  gain * DACPerADC
 668   //
 669   // Parameters:
 670   //    det     Detector ID
 671   //    rng     Ring ID
 672   //    sec     Sector ID
 673   //    rng     Strip ID
 674   //    eta     Psuedo-rapidity
 675   //    counts  Number of ADC counts over pedestal
 676   // Return
 677   //    The energy deposited in a single strip, or -1 in case of problems
 678   //
 679   Double_t edep = Adc2Energy(det, rng, sec, str, count);
 680
 681   if (fDiagStep2) fDiagStep2->Fill(count, edep);
 682   if (fAngleCorrect) {
 683     Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
 684     Double_t corr  = TMath::Abs(TMath::Cos(theta));
 685     Double_t cedep = corr * edep;
 686     AliFMDDebug(10, ("correcting for path %f * %f = %f (eta=%f, theta=%f)",
 687                       edep, corr, cedep, eta, theta));
 688     if (fDiagStep3) fDiagStep3->Fill(edep, cedep);
 689     edep = cedep;
 690   }
 691   return edep;
 692 }
 693
 694 //____________________________________________________________________
 695 Float_t
 696 AliFMDReconstructor::Energy2Multiplicity(UShort_t /*det*/,
 697                                          Char_t   /*rng*/,
 698                                          UShort_t /*sec*/,
 699                                          UShort_t /*str*/,
 700                                          Float_t  edep) const
 701 {
 702   // Converts an energy signal to number of particles.
 703   // Note, that this member function can be overloaded by derived
 704   // classes to do strip-specific look-ups in databases or the like,
 705   // to find the proper gain for a strip.
 706   //
 707   // In this simple version, we calculate the multiplicity as
 708   //
 709   //   multiplicity = Energy_deposited / Energy_deposited_per_MIP
 710   //
 711   // where
 712   //
 713   //   Energy_deposited_per_MIP = 1.664 * SI_density * SI_thickness
 714   //
 715   // is constant and the same for all strips
 716   //
 717   // Parameters:
 718   //    det     Detector ID
 719   //    rng     Ring ID
 720   //    sec     Sector ID
 721   //    rng     Strip ID
 722   //    edep    Energy deposited in a single strip
 723   // Return
 724   //    The "bare" multiplicity corresponding to the energy deposited
 725   AliFMDParameters* param   = AliFMDParameters::Instance();
 726   Double_t          edepMIP = param->GetEdepMip();
 727   Float_t           mult    = edep / edepMIP;
 728 #if 0
 729   if (edep > 0)
 730     AliFMDDebug(15, ("Translating energy %f to multiplicity via "
 731                      "divider %f->%f", edep, edepMIP, mult));
 732 #endif
 733   if (fDiagStep4) fDiagStep4->Fill(edep, mult);
 734   return mult;
 735 }
 736
 737 //____________________________________________________________________
 738 void
 739 AliFMDReconstructor::PhysicalCoordinates(UShort_t det,
 740                                          Char_t   rng,
 741                                          UShort_t sec,
 742                                          UShort_t str,
 743                                          Float_t& eta,
 744                                          Float_t& phi) const
 745 {
 746   // Get the eta and phi of a digit
 747   //
 748   // Parameters:
 749   //    det     Detector ID
 750   //    rng     Ring ID
 751   //    sec     Sector ID
 752   //    rng     Strip ID
 753   //    eta     On return, contains the psuedo-rapidity of the strip
 754   //    phi     On return, contains the azimuthal angle of the strip
 755   //
 756   AliFMDGeometry* geom = AliFMDGeometry::Instance();
 757   Double_t x, y, z, r, theta;
 758   geom->Detector2XYZ(det, rng, sec, str, x, y, z);
 759   // Correct for vertex offset.
 760   z     -= fCurrentVertex;
 761   phi   =  TMath::ATan2(y, x);
 762   r     =  TMath::Sqrt(y * y + x * x);
 763   theta =  TMath::ATan2(r, z);
 764   eta   = -TMath::Log(TMath::Tan(theta / 2));
 765 }
 766
 767
 768
 769 //____________________________________________________________________
 770 void
 771 AliFMDReconstructor::FillESD(TTree*  /* digitsTree */,
 772                              TTree*  /* clusterTree */,
 773                              AliESDEvent* esd) const
 774 {
 775   // nothing to be done
 776   // FIXME: The vertex may not be known when Reconstruct is executed,
 777   // so we may have to move some of that member function here.
 778   AliFMDDebug(2, ("Calling FillESD with two trees and one ESD"));
 779   // fESDObj->Print();
 780
 781   Double_t oldVz = fCurrentVertex;
 782   GetVertex(esd);
 783   if (fVertexType != kNoVertex) {
 784     AliFMDDebug(2, ("Revertexing the ESD data to vz=%f (was %f)",
 785                     fCurrentVertex, oldVz));
 786     AliFMDESDRevertexer revertexer;
 787     revertexer.Revertex(fESDObj, fCurrentVertex);
 788   }
 789
 790   if (esd) {
 791     AliFMDDebug(2, ("Writing FMD data to ESD tree"));
 792     esd->SetFMDData(fESDObj);
 793   }
 794
 795   if (!fDiagnostics || !esd) return;
 796   static bool first = true;
 797   // This is most likely NOT the event number you'd like to use. It
 798   // has nothing to do with the 'real' event number.
 799   // - That's OK.  We just use it for the name of the directory -
 800   // nothing else.  Christian
 801   Int_t evno = esd->GetEventNumberInFile();
 802   AliFMDDebug(3, ("Writing diagnostics histograms to FMD.Diag.root/%03d",evno));
 803   TFile f("FMD.Diag.root", (first ? "RECREATE" : "UPDATE"));
 804   first = false;
 805   f.cd();
 806   TDirectory* d = f.mkdir(Form("%03d", evno),
 807                           Form("Diagnostics histograms for event # %d", evno));
 808   d->cd();
 809   if (fDiagStep1) fDiagStep1->Write();
 810   if (fDiagStep2) fDiagStep2->Write();
 811   if (fDiagStep3) fDiagStep3->Write();
 812   if (fDiagStep4) fDiagStep4->Write();
 813   if (fDiagAll)   fDiagAll->Write();
 814   d->Write();
 815   f.Write();
 816   f.Close();
 817
 818   if (fDiagStep1) fDiagStep1->Reset();
 819   if (fDiagStep2) fDiagStep2->Reset();
 820   if (fDiagStep3) fDiagStep3->Reset();
 821   if (fDiagStep4) fDiagStep4->Reset();
 822   if (fDiagAll)   fDiagAll->Reset();
 823 }
 824
 825 //____________________________________________________________________
 826 void
 827 AliFMDReconstructor::FillESD(AliRawReader*, TTree* clusterTree,
 828                              AliESDEvent* esd) const
 829 {
 830   TTree* dummy = 0;
 831   FillESD(dummy, clusterTree, esd);
 832 }
 833 //____________________________________________________________________
 834 //
 835 // EOF
 836 //