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
  16 /* $Id$ */
  17
  18 //____________________________________________________________________
  19 //
  20 // This is a class that constructs AliFMDMult (reconstructed
  21 // multiplicity) from of Digits
  22 //
  23 // This class reads either digits from a TClonesArray or raw data from
  24 // a DDL file (or similar), and stores the read ADC counts in an
  25 // internal cache (fAdcs).
  26 //
  27 //-- Authors: Evgeny Karpechev(INR) and Alla Maevsksia
  28 //  Latest changes by Christian Holm Christensen <cholm@nbi.dk>
  29 //
  30 //
  31 //____________________________________________________________________
  32
  33 #include <AliLog.h>                        // ALILOG_H
  34 #include <AliRun.h>                        // ALIRUN_H
  35 #include <AliRunLoader.h>                  // ALIRUNLOADER_H
  36 #include <AliLoader.h>                     // ALILOADER_H
  37 #include <AliHeader.h>                     // ALIHEADER_H
  38 #include <AliRawReader.h>                  // ALIRAWREADER_H
  39 #include <AliGenEventHeader.h>             // ALIGENEVENTHEADER_H
  40 #include "AliFMD.h"                        // ALIFMD_H
  41 #include "AliFMDGeometry.h"                // ALIFMDGEOMETRY_H
  42 #include "AliFMDParameters.h"              // ALIFMDPARAMETERS_H
  43 #include "AliFMDDetector.h"                // ALIFMDDETECTOR_H
  44 #include "AliFMDRing.h"                    // ALIFMDRING_H
  45 #include "AliFMDDigit.h"                   // ALIFMDDIGIT_H
  46 #include "AliFMDReconstructor.h"           // ALIFMDRECONSTRUCTOR_H
  47 #include "AliFMDRawStream.h"               // ALIFMDRAWSTREAM_H
  48 #include "AliFMDRawReader.h"               // ALIFMDRAWREADER_H
  49 #include "AliFMDRecPoint.h"                // ALIFMDMULTNAIIVE_H
  50 #include "AliESD.h"                        // ALIESD_H
  51 #include <AliESDFMD.h>                     // ALIESDFMD_H
  52 #include <TFile.h>
  53
  54 //____________________________________________________________________
  55 ClassImp(AliFMDReconstructor)
  56 #if 0
  57   ; // This is here to keep Emacs for indenting the next line
  58 #endif
  59
  60 //____________________________________________________________________
  61 AliFMDReconstructor::AliFMDReconstructor()
  62   : AliReconstructor(),
  63     fMult(0),
  64     fESDObj(0)
  65 {
  66   // Make a new FMD reconstructor object - default CTOR.
  67 }
  68
  69
  70 //____________________________________________________________________
  71 AliFMDReconstructor::AliFMDReconstructor(const AliFMDReconstructor& other)
  72   : AliReconstructor(),
  73     fMult(other.fMult),
  74     fESDObj(other.fESDObj)
  75 {
  76   // Copy constructor
  77 }
  78
  79
  80 //____________________________________________________________________
  81 AliFMDReconstructor&
  82 AliFMDReconstructor::operator=(const AliFMDReconstructor& other)
  83 {
  84   // Assignment operator
  85   fMult   = other.fMult;
  86   fESDObj = other.fESDObj;
  87   return *this;
  88 }
  89
  90 //____________________________________________________________________
  91 AliFMDReconstructor::~AliFMDReconstructor()
  92 {
  93   // Destructor
  94   if (fMult)   fMult->Delete();
  95   if (fMult)   delete fMult;
  96   if (fESDObj) delete fESDObj;
  97 }
  98
  99 //____________________________________________________________________
 100 void
 101 AliFMDReconstructor::Init(AliRunLoader* runLoader)
 102 {
 103   // Initialize the reconstructor
 104   AliDebug(1, Form("Init called with runloader 0x%x", runLoader));
 105   // Initialize the geometry
 106   AliFMDGeometry* fmd = AliFMDGeometry::Instance();
 107   fmd->Init();
 108   fmd->InitTransformations();
 109
 110   // Current vertex position
 111   fCurrentVertex = 0;
 112   // Create array of reconstructed strip multiplicities
 113   fMult = new TClonesArray("AliFMDRecPoint", 51200);
 114   // Create ESD output object
 115   fESDObj = new AliESDFMD;
 116
 117   // Check that we have a run loader
 118   if (!runLoader) {
 119     Warning("Init", "No run loader");
 120     return;
 121   }
 122
 123   // Check if we can get the header tree
 124   AliHeader* header = runLoader->GetHeader();
 125   if (!header) {
 126     Warning("Init", "No header");
 127     return;
 128   }
 129
 130   // Check if we can get a simulation header
 131   AliGenEventHeader* eventHeader = header->GenEventHeader();
 132   if (eventHeader) {
 133     TArrayF vtx;
 134     eventHeader->PrimaryVertex(vtx);
 135     fCurrentVertex = vtx[2];
 136     AliDebug(1, Form("Primary vertex Z coordinate for event # %d/%d is %f",
 137                      header->GetRun(), header->GetEvent(), fCurrentVertex));
 138     Warning("Init", "no generator event header");
 139   }
 140   else {
 141     Warning("Init", "No generator event header - "
 142             "perhaps we get the vertex from ESD?");
 143   }
 144 }
 145
 146 //____________________________________________________________________
 147 void
 148 AliFMDReconstructor::ConvertDigits(AliRawReader* reader,
 149                                    TTree* digitsTree) const
 150 {
 151   // Convert Raw digits to AliFMDDigit's in a tree
 152   AliDebug(1, "Reading raw data into digits tree");
 153   AliFMDRawReader rawRead(reader, digitsTree);
 154   // rawRead.SetSampleRate(fFMD->GetSampleRate());
 155   rawRead.Exec();
 156 }
 157
 158 //____________________________________________________________________
 159 void
 160 AliFMDReconstructor::Reconstruct(TTree* digitsTree,
 161                                  TTree* clusterTree) const
 162 {
 163   // Reconstruct event from digits in tree
 164   // Get the FMD branch holding the digits.
 165   // FIXME: The vertex may not be known yet, so we may have to move
 166   // some of this to FillESD.
 167   AliDebug(1, "Reconstructing from digits in a tree");
 168 #if 1
 169   if (fESD) {
 170     const AliESDVertex* vertex = fESD->GetVertex();
 171     if (vertex) {
 172       AliDebug(1, Form("Got vertex from ESD: %f", vertex->GetZv()));
 173       fCurrentVertex = vertex->GetZv();
 174     }
 175   }
 176 #endif
 177   TBranch *digitBranch = digitsTree->GetBranch("FMD");
 178   if (!digitBranch) {
 179     Error("Exec", "No digit branch for the FMD found");
 180     return;
 181   }
 182   TClonesArray* digits = new TClonesArray("AliFMDDigit");
 183   digitBranch->SetAddress(&digits);
 184
 185   // TIter next(&fAlgorithms);
 186   // AliFMDMultAlgorithm* algorithm = 0;
 187   // while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) {
 188   //   AliDebug(10, Form("PreEvent called for algorithm %s",
 189   //                     algorithm->GetName()));
 190   //   algorithm->PreEvent(clusterTree, fCurrentVertex);
 191   // }
 192   if (fMult)   fMult->Clear();
 193   if (fESDObj) fESDObj->Clear();
 194
 195   fNMult = 0;
 196   fTreeR = clusterTree;
 197   fTreeR->Branch("FMD", &fMult);
 198
 199   AliDebug(5, "Getting entry 0 from digit branch");
 200   digitBranch->GetEntry(0);
 201
 202   AliDebug(5, "Processing digits");
 203   ProcessDigits(digits);
 204
 205   // next.Reset();
 206   // algorithm = 0;
 207   // while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) {
 208   //   AliDebug(10, Form("PostEvent called for algorithm %s",
 209   //                     algorithm->GetName()));
 210   // algorithm->PostEvent();
 211   // }
 212   Int_t written = clusterTree->Fill();
 213   AliDebug(10, Form("Filled %d bytes into cluster tree", written));
 214   digits->Delete();
 215   delete digits;
 216 }
 217
 218
 219 //____________________________________________________________________
 220 void
 221 AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
 222 {
 223   // For each digit, find the pseudo rapdity, azimuthal angle, and
 224   // number of corrected ADC counts, and pass it on to the algorithms
 225   // used.
 226   Int_t nDigits = digits->GetEntries();
 227   AliDebug(1, Form("Got %d digits", nDigits));
 228   for (Int_t i = 0; i < nDigits; i++) {
 229     AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
 230     AliFMDParameters* param  = AliFMDParameters::Instance();
 231     // Check that the strip is not marked as dead
 232     if (param->IsDead(digit->Detector(), digit->Ring(),
 233                       digit->Sector(), digit->Strip())) continue;
 234
 235     // digit->Print();
 236     // Get eta and phi
 237     Float_t eta, phi;
 238     PhysicalCoordinates(digit, eta, phi);
 239
 240     // Substract pedestal.
 241     UShort_t counts   = SubtractPedestal(digit);
 242
 243     // Gain match digits.
 244     Double_t edep     = Adc2Energy(digit, eta, counts);
 245
 246     // Make rough multiplicity
 247     Double_t mult     = Energy2Multiplicity(digit, edep);
 248
 249     AliDebug(10, Form("FMD%d%c[%2d,%3d]: "
 250                       "ADC: %d, Counts: %d, Energy: %f, Mult: %f",
 251                       digit->Detector(), digit->Ring(), digit->Sector(),
 252                       digit->Strip(), digit->Counts(), counts, edep, mult));
 253
 254     // Create a `RecPoint' on the output branch.
 255     AliFMDRecPoint* m =
 256       new ((*fMult)[fNMult]) AliFMDRecPoint(digit->Detector(),
 257                                             digit->Ring(),
 258                                             digit->Sector(),
 259                                             digit->Strip(),
 260                                             eta, phi,
 261                                             edep, mult);
 262     (void)m; // Suppress warnings about unused variables.
 263     fNMult++;
 264
 265     fESDObj->SetMultiplicity(digit->Detector(), digit->Ring(),
 266                              digit->Sector(),  digit->Strip(), mult);
 267     fESDObj->SetEta(digit->Detector(), digit->Ring(),
 268                     digit->Sector(),  digit->Strip(), eta);
 269   }
 270 }
 271
 272 //____________________________________________________________________
 273 UShort_t
 274 AliFMDReconstructor::SubtractPedestal(AliFMDDigit* digit) const
 275 {
 276   // Member function to subtract the pedestal from a digit
 277   // This implementation does nothing, but a derived class could over
 278   // load this to subtract a pedestal that was given in a database or
 279   // something like that.
 280
 281   Int_t             counts = 0;
 282   AliFMDParameters* param  = AliFMDParameters::Instance();
 283   Float_t           pedM   = param->GetPedestal(digit->Detector(),
 284                                                 digit->Ring(),
 285                                                 digit->Sector(),
 286                                                 digit->Strip());
 287   AliDebug(10, Form("Subtracting pedestal %f from signal %d",
 288                    pedM, digit->Counts()));
 289   if (digit->Count3() > 0)      counts = digit->Count3();
 290   else if (digit->Count2() > 0) counts = digit->Count2();
 291   else                          counts = digit->Count1();
 292   counts = TMath::Max(Int_t(counts - pedM), 0);
 293   if (counts > 0) AliDebug(10, "Got a hit strip");
 294
 295   return  UShort_t(counts);
 296 }
 297
 298 //____________________________________________________________________
 299 Float_t
 300 AliFMDReconstructor::Adc2Energy(AliFMDDigit* digit,
 301                                 Float_t      /* eta */,
 302                                 UShort_t     count) const
 303 {
 304   // Converts number of ADC counts to energy deposited.
 305   // Note, that this member function can be overloaded by derived
 306   // classes to do strip-specific look-ups in databases or the like,
 307   // to find the proper gain for a strip.
 308   //
 309   // In this simple version, we calculate the energy deposited as
 310   //
 311   //    EnergyDeposited = cos(theta) * gain * count
 312   //
 313   // where
 314   //
 315   //           Pre_amp_MIP_Range
 316   //    gain = ----------------- * Energy_deposited_per_MIP
 317   //           ADC_channel_size
 318   //
 319   // is constant and the same for all strips.
 320
 321   // Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
 322   // Double_t edep  = TMath::Abs(TMath::Cos(theta)) * fGain * count;
 323   AliFMDParameters* param = AliFMDParameters::Instance();
 324   Float_t           gain  = param->GetPulseGain(digit->Detector(),
 325                                                 digit->Ring(),
 326                                                 digit->Sector(),
 327                                                 digit->Strip());
 328   Double_t          edep  = count * gain;
 329   AliDebug(10, Form("Converting counts %d to energy via factor %f",
 330                     count, gain));
 331   return edep;
 332 }
 333
 334 //____________________________________________________________________
 335 Float_t
 336 AliFMDReconstructor::Energy2Multiplicity(AliFMDDigit* /* digit */,
 337                                          Float_t      edep) const
 338 {
 339   // Converts an energy signal to number of particles.
 340   // Note, that this member function can be overloaded by derived
 341   // classes to do strip-specific look-ups in databases or the like,
 342   // to find the proper gain for a strip.
 343   //
 344   // In this simple version, we calculate the multiplicity as
 345   //
 346   //   multiplicity = Energy_deposited / Energy_deposited_per_MIP
 347   //
 348   // where
 349   //
 350   //   Energy_deposited_per_MIP = 1.664 * SI_density * SI_thickness
 351   //
 352   // is constant and the same for all strips
 353   AliFMDParameters* param   = AliFMDParameters::Instance();
 354   Double_t          edepMIP = param->GetEdepMip();
 355   Float_t           mult    = edep / edepMIP;
 356   if (edep > 0)
 357     AliDebug(10, Form("Translating energy %f to multiplicity via "
 358                      "divider %f->%f", edep, edepMIP, mult));
 359   return mult;
 360 }
 361
 362 //____________________________________________________________________
 363 void
 364 AliFMDReconstructor::PhysicalCoordinates(AliFMDDigit* digit,
 365                                          Float_t& eta,
 366                                          Float_t& phi) const
 367 {
 368   // Get the eta and phi of a digit
 369   //
 370   // Get geometry.
 371   AliFMDGeometry* fmd = AliFMDGeometry::Instance();
 372 #if 0
 373   AliFMDDetector* subDetector = fmd->GetDetector(digit->Detector());
 374   if (!subDetector) {
 375     Warning("ProcessDigits", "Unknown detector: FMD%d" , digit->Detector());
 376     return;
 377   }
 378
 379   // Get the ring - we should really use geometry->Detector2XYZ(...)
 380   // here.
 381   AliFMDRing* ring  = subDetector->GetRing(digit->Ring());
 382   Float_t     ringZ = subDetector->GetRingZ(digit->Ring());
 383   if (!ring) {
 384     Warning("ProcessDigits", "Unknown ring: FMD%d%c", digit->Detector(),
 385             digit->Ring());
 386     return;
 387   }
 388   Float_t  realZ    = fCurrentVertex + ringZ;
 389   Float_t  stripR   = ((ring->GetHighR() - ring->GetLowR())
 390                        / ring->GetNStrips() * (digit->Strip() + .5)
 391                        + ring->GetLowR());
 392   Float_t  theta    = TMath::ATan2(stripR, realZ);
 393 #endif
 394   Double_t x, y, z, r, theta;
 395   fmd->Detector2XYZ(digit->Detector(), digit->Ring(), digit->Sector(),
 396                     digit->Strip(), x, y, z);
 397   // Correct for vertex offset.
 398   z     += fCurrentVertex;
 399   phi   =  TMath::ATan2(y, x);
 400   r     =  TMath::Sqrt(y * y + x * x);
 401   theta =  TMath::ATan2(r, z);
 402   eta   = -TMath::Log(TMath::Tan(theta / 2));
 403 }
 404
 405
 406
 407 //____________________________________________________________________
 408 void
 409 AliFMDReconstructor::FillESD(TTree*  /* digitsTree */,
 410                              TTree*  /* clusterTree */,
 411                              AliESD* esd) const
 412 {
 413   // nothing to be done
 414   // FIXME: The vertex may not be known when Reconstruct is executed,
 415   // so we may have to move some of that member function here.
 416   AliDebug(1, Form("Calling FillESD with two trees and one ESD"));
 417   // fESDObj->Print();
 418
 419   if (esd) {
 420     AliDebug(1, Form("Writing FMD data to ESD tree"));
 421     esd->SetFMDData(fESDObj);
 422     // Let's check the data in the ESD
 423 #if 0
 424     AliESDFMD* fromEsd = esd->GetFMDData();
 425     if (!fromEsd) {
 426       AliWarning("No FMD object in ESD!");
 427       return;
 428     }
 429     for (UShort_t det = 1; det <= fESDObj->MaxDetectors(); det++) {
 430       for (UShort_t ir = 0; ir < fESDObj->MaxRings(); ir++) {
 431         Char_t ring = (ir == 0 ? 'I' : 'O');
 432         for (UShort_t sec = 0; sec < fESDObj->MaxSectors(); sec++) {
 433           for (UShort_t str = 0; str < fESDObj->MaxStrips(); str++) {
 434             if (fESDObj->Multiplicity(det, ring, sec, str) !=
 435                 fromEsd->Multiplicity(det, ring, sec, str))
 436               AliWarning(Form("Mult for FMD%d%c[%2d,%3d]",det,ring,sec,str));
 437             if (fESDObj->Eta(det, ring, sec, str) !=
 438                 fromEsd->Eta(det, ring, sec, str))
 439               AliWarning(Form("Eta for FMD%d%c[%2d,%3d]", det,ring,sec,str));
 440             if (fESDObj->Multiplicity(det, ring, sec, str) > 0 &&
 441                 fESDObj->Multiplicity(det, ring, sec, str)
 442                 != AliESDFMD::kInvalidMult)
 443               AliInfo(Form("Mult in FMD%d%c[%2d,%3d] is %f", det,ring,sec,str,
 444                            fESDObj->Multiplicity(det, ring, sec, str)));
 445           }
 446         }
 447       }
 448     }
 449 #endif
 450   }
 451
 452 #if 0
 453   static Int_t evNo = -1;
 454   evNo++;
 455   if (esd) evNo = esd->GetEventNumber();
 456   TString fname(Form("FMD.ESD.%03d.root", evNo));
 457   TFile* file = TFile::Open(fname.Data(), "RECREATE");
 458   if (!file) {
 459     AliError(Form("Failed to open file %s", fname.Data()));
 460     return;
 461   }
 462   fESDObj->Write();
 463   file->Close();
 464 #endif
 465
 466 #if 0
 467   TClonesArray* multStrips  = 0;
 468   TClonesArray* multRegions = 0;
 469   TTree*        treeR  = fmdLoader->TreeR();
 470   TBranch*      branchRegions = treeR->GetBranch("FMDPoisson");
 471   TBranch*      branchStrips  = treeR->GetBranch("FMDNaiive");
 472   branchRegions->SetAddress(&multRegions);
 473   branchStrips->SetAddress(&multStrips);
 474
 475   Int_t total = 0;
 476   Int_t nEntries  = clusterTree->GetEntries();
 477   for (Int_t entry = 0; entry < nEntries; entry++) {
 478     AliDebug(5, Form("Entry # %d in cluster tree", entry));
 479     treeR->GetEntry(entry);
 480
 481
 482     Int_t nMults = multRegions->GetLast();
 483     for (Int_t i = 0; i <= nMults; i++) {
 484       AliFMDMultRegion* multR =
 485         static_cast<AliFMDMultRegion*>(multRegions->UncheckedAt(i));
 486       Int_t nParticles=multR->Particles();
 487       if (i>=0 && i<=13)   hEtaPoissonI1->AddBinContent(i+1,nParticles);
 488       if (i>=14 && i<=27 ) hEtaPoissonI2->AddBinContent(i-13,nParticles);
 489       if (i>=28 && i<=33 );
 490       if (i>=34 && i<=47 ) hEtaPoissonI3->AddBinContent(48-i,nParticles);
 491       if (i>=48 && i<=53)  hEtaPoissonO3->AddBinContent(54-i,nParticles);
 492     }
 493   }
 494 #endif
 495 }
 496
 497
 498 //____________________________________________________________________
 499 void
 500 AliFMDReconstructor::Reconstruct(AliRawReader*,TTree*) const
 501 {
 502   // Cannot be used.  See member function with same name but with 2
 503   // TTree arguments.   Make sure you do local reconstrucion
 504   AliDebug(1, Form("Calling FillESD with loader and tree"));
 505   AliError("MayNotUse");
 506 }
 507 //____________________________________________________________________
 508 void
 509 AliFMDReconstructor::Reconstruct(AliRunLoader*) const
 510 {
 511   // Cannot be used.  See member function with same name but with 2
 512   // TTree arguments.   Make sure you do local reconstrucion
 513   AliDebug(1, Form("Calling FillESD with loader"));
 514   AliError("MayNotUse");
 515 }
 516 //____________________________________________________________________
 517 void
 518 AliFMDReconstructor::Reconstruct(AliRunLoader*, AliRawReader*) const
 519 {
 520   // Cannot be used.  See member function with same name but with 2
 521   // TTree arguments.   Make sure you do local reconstrucion
 522   AliDebug(1, Form("Calling FillESD with loader and raw reader"));
 523   AliError("MayNotUse");
 524 }
 525 //____________________________________________________________________
 526 void
 527 AliFMDReconstructor::FillESD(AliRawReader*,TTree*,AliESD*) const
 528 {
 529   // Cannot be used.  See member function with same name but with 2
 530   // TTree arguments.   Make sure you do local reconstrucion
 531   AliDebug(1, Form("Calling FillESD with raw reader, tree, and ESD"));
 532   AliError("MayNotUse");
 533 }
 534 //____________________________________________________________________
 535 void
 536 AliFMDReconstructor::FillESD(AliRunLoader*,AliESD*) const
 537 {
 538   // Cannot be used.  See member function with same name but with 2
 539   // TTree arguments.   Make sure you do local reconstrucion
 540   AliDebug(1, Form("Calling FillESD with loader and ESD"));
 541   AliError("MayNotUse");
 542 }
 543 //____________________________________________________________________
 544 void
 545 AliFMDReconstructor::FillESD(AliRunLoader*,AliRawReader*,AliESD*) const
 546 {
 547   // Cannot be used.  See member function with same name but with 2
 548   // TTree arguments.   Make sure you do local reconstrucion
 549   AliDebug(1, Form("Calling FillESD with loader, raw reader, and ESD"));
 550   AliError("MayNotUse");
 551 }
 552
 553 //____________________________________________________________________
 554 //
 555 // EOF
 556 //