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