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 <AliRunLoader.h>                  // ALIRUNLOADER_H
  38 #include <AliHeader.h>                     // ALIHEADER_H
  39 #include <AliGenEventHeader.h>             // ALIGENEVENTHEADER_H
  40 #include "AliFMDGeometry.h"                // ALIFMDGEOMETRY_H
  41 #include "AliFMDParameters.h"              // ALIFMDPARAMETERS_H
  42 #include "AliFMDDigit.h"                   // ALIFMDDIGIT_H
  43 #include "AliFMDReconstructor.h"           // ALIFMDRECONSTRUCTOR_H
  44 #include "AliFMDRawReader.h"               // ALIFMDRAWREADER_H
  45 #include "AliFMDRecPoint.h"                // ALIFMDMULTNAIIVE_H
  46 #include "AliESD.h"                        // ALIESD_H
  47 #include <AliESDFMD.h>                     // ALIESDFMD_H
  48 class AliRawReader;
  49
  50 //____________________________________________________________________
  51 ClassImp(AliFMDReconstructor)
  52 #if 0
  53   ; // This is here to keep Emacs for indenting the next line
  54 #endif
  55
  56 //____________________________________________________________________
  57 AliFMDReconstructor::AliFMDReconstructor()
  58   : AliReconstructor(),
  59     fMult(0x0),
  60     fNMult(0),
  61     fTreeR(0x0),
  62     fCurrentVertex(0),
  63     fESDObj(0x0),
  64     fESD(0x0)
  65 {
  66   // Make a new FMD reconstructor object - default CTOR.
  67   SetNoiseFactor();
  68   SetAngleCorrect();
  69 }
  70
  71
  72 //____________________________________________________________________
  73 AliFMDReconstructor::AliFMDReconstructor(const AliFMDReconstructor& other)
  74   : AliReconstructor(),
  75     fMult(other.fMult),
  76     fNMult(other.fNMult),
  77     fTreeR(other.fTreeR),
  78     fCurrentVertex(other.fCurrentVertex),
  79     fESDObj(other.fESDObj),
  80     fESD(other.fESD),
  81     fNoiseFactor(other.fNoiseFactor),
  82     fAngleCorrect(other.fAngleCorrect)
  83 {
  84   // Copy constructor
  85 }
  86
  87
  88 //____________________________________________________________________
  89 AliFMDReconstructor&
  90 AliFMDReconstructor::operator=(const AliFMDReconstructor& other)
  91 {
  92   // Assignment operator
  93   fMult          = other.fMult;
  94   fNMult         = other.fNMult;
  95   fTreeR         = other.fTreeR;
  96   fCurrentVertex = other.fCurrentVertex;
  97   fESDObj        = other.fESDObj;
  98   fESD           = other.fESD;
  99   fNoiseFactor   = other.fNoiseFactor;
 100   fAngleCorrect  = other.fAngleCorrect;
 101
 102   return *this;
 103 }
 104
 105 //____________________________________________________________________
 106 AliFMDReconstructor::~AliFMDReconstructor()
 107 {
 108   // Destructor
 109   if (fMult)   fMult->Delete();
 110   if (fMult)   delete fMult;
 111   if (fESDObj) delete fESDObj;
 112 }
 113
 114 //____________________________________________________________________
 115 void
 116 AliFMDReconstructor::Init(AliRunLoader* runLoader)
 117 {
 118   // Initialize the reconstructor
 119   AliDebug(1, Form("Init called with runloader 0x%x", runLoader));
 120   // Initialize the geometry
 121   AliFMDGeometry* geom = AliFMDGeometry::Instance();
 122   geom->Init();
 123   geom->InitTransformations();
 124
 125   // Initialize the parameters
 126   AliFMDParameters* param = AliFMDParameters::Instance();
 127   param->Init();
 128
 129   // Current vertex position
 130   fCurrentVertex = 0;
 131   // Create array of reconstructed strip multiplicities
 132   fMult = new TClonesArray("AliFMDRecPoint", 51200);
 133   // Create ESD output object
 134   fESDObj = new AliESDFMD;
 135
 136   // Check that we have a run loader
 137   if (!runLoader) {
 138     Warning("Init", "No run loader");
 139     return;
 140   }
 141
 142   // Check if we can get the header tree
 143   AliHeader* header = runLoader->GetHeader();
 144   if (!header) {
 145     Warning("Init", "No header");
 146     return;
 147   }
 148
 149   // Check if we can get a simulation header
 150   AliGenEventHeader* eventHeader = header->GenEventHeader();
 151   if (eventHeader) {
 152     TArrayF vtx;
 153     eventHeader->PrimaryVertex(vtx);
 154     fCurrentVertex = vtx[2];
 155     AliDebug(1, Form("Primary vertex Z coordinate for event # %d/%d is %f",
 156                      header->GetRun(), header->GetEvent(), fCurrentVertex));
 157     Warning("Init", "no generator event header");
 158   }
 159   else {
 160     Warning("Init", "No generator event header - "
 161             "perhaps we get the vertex from ESD?");
 162   }
 163 }
 164
 165 //____________________________________________________________________
 166 void
 167 AliFMDReconstructor::ConvertDigits(AliRawReader* reader,
 168                                    TTree* digitsTree) const
 169 {
 170   // Convert Raw digits to AliFMDDigit's in a tree
 171   AliDebug(1, "Reading raw data into digits tree");
 172   AliFMDRawReader rawRead(reader, digitsTree);
 173   // rawRead.SetSampleRate(fFMD->GetSampleRate());
 174   rawRead.Exec();
 175 }
 176
 177 //____________________________________________________________________
 178 void
 179 AliFMDReconstructor::Reconstruct(TTree* digitsTree,
 180                                  TTree* clusterTree) const
 181 {
 182   // Reconstruct event from digits in tree
 183   // Get the FMD branch holding the digits.
 184   // FIXME: The vertex may not be known yet, so we may have to move
 185   // some of this to FillESD.
 186   AliDebug(1, "Reconstructing from digits in a tree");
 187 #if 1
 188   if (fESD) {
 189     const AliESDVertex* vertex = fESD->GetVertex();
 190     if (vertex) {
 191       AliDebug(1, Form("Got vertex from ESD: %f", vertex->GetZv()));
 192       fCurrentVertex = vertex->GetZv();
 193     }
 194   }
 195 #endif
 196   TBranch *digitBranch = digitsTree->GetBranch("FMD");
 197   if (!digitBranch) {
 198     Error("Exec", "No digit branch for the FMD found");
 199     return;
 200   }
 201   TClonesArray* digits = new TClonesArray("AliFMDDigit");
 202   digitBranch->SetAddress(&digits);
 203
 204   if (fMult)   fMult->Clear();
 205   if (fESDObj) fESDObj->Clear();
 206
 207   fNMult = 0;
 208   fTreeR = clusterTree;
 209   fTreeR->Branch("FMD", &fMult);
 210
 211   AliDebug(5, "Getting entry 0 from digit branch");
 212   digitBranch->GetEntry(0);
 213
 214   AliDebug(5, "Processing digits");
 215   ProcessDigits(digits);
 216
 217   Int_t written = clusterTree->Fill();
 218   AliDebug(10, Form("Filled %d bytes into cluster tree", written));
 219   digits->Delete();
 220   delete digits;
 221 }
 222
 223
 224 //____________________________________________________________________
 225 void
 226 AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
 227 {
 228   // For each digit, find the pseudo rapdity, azimuthal angle, and
 229   // number of corrected ADC counts, and pass it on to the algorithms
 230   // used.
 231   Int_t nDigits = digits->GetEntries();
 232   AliDebug(1, Form("Got %d digits", nDigits));
 233   fESDObj->SetNoiseFactor(fNoiseFactor);
 234   fESDObj->SetAngleCorrected(fAngleCorrect);
 235   for (Int_t i = 0; i < nDigits; i++) {
 236     AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
 237     AliFMDParameters* param  = AliFMDParameters::Instance();
 238     // Check that the strip is not marked as dead
 239     if (param->IsDead(digit->Detector(), digit->Ring(),
 240                       digit->Sector(), digit->Strip())) {
 241       AliDebug(10, Form("FMD%d%c[%2d,%3d] is dead", digit->Detector(),
 242                         digit->Ring(), digit->Sector(), digit->Strip()));
 243       continue;
 244     }
 245
 246     // digit->Print();
 247     // Get eta and phi
 248     Float_t eta, phi;
 249     PhysicalCoordinates(digit, eta, phi);
 250
 251     // Substract pedestal.
 252     UShort_t counts   = SubtractPedestal(digit);
 253
 254     // Gain match digits.
 255     Double_t edep     = Adc2Energy(digit, eta, counts);
 256
 257     // Make rough multiplicity
 258     Double_t mult     = Energy2Multiplicity(digit, edep);
 259
 260     AliDebug(10, Form("FMD%d%c[%2d,%3d]: "
 261                       "ADC: %d, Counts: %d, Energy: %f, Mult: %f",
 262                       digit->Detector(), digit->Ring(), digit->Sector(),
 263                       digit->Strip(), digit->Counts(), counts, edep, mult));
 264
 265     // Create a `RecPoint' on the output branch.
 266     AliFMDRecPoint* m =
 267       new ((*fMult)[fNMult]) AliFMDRecPoint(digit->Detector(),
 268                                             digit->Ring(),
 269                                             digit->Sector(),
 270                                             digit->Strip(),
 271                                             eta, phi,
 272                                             edep, mult);
 273     (void)m; // Suppress warnings about unused variables.
 274     fNMult++;
 275
 276     fESDObj->SetMultiplicity(digit->Detector(), digit->Ring(),
 277                              digit->Sector(),  digit->Strip(), mult);
 278     fESDObj->SetEta(digit->Detector(), digit->Ring(),
 279                     digit->Sector(),  digit->Strip(), eta);
 280   }
 281 }
 282
 283 //____________________________________________________________________
 284 UShort_t
 285 AliFMDReconstructor::SubtractPedestal(AliFMDDigit* digit) const
 286 {
 287   // Member function to subtract the pedestal from a digit
 288   // This implementation does nothing, but a derived class could over
 289   // load this to subtract a pedestal that was given in a database or
 290   // something like that.
 291
 292   Int_t             counts = 0;
 293   AliFMDParameters* param  = AliFMDParameters::Instance();
 294   Float_t           ped    = param->GetPedestal(digit->Detector(),
 295                                                 digit->Ring(),
 296                                                 digit->Sector(),
 297                                                 digit->Strip());
 298   Float_t           noise  = param->GetPedestalWidth(digit->Detector(),
 299                                                      digit->Ring(),
 300                                                      digit->Sector(),
 301                                                      digit->Strip());
 302   AliDebug(15, Form("Subtracting pedestal %f from signal %d",
 303                    ped, digit->Counts()));
 304   if (digit->Count3() > 0)      counts = digit->Count3();
 305   else if (digit->Count2() > 0) counts = digit->Count2();
 306   else                          counts = digit->Count1();
 307   counts = TMath::Max(Int_t(counts - ped), 0);
 308   if (counts < noise * fNoiseFactor) counts = 0;
 309   if (counts > 0) AliDebug(15, "Got a hit strip");
 310
 311   return  UShort_t(counts);
 312 }
 313
 314 //____________________________________________________________________
 315 Float_t
 316 AliFMDReconstructor::Adc2Energy(AliFMDDigit* digit,
 317                                 Float_t      eta,
 318                                 UShort_t     count) const
 319 {
 320   // Converts number of ADC counts to energy deposited.
 321   // Note, that this member function can be overloaded by derived
 322   // classes to do strip-specific look-ups in databases or the like,
 323   // to find the proper gain for a strip.
 324   //
 325   // In this simple version, we calculate the energy deposited as
 326   //
 327   //    EnergyDeposited = cos(theta) * gain * count
 328   //
 329   // where
 330   //
 331   //           Pre_amp_MIP_Range
 332   //    gain = ----------------- * Energy_deposited_per_MIP
 333   //           ADC_channel_size
 334   //
 335   // is constant and the same for all strips.
 336
 337   if (count <= 0) return 0;
 338   AliFMDParameters* param = AliFMDParameters::Instance();
 339   Float_t           gain  = param->GetPulseGain(digit->Detector(),
 340                                                 digit->Ring(),
 341                                                 digit->Sector(),
 342                                                 digit->Strip());
 343   AliDebug(15, Form("Converting counts %d to energy via factor %f",
 344                     count, gain));
 345
 346   Double_t edep  = count * gain;
 347   if (fAngleCorrect) {
 348     Double_t theta =  2 * TMath::ATan(TMath::Exp(-eta));
 349     edep           *= TMath::Abs(TMath::Cos(theta));
 350   }
 351   return edep;
 352 }
 353
 354 //____________________________________________________________________
 355 Float_t
 356 AliFMDReconstructor::Energy2Multiplicity(AliFMDDigit* /* digit */,
 357                                          Float_t      edep) const
 358 {
 359   // Converts an energy signal to number of particles.
 360   // Note, that this member function can be overloaded by derived
 361   // classes to do strip-specific look-ups in databases or the like,
 362   // to find the proper gain for a strip.
 363   //
 364   // In this simple version, we calculate the multiplicity as
 365   //
 366   //   multiplicity = Energy_deposited / Energy_deposited_per_MIP
 367   //
 368   // where
 369   //
 370   //   Energy_deposited_per_MIP = 1.664 * SI_density * SI_thickness
 371   //
 372   // is constant and the same for all strips
 373   AliFMDParameters* param   = AliFMDParameters::Instance();
 374   Double_t          edepMIP = param->GetEdepMip();
 375   Float_t           mult    = edep / edepMIP;
 376   if (edep > 0)
 377     AliDebug(15, Form("Translating energy %f to multiplicity via "
 378                      "divider %f->%f", edep, edepMIP, mult));
 379   return mult;
 380 }
 381
 382 //____________________________________________________________________
 383 void
 384 AliFMDReconstructor::PhysicalCoordinates(AliFMDDigit* digit,
 385                                          Float_t& eta,
 386                                          Float_t& phi) const
 387 {
 388   // Get the eta and phi of a digit
 389   //
 390   // Get geometry.
 391   AliFMDGeometry* geom = AliFMDGeometry::Instance();
 392   Double_t x, y, z, r, theta;
 393   geom->Detector2XYZ(digit->Detector(), digit->Ring(), digit->Sector(),
 394                     digit->Strip(), x, y, z);
 395   // Correct for vertex offset.
 396   z     += fCurrentVertex;
 397   phi   =  TMath::ATan2(y, x);
 398   r     =  TMath::Sqrt(y * y + x * x);
 399   theta =  TMath::ATan2(r, z);
 400   eta   = -TMath::Log(TMath::Tan(theta / 2));
 401 }
 402
 403
 404
 405 //____________________________________________________________________
 406 void
 407 AliFMDReconstructor::FillESD(TTree*  /* digitsTree */,
 408                              TTree*  /* clusterTree */,
 409                              AliESD* esd) const
 410 {
 411   // nothing to be done
 412   // FIXME: The vertex may not be known when Reconstruct is executed,
 413   // so we may have to move some of that member function here.
 414   AliDebug(1, Form("Calling FillESD with two trees and one ESD"));
 415   // fESDObj->Print();
 416
 417   if (esd) {
 418     AliDebug(1, Form("Writing FMD data to ESD tree"));
 419     esd->SetFMDData(fESDObj);
 420   }
 421 }
 422
 423
 424 //____________________________________________________________________
 425 void
 426 AliFMDReconstructor::Reconstruct(AliRawReader*,TTree*) const
 427 {
 428   // Cannot be used.  See member function with same name but with 2
 429   // TTree arguments.   Make sure you do local reconstrucion
 430   AliDebug(1, Form("Calling FillESD with loader and tree"));
 431   AliError("MayNotUse");
 432 }
 433 //____________________________________________________________________
 434 void
 435 AliFMDReconstructor::Reconstruct(AliRunLoader*) const
 436 {
 437   // Cannot be used.  See member function with same name but with 2
 438   // TTree arguments.   Make sure you do local reconstrucion
 439   AliDebug(1, Form("Calling FillESD with loader"));
 440   AliError("MayNotUse");
 441 }
 442 //____________________________________________________________________
 443 void
 444 AliFMDReconstructor::Reconstruct(AliRunLoader*, AliRawReader*) const
 445 {
 446   // Cannot be used.  See member function with same name but with 2
 447   // TTree arguments.   Make sure you do local reconstrucion
 448   AliDebug(1, Form("Calling FillESD with loader and raw reader"));
 449   AliError("MayNotUse");
 450 }
 451 //____________________________________________________________________
 452 void
 453 AliFMDReconstructor::FillESD(AliRawReader*,TTree*,AliESD*) const
 454 {
 455   // Cannot be used.  See member function with same name but with 2
 456   // TTree arguments.   Make sure you do local reconstrucion
 457   AliDebug(1, Form("Calling FillESD with raw reader, tree, and ESD"));
 458   AliError("MayNotUse");
 459 }
 460 //____________________________________________________________________
 461 void
 462 AliFMDReconstructor::FillESD(AliRunLoader*,AliESD*) const
 463 {
 464   // Cannot be used.  See member function with same name but with 2
 465   // TTree arguments.   Make sure you do local reconstrucion
 466   AliDebug(1, Form("Calling FillESD with loader and ESD"));
 467   AliError("MayNotUse");
 468 }
 469 //____________________________________________________________________
 470 void
 471 AliFMDReconstructor::FillESD(AliRunLoader*,AliRawReader*,AliESD*) const
 472 {
 473   // Cannot be used.  See member function with same name but with 2
 474   // TTree arguments.   Make sure you do local reconstrucion
 475   AliDebug(1, Form("Calling FillESD with loader, raw reader, and ESD"));
 476   AliError("MayNotUse");
 477 }
 478
 479 //____________________________________________________________________
 480 //
 481 // EOF
 482 //