]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - FMD/AliFMDReconstructor.cxx
Whoops - stupid mistake - sorry
[u/mrichter/AliRoot.git] / FMD / AliFMDReconstructor.cxx
index 224dac0fcf9fd91e2442c52f30c3605d9071f735..846a59a69b16f3b5bc65edb580f87fc46b3147c8 100644 (file)
@@ -1,3 +1,16 @@
+//____________________________________________________________________
+//
+// This is a class that constructs AliFMDRecPoint objects from of Digits
+// This class reads either digits from a TClonesArray or raw data from 
+// a DDL file (or similar), and stores the read ADC counts in an
+// internal cache (fAdcs).   The rec-points are made via the naiive
+// method. 
+//
+//-- Authors: Evgeny Karpechev(INR) and Alla Maevsksia
+//  Latest changes by Christian Holm Christensen <cholm@nbi.dk>
+//
+//
+//____________________________________________________________________
 /**************************************************************************
  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  *                                                                        *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 /* $Id$ */
-/** @file    AliFMDReconstructor.cxx
-    @author  Christian Holm Christensen <cholm@nbi.dk>
-    @date    Mon Mar 27 12:47:09 2006
-    @brief   FMD reconstruction 
+/** 
+ * @file    AliFMDReconstructor.cxx
+ * @author  Christian Holm Christensen <cholm@nbi.dk>
+ * @date    Mon Mar 27 12:47:09 2006
+ * @brief   FMD reconstruction 
 */
-//____________________________________________________________________
-//
-// This is a class that constructs AliFMDRecPoint objects from of Digits
-// This class reads either digits from a TClonesArray or raw data from 
-// a DDL file (or similar), and stores the read ADC counts in an
-// internal cache (fAdcs).   The rec-points are made via the naiive
-// method. 
-//
-//-- Authors: Evgeny Karpechev(INR) and Alla Maevsksia
-//  Latest changes by Christian Holm Christensen <cholm@nbi.dk>
-//
-//
-//____________________________________________________________________
 
 // #include <AliLog.h>                        // ALILOG_H
 // #include <AliRun.h>                        // ALIRUN_H
@@ -76,12 +77,14 @@ AliFMDReconstructor::AliFMDReconstructor()
     fAngleCorrect(kTRUE),
     fVertexType(kNoVertex),
     fESD(0x0),
-    fDiagnostics(kTRUE),
+    fDiagnostics(kFALSE),
     fDiagStep1(0), 
     fDiagStep2(0),
     fDiagStep3(0),
     fDiagStep4(0),
-    fDiagAll(0)
+    fDiagAll(0),
+    fBad(0), 
+    fZombie(false)
 {
   // Make a new FMD reconstructor object - default CTOR.  
   SetNoiseFactor();
@@ -115,44 +118,59 @@ AliFMDReconstructor::Init()
 
   // Initialize the parameters
   AliFMDParameters* param = AliFMDParameters::Instance();
-  param->Init();
+  if (param->Init() != 0) {
+    AliError("Failed to initialize parameters, making zombie");
+    fZombie = true;
+  }
+  else 
+    fZombie = false;
   
   // Current vertex position
   fCurrentVertex = 0;
   // Create array of reconstructed strip multiplicities 
-  fMult = new TClonesArray("AliFMDRecPoint", 51200);
+  // fMult = new TClonesArray("AliFMDRecPoint", 51200);
   // Create ESD output object 
   fESDObj = new AliESDFMD;
   
   // Check if we need diagnostics histograms 
   if (!fDiagnostics) return;
   AliInfo("Making diagnostics histograms");
-  fDiagStep1   = new TH2I("diagStep1", "Read ADC vs. Noise surpressed ADC",
-                       1024, -.5, 1023.5, 1024, -.5, 1023.5);
-  fDiagStep1->SetDirectory(0);
-  fDiagStep1->GetXaxis()->SetTitle("ADC (read)");
-  fDiagStep1->GetYaxis()->SetTitle(Form("ADC (noise surpressed %4.f)", 
-                                       fNoiseFactor));
-  fDiagStep2  = new TH2F("diagStep2",  "ADC vs Edep deduced",
-                       1024, -.5, 1023.5, 100, 0, 2);
-  fDiagStep2->SetDirectory(0);
-  fDiagStep2->GetXaxis()->SetTitle("ADC (noise surpressed)");
-  fDiagStep2->GetYaxis()->SetTitle("#Delta E [GeV]");
-  fDiagStep3  = new TH2F("diagStep3",  "Edep vs Edep path corrected",
-                       100, 0., 2., 100, 0., 2.);
-  fDiagStep3->SetDirectory(0);
-  fDiagStep3->GetXaxis()->SetTitle("#Delta E [GeV]");
-  fDiagStep3->GetYaxis()->SetTitle("#Delta E/#Delta x #times #delta x [GeV]");
-  fDiagStep4  = new TH2F("diagStep4",  "Edep vs Multiplicity deduced", 
-                       100, 0., 2., 100, -.1, 19.9);
-  fDiagStep4->SetDirectory(0);
-  fDiagStep4->GetXaxis()->SetTitle("#Delta E/#Delta x #times #delta x [GeV]");
-  fDiagStep4->GetYaxis()->SetTitle("Multiplicity");
-  fDiagAll    = new TH2F("diagAll",    "Read ADC vs Multiplicity deduced", 
-                        1024, -.5, 1023.5, 100, -.1, 19.9);
-  fDiagAll->SetDirectory(0);
-  fDiagAll->GetXaxis()->SetTitle("ADC (read)");
-  fDiagAll->GetYaxis()->SetTitle("Multiplicity");
+  if (!fDiagStep1) {
+    fDiagStep1   = new TH2I("diagStep1", "Read ADC vs. Noise surpressed ADC",
+                           1024, -.5, 1023.5, 1024, -.5, 1023.5);
+    fDiagStep1->SetDirectory(0);
+    fDiagStep1->GetXaxis()->SetTitle("ADC (read)");
+    fDiagStep1->GetYaxis()->SetTitle(Form("ADC (noise surpressed %4.f)", 
+                                         fNoiseFactor));
+  }
+  if (!fDiagStep2) {
+    fDiagStep2  = new TH2F("diagStep2",  "ADC vs Edep deduced",
+                          1024, -.5, 1023.5, 100, 0, 2);
+    fDiagStep2->SetDirectory(0);
+    fDiagStep2->GetXaxis()->SetTitle("ADC (noise surpressed)");
+    fDiagStep2->GetYaxis()->SetTitle("#Delta E [GeV]");
+  }
+  if (!fDiagStep3) {
+    fDiagStep3  = new TH2F("diagStep3",  "Edep vs Edep path corrected",
+                          100, 0., 2., 100, 0., 2.);
+    fDiagStep3->SetDirectory(0);
+    fDiagStep3->GetXaxis()->SetTitle("#Delta E [GeV]");
+    fDiagStep3->GetYaxis()->SetTitle("#Delta E/#Delta x #times #delta x [GeV]");
+  }
+  if (!fDiagStep4) {
+    fDiagStep4  = new TH2F("diagStep4",  "Edep vs Multiplicity deduced", 
+                          100, 0., 2., 100, -.1, 19.9);
+    fDiagStep4->SetDirectory(0);
+    fDiagStep4->GetXaxis()->SetTitle("#Delta E/#Delta x #times #delta x [GeV]");
+    fDiagStep4->GetYaxis()->SetTitle("Multiplicity");
+    fDiagAll    = new TH2F("diagAll",    "Read ADC vs Multiplicity deduced", 
+                          1024, -.5, 1023.5, 100, -.1, 19.9);
+  }
+  if (!fDiagAll) {
+    fDiagAll->SetDirectory(0);
+    fDiagAll->GetXaxis()->SetTitle("ADC (read)");
+    fDiagAll->GetYaxis()->SetTitle("Multiplicity");
+  }
 }
 
 //____________________________________________________________________
@@ -161,10 +179,29 @@ AliFMDReconstructor::ConvertDigits(AliRawReader* reader,
                                   TTree* digitsTree) const
 {
   // Convert Raw digits to AliFMDDigit's in a tree 
+  if (fZombie) { 
+    AliWarning("I'm a zombie - cannot do anything");
+    return;
+  }
   AliFMDDebug(1, ("Reading raw data into digits tree"));
+  if (!digitsTree) { 
+    AliError("No digits tree passed");
+    return;
+  }
+  static TClonesArray* array = new TClonesArray("AliFMDDigit");
+  digitsTree->Branch("FMD", &array);
+  array->Clear();
+  
   AliFMDRawReader rawRead(reader, digitsTree);
   // rawRead.SetSampleRate(fFMD->GetSampleRate());
-  rawRead.Exec();
+  // rawRead.Exec();
+  rawRead.ReadAdcs(array);
+
+  Int_t nWrite = digitsTree->Fill();
+  AliFMDDebug(1, ("Got a grand total of %d digits, wrote %d bytes to tree", 
+                  array->GetEntriesFast(), nWrite));
+
+  
   AliFMDAltroMapping* map = AliFMDParameters::Instance()->GetAltroMap();
   for (size_t i = 1; i <= 3; i++) { 
     fZS[i-1]       = rawRead.IsZeroSuppressed(map->Detector2DDL(i));
@@ -208,14 +245,22 @@ AliFMDReconstructor::GetVertex(AliESDEvent* esd) const
 Int_t
 AliFMDReconstructor::GetIdentifier() const
 {
-  return AliReconstruction::GetDetIndex(GetDetectorName());
+  // Get the detector identifier. 
+  // Note the actual value is cached so that we do not 
+  // need to do many expensive string comparisons. 
+  static Int_t idx = AliReconstruction::GetDetIndex(GetDetectorName());
+  return idx;
 }
 
 //____________________________________________________________________
 const AliFMDRecoParam*
 AliFMDReconstructor::GetParameters() const
 {
-  Int_t iDet = 12; // GetIdentifier();
+  // Get the reconstruction parameters. 
+  // 
+  // Return: 
+  //   Pointer to reconstruction parameters or null if not found or wrong type
+  Int_t iDet = GetIdentifier(); // Was 12 - but changed on Cvetans request
   const AliDetectorRecoParam* params = AliReconstructor::GetRecoParam(iDet);
   if (!params || params->IsA() != AliFMDRecoParam::Class()) return 0;
   return static_cast<const AliFMDRecoParam*>(params);
@@ -225,6 +270,14 @@ AliFMDReconstructor::GetParameters() const
 void 
 AliFMDReconstructor::UseRecoParam(Bool_t set) const
 {
+  // 
+  // Set-up reconstructor to use values from reconstruction
+  // parameters, if present, for this event.   If the argument @a set
+  // is @c false, then restore preset values. 
+  // 
+  // Parameters:
+  //    set 
+  //  
   static Float_t savedNoiseFactor  = fNoiseFactor;
   static Bool_t  savedAngleCorrect = fAngleCorrect;
   if (set) { 
@@ -240,20 +293,55 @@ AliFMDReconstructor::UseRecoParam(Bool_t set) const
 }
   
   
-
 //____________________________________________________________________
 void 
-AliFMDReconstructor::Reconstruct(AliRawReader* reader, TTree*) const
+AliFMDReconstructor::MarkDeadChannels(AliESDFMD* esd) const
 {
-  // Reconstruct directly from raw data (no intermediate output on
-  // digit tree or rec point tree).  
-  // 
-  // Parameters: 
-  //   reader  Raw event reader 
-  //   ctree    Not used. 
-  AliFMDDebug(1, ("Reconstructing from raw reader"));
-  AliFMDRawReader rawReader(reader, 0);
+  // Loop over all entries of the ESD and mark 
+  // those that are dead as such 
+  // - otherwise put in the zero signal. 
+  AliFMDParameters* param = AliFMDParameters::Instance();
+
+  for (UShort_t d = 1; d <= 3; d++) { 
+    UShort_t nR = (d == 1 ? 1 : 2);
+
+    for (UShort_t q = 0; q < nR; q++) {
+      Char_t   r  = (q == 0 ? 'I' : 'O');
+      UShort_t nS = (q == 0 ?  20 :  40);
+      UShort_t nT = (q == 0 ? 512 : 256);
+
+      for (UShort_t s = 0; s < nS; s++) { 
+       for (UShort_t t = 0; t < nT; t++) {
+         if (fBad(d, r, s, t)) { 
+           AliDebug(5, Form("Marking FMD%d%c[%2d,%3d] as bad", d, r, s, t));
+           esd->SetMultiplicity(d, r, s, t, AliESDFMD::kInvalidMult);
+         }         
+         if (param->IsDead(d, r, s, t)) { 
+           AliDebug(5, Form("Marking FMD%d%c[%2d,%3d] as dead", d, r, s, t));
+           esd->SetMultiplicity(d, r, s, t, AliESDFMD::kInvalidMult);
+           // esd->SetEta(d, r, s, t, AliESDFMD::kInvalidEta);
+         }
+         else if (esd->Multiplicity(d, r, s, t) == AliESDFMD::kInvalidMult) {
+           AliDebug(10, Form("Setting null signal in FMD%d%c[%2d,%3d]", 
+                             d, r, s, t));
+           esd->SetMultiplicity(d, r, s, t, 0);
+         }
+       }
+      }
+    }
+  }
+}
 
+//____________________________________________________________________
+void 
+AliFMDReconstructor::Reconstruct(AliFMDRawReader& rawReader) const
+{
+  AliFMDDebug(1, ("Reconstructing from FMD raw reader"));
+  if (fZombie) { 
+    AliWarning("I'm a zombie - cannot do anything");
+    return;
+  }
+  fBad.Reset(false);
   UShort_t det, sec, str, fac;
   Short_t  adc, oldDet = -1;
   Bool_t   zs;
@@ -269,7 +357,26 @@ AliFMDReconstructor::Reconstruct(AliRawReader* reader, TTree*) const
     ProcessSignal(det, rng, sec, str, adc);
   }
   UseRecoParam(kFALSE);
-  
+
+}
+
+//____________________________________________________________________
+void 
+AliFMDReconstructor::Reconstruct(AliRawReader* reader, TTree*) const
+{
+  // Reconstruct directly from raw data (no intermediate output on
+  // digit tree or rec point tree).  
+  // 
+  // Parameters: 
+  //   reader  Raw event reader 
+  //   ctree    Not used - 'cluster tree' to store rec-points in. 
+  AliFMDDebug(1, ("Reconstructing from raw reader"));
+  if (fZombie) { 
+    AliWarning("I'm a zombie - cannot do anything");
+    return;
+  }
+  AliFMDRawReader rawReader(reader, 0);
+  Reconstruct(rawReader);
 }
 
 //____________________________________________________________________
@@ -282,6 +389,10 @@ AliFMDReconstructor::Digitize(AliRawReader* reader, TClonesArray* sdigits) const
   // Parameters: 
   //   reader  Raw event reader 
   //   ctree    Not used. 
+  if (fZombie) { 
+    AliWarning("I'm a zombie - cannot do anything");
+    return;
+  }
   AliFMDRawReader rawReader(reader, 0);
 
   UShort_t det, sec, str, sam, rat, fac;
@@ -316,18 +427,26 @@ AliFMDReconstructor::Reconstruct(TTree* digitsTree,
   //   digitsTree      Pointer to a tree containing digits 
   //   clusterTree     Pointer to output tree 
   // 
+  if (fZombie) { 
+    AliWarning("I'm a zombie - cannot do anything");
+    return;
+  }
+  if (!fMult) fMult = new TClonesArray("AliFMDRecPoint");
+
   AliFMDDebug(1, ("Reconstructing from digits in a tree"));
   GetVertex(fESD);
 
+
+
   static TClonesArray* digits = new TClonesArray("AliFMDDigit");
-  TBranch *digitBranch = digitsTree->GetBranch("FMD");
+  TBranch*      digitBranch   = digitsTree->GetBranch("FMD");
   if (!digitBranch) {
     Error("Exec", "No digit branch for the FMD found");
     return;
   }
-  // TClonesArray* digits = new TClonesArray("AliFMDDigit");
   digitBranch->SetAddress(&digits);
 
+  if (digits)  digits->Clear();
   if (fMult)   fMult->Clear();
   if (fESDObj) fESDObj->Clear();
   
@@ -345,11 +464,37 @@ AliFMDReconstructor::Reconstruct(TTree* digitsTree,
 
   Int_t written = clusterTree->Fill();
   AliFMDDebug(10, ("Filled %d bytes into cluster tree", written));
-  digits->Delete();
+  // digits->Delete();
   // delete digits;
 }
  
 
+//____________________________________________________________________
+void
+AliFMDReconstructor::ProcessDigits(TClonesArray* digits, 
+                                  const AliFMDRawReader& rawRead) const
+{
+  // For each digit, find the pseudo rapdity, azimuthal angle, and
+  // number of corrected ADC counts, and pass it on to the algorithms
+  // used. 
+  // 
+  // Parameters: 
+  //    digits Array of digits
+  // 
+  if (fZombie) { 
+    AliWarning("I'm a zombie - cannot do anything");
+    return;
+  }
+  AliFMDAltroMapping* map = AliFMDParameters::Instance()->GetAltroMap();
+  for (size_t i = 1; i <= 3; i++) { 
+    fZS[i-1]       = rawRead.IsZeroSuppressed(map->Detector2DDL(i));
+    fZSFactor[i-1] = rawRead.NoiseFactor(map->Detector2DDL(i));
+  }
+  UseRecoParam(kTRUE);
+  ProcessDigits(digits);
+  UseRecoParam(kFALSE);
+}
+
 //____________________________________________________________________
 void
 AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
@@ -365,6 +510,7 @@ AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
   AliFMDDebug(2, ("Got %d digits", nDigits));
   fESDObj->SetNoiseFactor(fNoiseFactor);
   fESDObj->SetAngleCorrected(fAngleCorrect);
+  fBad.Reset(false);
   for (Int_t i = 0; i < nDigits; i++) {
     AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
     if (!digit) continue;
@@ -376,6 +522,12 @@ AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
 void
 AliFMDReconstructor::ProcessDigit(AliFMDDigit* digit) const
 {
+  // 
+  // Process a single digit 
+  // 
+  // Parameters:
+  //    digit Digiti to process
+  // 
   UShort_t det = digit->Detector();
   Char_t   rng = digit->Ring();
   UShort_t sec = digit->Sector();
@@ -402,10 +554,17 @@ AliFMDReconstructor::ProcessSignal(UShort_t det,
   //    rng    Strip ID
   //    adc     ADC counts
   // 
-  AliFMDParameters* param  = AliFMDParameters::Instance();
+  if (adc >= AliFMDRawReader::kBadSignal) { 
+    AliFMDDebug(1, ("FMD%d%c[%2d,%3d] is marked bad", det, rng, sec, str));
+    fBad(det,rng,sec,str) = true;
+    return;
+  }
+
   // Check that the strip is not marked as dead 
+  AliFMDParameters* param  = AliFMDParameters::Instance();
   if (param->IsDead(det, rng, sec, str)) {
     AliFMDDebug(1, ("FMD%d%c[%2d,%3d] is dead", det, rng, sec, str));
+    fBad(det,rng,sec,str) = true;
     return;
   }
   
@@ -505,6 +664,25 @@ AliFMDReconstructor::SubtractPedestal(UShort_t det,
                                      Bool_t   zsEnabled, 
                                      UShort_t zsNoiseFactor) const
 {
+  // 
+  // Subtract the pedestal off the ADC counts. 
+  // 
+  // Parameters:
+  //    det           Detector number
+  //    rng           Ring identifier
+  //    sec           Sector number
+  //    str           Strip number
+  //    adc           ADC counts
+  //    noiseFactor   If pedestal substracted pedestal is less then
+  //        this times the noise, then consider this to be 0. 
+  //    zsEnabled     Whether zero-suppression is on.
+  //    zsNoiseFactor Noise factor used in on-line pedestal
+  //        subtraction. 
+  // 
+  // Return:
+  //    The pedestal subtracted ADC counts (possibly 0), or @c
+  //         USHRT_MAX in case of problems.
+  //  
   AliFMDParameters* param  = AliFMDParameters::Instance();
   Float_t           ped    = (zsEnabled ? 0 : 
                                param->GetPedestal(det, rng, sec, str));
@@ -758,17 +936,31 @@ AliFMDReconstructor::PhysicalCoordinates(UShort_t det,
   //    phi     On return, contains the azimuthal angle of the strip
   // 
   AliFMDGeometry* geom = AliFMDGeometry::Instance();
-  Double_t x, y, z, r, theta;
+  Double_t x, y, z, r, theta, deta, dphi;
   geom->Detector2XYZ(det, rng, sec, str, x, y, z);
+
   // Correct for vertex offset. 
   z     -= fCurrentVertex;
-  phi   =  TMath::ATan2(y, x);
-  r     =  TMath::Sqrt(y * y + x * x);
-  theta =  TMath::ATan2(r, z);
-  eta   = -TMath::Log(TMath::Tan(theta / 2));
+  AliFMDGeometry::XYZ2REtaPhiTheta(x, y, z, r, deta, dphi, theta);
+  eta = deta;
+  phi = dphi;
+}
+
+namespace { 
+  class ESDPrinter : public AliESDFMD::ForOne
+  {
+  public:
+    ESDPrinter() {}
+    Bool_t operator()(UShort_t d, Char_t r, UShort_t s, UShort_t t, 
+                     Float_t m, Float_t e)
+    {
+      if (m > 0 && m != AliESDFMD::kInvalidMult) 
+       printf("  FMD%d%c[%2d,%3d] = %6.3f / %6.3f\n", d, r, s, t, m, e);
+      return kTRUE;
+    }
+  };
 }
 
-      
 
 //____________________________________________________________________
 void 
@@ -780,8 +972,15 @@ AliFMDReconstructor::FillESD(TTree*  /* digitsTree */,
   // FIXME: The vertex may not be known when Reconstruct is executed,
   // so we may have to move some of that member function here. 
   AliFMDDebug(2, ("Calling FillESD with two trees and one ESD"));
+  if (fZombie) { 
+    AliWarning("I'm a zombie - cannot do anything");
+    return;
+  }
   // fESDObj->Print();
 
+  // Fix up ESD so that only truely dead channels get the kInvalidMult flag. 
+  MarkDeadChannels(fESDObj);
+
   Double_t oldVz = fCurrentVertex;
   GetVertex(esd);
   if (fVertexType != kNoVertex) { 
@@ -790,6 +989,11 @@ AliFMDReconstructor::FillESD(TTree*  /* digitsTree */,
     AliFMDESDRevertexer revertexer;
     revertexer.Revertex(fESDObj, fCurrentVertex);
   }
+  
+  if (AliDebugLevel() > 10) { 
+    ESDPrinter p;
+    fESDObj->ForEach(p);
+  }
 
   if (esd) { 
     AliFMDDebug(2, ("Writing FMD data to ESD tree"));
@@ -831,6 +1035,13 @@ void
 AliFMDReconstructor::FillESD(AliRawReader*, TTree* clusterTree, 
                             AliESDEvent* esd) const
 {
+  // 
+  // Forwards to above member function 
+  //
+  if (fZombie) { 
+    AliWarning("I'm a zombie - cannot do anything");
+    return;
+  }
   TTree* dummy = 0;
   FillESD(dummy, clusterTree, esd);
 }