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