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   AliFMDDebug(1, ("Reconstructing from raw reader"));
 305   AliFMDRawReader rawReader(reader, 0);
 306
 307   UShort_t det, sec, str, fac;
 308   Short_t  adc, oldDet = -1;
 309   Bool_t   zs;
 310   Char_t   rng;
 311
 312   UseRecoParam(kTRUE);
 313   while (rawReader.NextSignal(det, rng, sec, str, adc, zs, fac)) {
 314     if (det != oldDet) {
 315       fZS[det-1]       = zs;
 316       fZSFactor[det-1] = fac;
 317       oldDet           = det;
 318     }
 319     ProcessSignal(det, rng, sec, str, adc);
 320   }
 321   UseRecoParam(kFALSE);
 322
 323 }
 324
 325 //____________________________________________________________________
 326 void
 327 AliFMDReconstructor::Digitize(AliRawReader* reader, TClonesArray* sdigits) const
 328 {
 329   // Reconstruct directly from raw data (no intermediate output on
 330   // digit tree or rec point tree).
 331   //
 332   // Parameters:
 333   //   reader   Raw event reader
 334   //   ctree    Not used.
 335   AliFMDRawReader rawReader(reader, 0);
 336
 337   UShort_t det, sec, str, sam, rat, fac;
 338   Short_t  adc, oldDet = -1;
 339   Bool_t   zs;
 340   Char_t   rng;
 341
 342   UseRecoParam(kTRUE);
 343   while (rawReader.NextSample(det, rng, sec, str, sam, rat, adc, zs, fac)) {
 344     if (!rawReader.SelectSample(sam, rat)) continue;
 345     if (det != oldDet) {
 346       fZS[det-1]       = zs;
 347       fZSFactor[det-1] = fac;
 348       oldDet           = det;
 349     }
 350     DigitizeSignal(sdigits, det, rng, sec, str, sam, adc);
 351   }
 352   UseRecoParam(kFALSE);
 353 }
 354
 355 //____________________________________________________________________
 356 void
 357 AliFMDReconstructor::Reconstruct(TTree* digitsTree,
 358                                  TTree* clusterTree) const
 359 {
 360   // Reconstruct event from digits in tree
 361   // Get the FMD branch holding the digits.
 362   // FIXME: The vertex may not be known yet, so we may have to move
 363   // some of this to FillESD.
 364   //
 365   // Parameters:
 366   //   digitsTree       Pointer to a tree containing digits
 367   //   clusterTree      Pointer to output tree
 368   //
 369   AliFMDDebug(1, ("Reconstructing from digits in a tree"));
 370   GetVertex(fESD);
 371
 372   TBranch *digitBranch = digitsTree->GetBranch("FMD");
 373   if (!digitBranch) {
 374     Error("Exec", "No digit branch for the FMD found");
 375     return;
 376   }
 377   TClonesArray* digits = new TClonesArray("AliFMDDigit");
 378   digitBranch->SetAddress(&digits);
 379
 380   if (fMult)   fMult->Clear();
 381   if (fESDObj) fESDObj->Clear();
 382
 383   fNMult = 0;
 384   fTreeR = clusterTree;
 385   fTreeR->Branch("FMD", &fMult);
 386
 387   AliFMDDebug(5, ("Getting entry 0 from digit branch"));
 388   digitBranch->GetEntry(0);
 389
 390   AliFMDDebug(5, ("Processing digits"));
 391   UseRecoParam(kTRUE);
 392   ProcessDigits(digits);
 393   UseRecoParam(kFALSE);
 394
 395   Int_t written = clusterTree->Fill();
 396   AliFMDDebug(10, ("Filled %d bytes into cluster tree", written));
 397   digits->Delete();
 398   delete digits;
 399 }
 400
 401
 402 //____________________________________________________________________
 403 void
 404 AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
 405 {
 406   // For each digit, find the pseudo rapdity, azimuthal angle, and
 407   // number of corrected ADC counts, and pass it on to the algorithms
 408   // used.
 409   //
 410   // Parameters:
 411   //    digits  Array of digits
 412   //
 413   Int_t nDigits = digits->GetEntries();
 414   AliFMDDebug(2, ("Got %d digits", nDigits));
 415   fESDObj->SetNoiseFactor(fNoiseFactor);
 416   fESDObj->SetAngleCorrected(fAngleCorrect);
 417   for (Int_t i = 0; i < nDigits; i++) {
 418     AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
 419     if (!digit) continue;
 420     ProcessDigit(digit);
 421   }
 422 }
 423
 424 //____________________________________________________________________
 425 void
 426 AliFMDReconstructor::ProcessDigit(AliFMDDigit* digit) const
 427 {
 428   UShort_t det = digit->Detector();
 429   Char_t   rng = digit->Ring();
 430   UShort_t sec = digit->Sector();
 431   UShort_t str = digit->Strip();
 432   Short_t  adc = digit->Counts();
 433
 434   ProcessSignal(det, rng, sec, str, adc);
 435 }
 436
 437 //____________________________________________________________________
 438 void
 439 AliFMDReconstructor::ProcessSignal(UShort_t det,
 440                                    Char_t   rng,
 441                                    UShort_t sec,
 442                                    UShort_t str,
 443                                    Short_t  adc) const
 444 {
 445   // Process the signal from a single strip
 446   //
 447   // Parameters:
 448   //    det     Detector ID
 449   //    rng     Ring ID
 450   //    sec     Sector ID
 451   //    rng     Strip ID
 452   //    adc     ADC counts
 453   //
 454   AliFMDParameters* param  = AliFMDParameters::Instance();
 455   // Check that the strip is not marked as dead
 456   if (param->IsDead(det, rng, sec, str)) {
 457     AliFMDDebug(3, ("FMD%d%c[%2d,%3d] is dead", det, rng, sec, str));
 458     return;
 459   }
 460
 461   // digit->Print();
 462   // Get eta and phi
 463   Float_t eta, phi;
 464   PhysicalCoordinates(det, rng, sec, str, eta, phi);
 465
 466   // Substract pedestal.
 467   UShort_t counts   = SubtractPedestal(det, rng, sec, str, adc);
 468   if(counts == USHRT_MAX) return;
 469
 470     // Gain match digits.
 471   Double_t edep     = Adc2Energy(det, rng, sec, str, eta, counts);
 472   // Get rid of nonsense energy
 473   if(edep < 0)  return;
 474
 475   // Make rough multiplicity
 476   Double_t mult     = Energy2Multiplicity(det, rng, sec, str, edep);
 477   // Get rid of nonsense mult
 478   if (mult < 0)  return;
 479   AliFMDDebug(10, ("FMD%d%c[%2d,%3d]: "
 480                     "ADC: %d, Counts: %d, Energy: %f, Mult: %f",
 481                   det, rng, sec, str, adc, counts, edep, mult));
 482
 483   // Create a `RecPoint' on the output branch.
 484   if (fMult) {
 485     AliFMDRecPoint* m =
 486       new ((*fMult)[fNMult]) AliFMDRecPoint(det, rng, sec, str,
 487                                             eta, phi, edep, mult);
 488     (void)m; // Suppress warnings about unused variables.
 489     fNMult++;
 490   }
 491
 492   fESDObj->SetMultiplicity(det, rng, sec, str, mult);
 493   fESDObj->SetEta(det, rng, sec, str, eta);
 494
 495   if (fDiagAll) fDiagAll->Fill(adc, mult);
 496
 497 }
 498
 499 //____________________________________________________________________
 500 void
 501 AliFMDReconstructor::DigitizeSignal(TClonesArray* sdigits,
 502                                     UShort_t det,
 503                                     Char_t   rng,
 504                                     UShort_t sec,
 505                                     UShort_t str,
 506                                     UShort_t /* sam */,
 507                                     Short_t  adc) const
 508 {
 509   // Process the signal from a single strip
 510   //
 511   // Parameters:
 512   //    det     Detector ID
 513   //    rng     Ring ID
 514   //    sec     Sector ID
 515   //    rng     Strip ID
 516   //    adc     ADC counts
 517   //
 518   AliFMDParameters* param  = AliFMDParameters::Instance();
 519   // Check that the strip is not marked as dead
 520   if (param->IsDead(det, rng, sec, str)) {
 521     AliFMDDebug(10, ("FMD%d%c[%2d,%3d] is dead", det, rng, sec, str));
 522     return;
 523   }
 524
 525   // Substract pedestal.
 526   UShort_t counts   = SubtractPedestal(det, rng, sec, str, adc);
 527   if(counts == USHRT_MAX || counts == 0) return;
 528
 529     // Gain match digits.
 530   Double_t edep     = Adc2Energy(det, rng, sec, str, counts);
 531   // Get rid of nonsense energy
 532   if(edep < 0)  return;
 533
 534   Int_t n = sdigits->GetEntriesFast();
 535   // AliFMDSDigit* sdigit =
 536   new ((*sdigits)[n])
 537     AliFMDSDigit(det, rng, sec, str, edep, counts, counts, counts, counts);
 538   // sdigit->SetCount(sam, counts);
 539 }
 540
 541 //____________________________________________________________________
 542 UShort_t
 543 AliFMDReconstructor::SubtractPedestal(UShort_t det,
 544                                       Char_t   rng,
 545                                       UShort_t sec,
 546                                       UShort_t str,
 547                                       UShort_t adc,
 548                                       Float_t  noiseFactor,
 549                                       Bool_t   zsEnabled,
 550                                       UShort_t zsNoiseFactor) const
 551 {
 552   AliFMDParameters* param  = AliFMDParameters::Instance();
 553   Float_t           ped    = (zsEnabled ? 0 :
 554                                 param->GetPedestal(det, rng, sec, str));
 555   Float_t           noise  = param->GetPedestalWidth(det, rng, sec, str);
 556   if(ped < 0 || noise < 0) {
 557     AliWarningClass(Form("Invalid pedestal (%f) or noise (%f) "
 558                          "for FMD%d%c[%02d,%03d]",
 559                     ped, noise, det, rng, sec, str));
 560     return USHRT_MAX;
 561   }
 562   AliDebugClass(15, Form("Subtracting pedestal for FMD%d%c[%2d,%3d]=%4d "
 563                          "(%s w/factor %d, noise factor %f, "
 564                          "pedestal %8.2f+/-%8.2f)",
 565                          det, rng, sec, str, adc,
 566                          (zsEnabled ? "zs'ed" : "straight"),
 567                          zsNoiseFactor, noiseFactor, ped, noise));
 568
 569   Int_t counts = adc + Int_t(zsEnabled ? zsNoiseFactor * noise : - ped);
 570   counts =  TMath::Max(Int_t(counts), 0);
 571   // Calculate the noise factor for suppressing remenants of the noise
 572   // peak.  If we have done on-line zero suppression, we only check
 573   // for noise signals that are larger than the suppressed noise.  If
 574   // the noise factor used on line is larger than the factor used
 575   // here, we do not do this check at all.
 576   //
 577   // For example:
 578   //    Online factor  |  Read factor |  Result
 579   //    ---------------+--------------+-------------------------------
 580   //           2       |      3       | Check if signal > 1 * noise
 581   //           3       |      3       | Check if signal > 0
 582   //           3       |      2       | Check if signal > 0
 583   //
 584   // In this way, we make sure that we do not suppress away too much
 585   // data, and that the read-factor is the most stringent cut.
 586   Float_t nf = TMath::Max(0.F, noiseFactor - (zsEnabled ? zsNoiseFactor : 0));
 587   if (counts < noise * nf) counts = 0;
 588   if (counts > 0) AliDebugClass(15, "Got a hit strip");
 589
 590   return counts;
 591 }
 592
 593
 594 //____________________________________________________________________
 595 UShort_t
 596 AliFMDReconstructor::SubtractPedestal(UShort_t det,
 597                                       Char_t   rng,
 598                                       UShort_t sec,
 599                                       UShort_t str,
 600                                       Short_t  adc) const
 601 {
 602   // Member function to subtract the pedestal from a digit
 603   //
 604   // Parameters:
 605   //    det     Detector ID
 606   //    rng     Ring ID
 607   //    sec     Sector ID
 608   //    rng     Strip ID
 609   //    adc     # of ADC counts
 610   // Return:
 611   //    Pedestal subtracted signal or USHRT_MAX in case of problems
 612   //
 613   UShort_t counts = SubtractPedestal(det, rng, sec, str, adc,
 614                                      fNoiseFactor, fZS[det-1],
 615                                      fZSFactor[det-1]);
 616   if (fDiagStep1) fDiagStep1->Fill(adc, counts);
 617
 618   return counts;
 619 }
 620
 621 //____________________________________________________________________
 622 Float_t
 623 AliFMDReconstructor::Adc2Energy(UShort_t det,
 624                                 Char_t   rng,
 625                                 UShort_t sec,
 626                                 UShort_t str,
 627                                 UShort_t count) const
 628 {
 629   // Converts number of ADC counts to energy deposited.
 630   // Note, that this member function can be overloaded by derived
 631   // classes to do strip-specific look-ups in databases or the like,
 632   // to find the proper gain for a strip.
 633   //
 634   // In the first simple version, we calculate the energy deposited as
 635   //
 636   //    EnergyDeposited = cos(theta) * gain * count
 637   //
 638   // where
 639   //
 640   //           Pre_amp_MIP_Range
 641   //    gain = ----------------- * Energy_deposited_per_MIP
 642   //           ADC_channel_size
 643   //
 644   // is constant and the same for all strips.
 645   //
 646   // For the production we use the conversion measured in the NBI lab.
 647   // The total conversion is then:
 648   //
 649   //    gain = ADC / DAC
 650   //
 651   //                  EdepMip * count
 652   //      => energy = ----------------
 653   //                  gain * DACPerADC
 654   //
 655   // Parameters:
 656   //    det     Detector ID
 657   //    rng     Ring ID
 658   //    sec     Sector ID
 659   //    rng     Strip ID
 660   //    counts  Number of ADC counts over pedestal
 661   // Return
 662   //    The energy deposited in a single strip, or -1 in case of problems
 663   //
 664   if (count <= 0) return 0;
 665   AliFMDParameters* param = AliFMDParameters::Instance();
 666   Float_t           gain  = param->GetPulseGain(det, rng, sec, str);
 667   // 'Tagging' bad gains as bad energy
 668   if (gain < 0) {
 669     AliWarning(Form("Invalid gain (%f) for FMD%d%c[%02d,%03d]",
 670                     gain, det, rng, sec, str));
 671     return -1;
 672   }
 673   AliFMDDebug(5, ("Converting counts %d to energy (factor=%f, DAC2MIP=%f)",
 674                   count, gain,param->GetDACPerMIP()));
 675
 676   Double_t edep  = ((count * param->GetEdepMip())
 677                     / (gain * param->GetDACPerMIP()));
 678   return edep;
 679 }
 680
 681 //____________________________________________________________________
 682 Float_t
 683 AliFMDReconstructor::Adc2Energy(UShort_t det,
 684                                 Char_t   rng,
 685                                 UShort_t sec,
 686                                 UShort_t str,
 687                                 Float_t  eta,
 688                                 UShort_t count) const
 689 {
 690   // Converts number of ADC counts to energy deposited.
 691   // Note, that this member function can be overloaded by derived
 692   // classes to do strip-specific look-ups in databases or the like,
 693   // to find the proper gain for a strip.
 694   //
 695   // In the first simple version, we calculate the energy deposited as
 696   //
 697   //    EnergyDeposited = cos(theta) * gain * count
 698   //
 699   // where
 700   //
 701   //           Pre_amp_MIP_Range
 702   //    gain = ----------------- * Energy_deposited_per_MIP
 703   //           ADC_channel_size
 704   //
 705   // is constant and the same for all strips.
 706   //
 707   // For the production we use the conversion measured in the NBI lab.
 708   // The total conversion is then:
 709   //
 710   //    gain = ADC / DAC
 711   //
 712   //                  EdepMip * count
 713   //      => energy = ----------------
 714   //                  gain * DACPerADC
 715   //
 716   // Parameters:
 717   //    det     Detector ID
 718   //    rng     Ring ID
 719   //    sec     Sector ID
 720   //    rng     Strip ID
 721   //    eta     Psuedo-rapidity
 722   //    counts  Number of ADC counts over pedestal
 723   // Return
 724   //    The energy deposited in a single strip, or -1 in case of problems
 725   //
 726   Double_t edep = Adc2Energy(det, rng, sec, str, count);
 727
 728   if (fDiagStep2) fDiagStep2->Fill(count, edep);
 729   if (fAngleCorrect) {
 730     Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
 731     Double_t corr  = TMath::Abs(TMath::Cos(theta));
 732     Double_t cedep = corr * edep;
 733     AliFMDDebug(10, ("correcting for path %f * %f = %f (eta=%f, theta=%f)",
 734                       edep, corr, cedep, eta, theta));
 735     if (fDiagStep3) fDiagStep3->Fill(edep, cedep);
 736     edep = cedep;
 737   }
 738   return edep;
 739 }
 740
 741 //____________________________________________________________________
 742 Float_t
 743 AliFMDReconstructor::Energy2Multiplicity(UShort_t /*det*/,
 744                                          Char_t   /*rng*/,
 745                                          UShort_t /*sec*/,
 746                                          UShort_t /*str*/,
 747                                          Float_t  edep) const
 748 {
 749   // Converts an energy signal to number of particles.
 750   // Note, that this member function can be overloaded by derived
 751   // classes to do strip-specific look-ups in databases or the like,
 752   // to find the proper gain for a strip.
 753   //
 754   // In this simple version, we calculate the multiplicity as
 755   //
 756   //   multiplicity = Energy_deposited / Energy_deposited_per_MIP
 757   //
 758   // where
 759   //
 760   //   Energy_deposited_per_MIP = 1.664 * SI_density * SI_thickness
 761   //
 762   // is constant and the same for all strips
 763   //
 764   // Parameters:
 765   //    det     Detector ID
 766   //    rng     Ring ID
 767   //    sec     Sector ID
 768   //    rng     Strip ID
 769   //    edep    Energy deposited in a single strip
 770   // Return
 771   //    The "bare" multiplicity corresponding to the energy deposited
 772   AliFMDParameters* param   = AliFMDParameters::Instance();
 773   Double_t          edepMIP = param->GetEdepMip();
 774   Float_t           mult    = edep / edepMIP;
 775 #if 0
 776   if (edep > 0)
 777     AliFMDDebug(15, ("Translating energy %f to multiplicity via "
 778                      "divider %f->%f", edep, edepMIP, mult));
 779 #endif
 780   if (fDiagStep4) fDiagStep4->Fill(edep, mult);
 781   return mult;
 782 }
 783
 784 //____________________________________________________________________
 785 void
 786 AliFMDReconstructor::PhysicalCoordinates(UShort_t det,
 787                                          Char_t   rng,
 788                                          UShort_t sec,
 789                                          UShort_t str,
 790                                          Float_t& eta,
 791                                          Float_t& phi) const
 792 {
 793   // Get the eta and phi of a digit
 794   //
 795   // Parameters:
 796   //    det     Detector ID
 797   //    rng     Ring ID
 798   //    sec     Sector ID
 799   //    rng     Strip ID
 800   //    eta     On return, contains the psuedo-rapidity of the strip
 801   //    phi     On return, contains the azimuthal angle of the strip
 802   //
 803   AliFMDGeometry* geom = AliFMDGeometry::Instance();
 804   Double_t x, y, z, r, theta;
 805   geom->Detector2XYZ(det, rng, sec, str, x, y, z);
 806   // Correct for vertex offset.
 807   z     -= fCurrentVertex;
 808   phi   =  TMath::ATan2(y, x);
 809   r     =  TMath::Sqrt(y * y + x * x);
 810   theta =  TMath::ATan2(r, z);
 811   eta   = -TMath::Log(TMath::Tan(theta / 2));
 812 }
 813
 814
 815
 816 //____________________________________________________________________
 817 void
 818 AliFMDReconstructor::FillESD(TTree*  /* digitsTree */,
 819                              TTree*  /* clusterTree */,
 820                              AliESDEvent* esd) const
 821 {
 822   // nothing to be done
 823   // FIXME: The vertex may not be known when Reconstruct is executed,
 824   // so we may have to move some of that member function here.
 825   AliFMDDebug(2, ("Calling FillESD with two trees and one ESD"));
 826   // fESDObj->Print();
 827
 828   Double_t oldVz = fCurrentVertex;
 829   GetVertex(esd);
 830   if (fVertexType != kNoVertex) {
 831     AliFMDDebug(2, ("Revertexing the ESD data to vz=%f (was %f)",
 832                     fCurrentVertex, oldVz));
 833     AliFMDESDRevertexer revertexer;
 834     revertexer.Revertex(fESDObj, fCurrentVertex);
 835   }
 836
 837   if (esd) {
 838     AliFMDDebug(2, ("Writing FMD data to ESD tree"));
 839     esd->SetFMDData(fESDObj);
 840   }
 841
 842   if (!fDiagnostics || !esd) return;
 843   static bool first = true;
 844   // This is most likely NOT the event number you'd like to use. It
 845   // has nothing to do with the 'real' event number.
 846   // - That's OK.  We just use it for the name of the directory -
 847   // nothing else.  Christian
 848   Int_t evno = esd->GetEventNumberInFile();
 849   AliFMDDebug(3, ("Writing diagnostics histograms to FMD.Diag.root/%03d",evno));
 850   TFile f("FMD.Diag.root", (first ? "RECREATE" : "UPDATE"));
 851   first = false;
 852   f.cd();
 853   TDirectory* d = f.mkdir(Form("%03d", evno),
 854                           Form("Diagnostics histograms for event # %d", evno));
 855   d->cd();
 856   if (fDiagStep1) fDiagStep1->Write();
 857   if (fDiagStep2) fDiagStep2->Write();
 858   if (fDiagStep3) fDiagStep3->Write();
 859   if (fDiagStep4) fDiagStep4->Write();
 860   if (fDiagAll)   fDiagAll->Write();
 861   d->Write();
 862   f.Write();
 863   f.Close();
 864
 865   if (fDiagStep1) fDiagStep1->Reset();
 866   if (fDiagStep2) fDiagStep2->Reset();
 867   if (fDiagStep3) fDiagStep3->Reset();
 868   if (fDiagStep4) fDiagStep4->Reset();
 869   if (fDiagAll)   fDiagAll->Reset();
 870 }
 871
 872 //____________________________________________________________________
 873 void
 874 AliFMDReconstructor::FillESD(AliRawReader*, TTree* clusterTree,
 875                              AliESDEvent* esd) const
 876 {
 877   TTree* dummy = 0;
 878   FillESD(dummy, clusterTree, esd);
 879 }
 880 //_____________________________________________________________________
 881 Bool_t AliFMDReconstructor::GetFMDAsideBit(AliESDFMD* fmd) {
 882
 883   AliFMDParameters* pars = AliFMDParameters::Instance();
 884   Float_t totalMult = 0;
 885   for(UShort_t det=1;det<=2;det++) {
 886     Int_t nRings = (det == 1 ? 1 : 2);
 887     for (UShort_t ir = 0; ir < nRings; ir++) {
 888       Char_t   ring = (ir == 0 ? 'I' : 'O');
 889       UShort_t nsec = (ir == 0 ? 20  : 40);
 890       UShort_t nstr = (ir == 0 ? 512 : 256);
 891       for(UShort_t sec =0; sec < nsec;  sec++)  {
 892         for(UShort_t strip = 0; strip < nstr; strip++) {
 893           Float_t mult = fmd->Multiplicity(det,ring,sec,strip);
 894           if(mult == AliESDFMD::kInvalidMult) continue;
 895
 896           if(mult > pars->GetOfflineTriggerLowCut())
 897             totalMult = totalMult + mult;
 898           else
 899             {
 900               if( totalMult > pars->GetOfflineTriggerHitCut()) {
 901                 return kTRUE;
 902               }
 903               else totalMult = 0 ;
 904             }
 905         }
 906       }
 907     }
 908   }
 909   return kFALSE;
 910
 911 }
 912 //_____________________________________________________________________
 913 Bool_t AliFMDReconstructor::GetFMDCsideBit(AliESDFMD* fmd) {
 914
 915   AliFMDParameters* pars = AliFMDParameters::Instance();
 916   Float_t totalMult = 0;
 917   UShort_t det = 3;
 918   Int_t nRings = 2;
 919   for (UShort_t ir = 0; ir < nRings; ir++) {
 920     Char_t   ring = (ir == 0 ? 'I' : 'O');
 921     UShort_t nsec = (ir == 0 ? 20  : 40);
 922     UShort_t nstr = (ir == 0 ? 512 : 256);
 923     for(UShort_t sec =0; sec < nsec;  sec++)  {
 924       for(UShort_t strip = 0; strip < nstr; strip++) {
 925         Float_t mult = fmd->Multiplicity(det,ring,sec,strip);
 926         if(mult == AliESDFMD::kInvalidMult) continue;
 927
 928         if(mult > pars->GetOfflineTriggerLowCut())
 929           totalMult = totalMult + mult;
 930         else
 931           {
 932             if( totalMult > pars->GetOfflineTriggerHitCut()) {
 933               return kTRUE;
 934             }
 935             else totalMult = 0 ;
 936           }
 937       }
 938     }
 939   }
 940
 941   return kFALSE;
 942 }
 943 //____________________________________________________________________
 944 //
 945 // EOF
 946 //