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