FMD/AliFMDReconstructor.cxx

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