FMD/AliFMDReconstructor.cxx

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