]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - FMD/AliFMDReconstructor.cxx
setting the period number also for ESDHeader
[u/mrichter/AliRoot.git] / FMD / AliFMDReconstructor.cxx
index 327cde400318e0dcd3da2e01b6c0a0a1535f4457..3befa0df25f789ae22ab5c9d023bba3f4e35c6c5 100644 (file)
@@ -41,6 +41,7 @@
 #include "AliFMDDigit.h"                   // ALIFMDDIGIT_H
 #include "AliFMDSDigit.h"                  // ALIFMDDIGIT_H
 #include "AliFMDReconstructor.h"           // ALIFMDRECONSTRUCTOR_H
+#include "AliFMDRecoParam.h"               // ALIFMDRECOPARAM_H
 #include "AliFMDRawReader.h"               // ALIFMDRAWREADER_H
 #include "AliFMDRecPoint.h"               // ALIFMDMULTNAIIVE_H
 #include "AliESDEvent.h"                  // ALIESDEVENT_H
 #include <TH2.h>
 #include <TFile.h>
 #include <climits>
-// Import revertexer into a private namespace (to prevent conflicts) 
-namespace { 
-# include "AliFMDESDRevertexer.h"
-}
+#include "AliFMDESDRevertexer.h"
 
 
 class AliRawReader;
@@ -78,7 +76,7 @@ AliFMDReconstructor::AliFMDReconstructor()
     fAngleCorrect(kTRUE),
     fVertexType(kNoVertex),
     fESD(0x0),
-    fDiagnostics(kTRUE),
+    fDiagnostics(kFALSE),
     fDiagStep1(0), 
     fDiagStep2(0),
     fDiagStep3(0),
@@ -89,53 +87,10 @@ AliFMDReconstructor::AliFMDReconstructor()
   SetNoiseFactor();
   SetAngleCorrect();
   if (AliDebugLevel() > 0) fDiagnostics = kTRUE;
-}
-  
-
-//____________________________________________________________________
-AliFMDReconstructor::AliFMDReconstructor(const AliFMDReconstructor& other) 
-  : AliReconstructor(), 
-    fMult(other.fMult),
-    fNMult(other.fNMult),
-    fTreeR(other.fTreeR),
-    fCurrentVertex(other.fCurrentVertex),
-    fESDObj(other.fESDObj),
-    fNoiseFactor(other.fNoiseFactor),
-    fAngleCorrect(other.fAngleCorrect),
-    fVertexType(other.fVertexType),
-    fESD(other.fESD),
-    fDiagnostics(other.fDiagnostics),
-    fDiagStep1(other.fDiagStep1), 
-    fDiagStep2(other.fDiagStep2),
-    fDiagStep3(other.fDiagStep3),
-    fDiagStep4(other.fDiagStep4),
-    fDiagAll(other.fDiagAll) 
-{
-  // Copy constructor 
-}
-  
-
-//____________________________________________________________________
-AliFMDReconstructor&
-AliFMDReconstructor::operator=(const AliFMDReconstructor& other) 
-{
-  // Assignment operator
-  fMult          = other.fMult;
-  fNMult         = other.fNMult;
-  fTreeR         = other.fTreeR;
-  fCurrentVertex = other.fCurrentVertex;
-  fESDObj        = other.fESDObj;
-  fNoiseFactor   = other.fNoiseFactor;
-  fAngleCorrect  = other.fAngleCorrect;
-  fVertexType    = other.fVertexType;
-  fESD           = other.fESD;
-  fDiagnostics   = other.fDiagnostics;
-  fDiagStep1     = other.fDiagStep1;
-  fDiagStep2     = other.fDiagStep2;
-  fDiagStep3     = other.fDiagStep3;
-  fDiagStep4     = other.fDiagStep4;
-  fDiagAll       = other.fDiagAll;
-  return *this;
+  for(Int_t det = 1; det<=3; det++) {
+    fZS[det-1]       = kFALSE;
+    fZSFactor[det-1] = 0;
+  }
 }
 
 //____________________________________________________________________
@@ -165,7 +120,7 @@ AliFMDReconstructor::Init()
   // 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;
   
@@ -207,13 +162,27 @@ AliFMDReconstructor::ConvertDigits(AliRawReader* reader,
 {
   // Convert Raw digits to AliFMDDigit's in a tree 
   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);
+  
   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]       = rawRead.IsZeroSuppressed(map->Detector2DDL(i));
-    fZSFactor[i] = rawRead.NoiseFactor(map->Detector2DDL(i));
+    fZS[i-1]       = rawRead.IsZeroSuppressed(map->Detector2DDL(i));
+    fZSFactor[i-1] = rawRead.NoiseFactor(map->Detector2DDL(i));
   }
 }
 
@@ -249,6 +218,70 @@ AliFMDReconstructor::GetVertex(AliESDEvent* esd) const
   AliWarning("Didn't get any vertex from ESD or generator");
 }
   
+//____________________________________________________________________
+Int_t
+AliFMDReconstructor::GetIdentifier() const
+{
+  return AliReconstruction::GetDetIndex(GetDetectorName());
+}
+
+//____________________________________________________________________
+const AliFMDRecoParam*
+AliFMDReconstructor::GetParameters() const
+{
+  Int_t iDet = 12; // GetIdentifier();
+  const AliDetectorRecoParam* params = AliReconstructor::GetRecoParam(iDet);
+  if (!params || params->IsA() != AliFMDRecoParam::Class()) return 0;
+  return static_cast<const AliFMDRecoParam*>(params);
+}
+
+//____________________________________________________________________
+void 
+AliFMDReconstructor::UseRecoParam(Bool_t set) const
+{
+  static Float_t savedNoiseFactor  = fNoiseFactor;
+  static Bool_t  savedAngleCorrect = fAngleCorrect;
+  if (set) { 
+    const AliFMDRecoParam* params  = GetParameters();
+    if (params) { 
+      fNoiseFactor  = params->NoiseFactor();
+      fAngleCorrect = params->AngleCorrect();
+    }
+    return;
+  }
+  fNoiseFactor  = savedNoiseFactor;
+  fAngleCorrect = savedAngleCorrect;
+}
+  
+  
+//____________________________________________________________________
+void 
+AliFMDReconstructor::MarkDeadChannels(AliESDFMD* esd) const
+{
+  // 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 (param->IsDead(d, r, s, t)) 
+           esd->SetMultiplicity(d, r, s, t, AliESDFMD::kInvalidMult);
+         else if (esd->Multiplicity(d, r, s, t) == AliESDFMD::kInvalidMult) 
+           esd->SetMultiplicity(d, r, s, t, 0);
+       }
+      }
+    }
+  }
+}
 
 //____________________________________________________________________
 void 
@@ -259,14 +292,16 @@ AliFMDReconstructor::Reconstruct(AliRawReader* reader, TTree*) const
   // 
   // Parameters: 
   //   reader  Raw event reader 
-  //   ctree    Not used. 
+  //   ctree    Not used - 'cluster tree' to store rec-points in. 
+  AliFMDDebug(1, ("Reconstructing from raw reader"));
   AliFMDRawReader rawReader(reader, 0);
 
   UShort_t det, sec, str, fac;
   Short_t  adc, oldDet = -1;
   Bool_t   zs;
   Char_t   rng;
-    
+  UseRecoParam(kTRUE);
   while (rawReader.NextSignal(det, rng, sec, str, adc, zs, fac)) { 
     if (det != oldDet) { 
       fZS[det-1]       = zs;
@@ -275,6 +310,8 @@ AliFMDReconstructor::Reconstruct(AliRawReader* reader, TTree*) const
     }
     ProcessSignal(det, rng, sec, str, adc);
   }
+  UseRecoParam(kFALSE);
+  
 }
 
 //____________________________________________________________________
@@ -294,6 +331,7 @@ AliFMDReconstructor::Digitize(AliRawReader* reader, TClonesArray* sdigits) const
   Bool_t   zs;
   Char_t   rng;
     
+  UseRecoParam(kTRUE);
   while (rawReader.NextSample(det, rng, sec, str, sam, rat, adc, zs, fac)) { 
     if (!rawReader.SelectSample(sam, rat)) continue;
     if (det != oldDet) { 
@@ -303,6 +341,7 @@ AliFMDReconstructor::Digitize(AliRawReader* reader, TClonesArray* sdigits) const
     }
     DigitizeSignal(sdigits, det, rng, sec, str, sam, adc);
   }
+  UseRecoParam(kFALSE);
 }
 
 //____________________________________________________________________
@@ -319,17 +358,22 @@ AliFMDReconstructor::Reconstruct(TTree* digitsTree,
   //   digitsTree      Pointer to a tree containing digits 
   //   clusterTree     Pointer to output tree 
   // 
-  AliFMDDebug(2, ("Reconstructing from digits in a tree"));
+  if (!fMult) fMult = new TClonesArray("AliFMDRecPoint");
+
+  AliFMDDebug(1, ("Reconstructing from digits in a tree"));
   GetVertex(fESD);
 
-  TBranch *digitBranch = digitsTree->GetBranch("FMD");
+
+
+  static TClonesArray* digits = new TClonesArray("AliFMDDigit");
+  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();
   
@@ -340,13 +384,15 @@ AliFMDReconstructor::Reconstruct(TTree* digitsTree,
   AliFMDDebug(5, ("Getting entry 0 from digit branch"));
   digitBranch->GetEntry(0);
   
-  AliFMDDebug(1, ("Processing digits"));
+  AliFMDDebug(5, ("Processing digits"));
+  UseRecoParam(kTRUE);
   ProcessDigits(digits);
+  UseRecoParam(kFALSE);
 
   Int_t written = clusterTree->Fill();
   AliFMDDebug(10, ("Filled %d bytes into cluster tree", written));
-  digits->Delete();
-  delete digits;
+  // digits->Delete();
+  // delete digits;
 }
  
 
@@ -362,7 +408,7 @@ AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
   //    digits Array of digits
   // 
   Int_t nDigits = digits->GetEntries();
-  AliFMDDebug(1, ("Got %d digits", nDigits));
+  AliFMDDebug(2, ("Got %d digits", nDigits));
   fESDObj->SetNoiseFactor(fNoiseFactor);
   fESDObj->SetAngleCorrected(fAngleCorrect);
   for (Int_t i = 0; i < nDigits; i++) {
@@ -381,7 +427,7 @@ AliFMDReconstructor::ProcessDigit(AliFMDDigit* digit) const
   UShort_t sec = digit->Sector();
   UShort_t str = digit->Strip();
   Short_t  adc = digit->Counts();
-  
   ProcessSignal(det, rng, sec, str, adc);
 }
 
@@ -405,7 +451,7 @@ AliFMDReconstructor::ProcessSignal(UShort_t det,
   AliFMDParameters* param  = AliFMDParameters::Instance();
   // Check that the strip is not marked as dead 
   if (param->IsDead(det, rng, sec, str)) {
-    AliFMDDebug(10, ("FMD%d%c[%2d,%3d] is dead", det, rng, sec, str));
+    AliFMDDebug(1, ("FMD%d%c[%2d,%3d] is dead", det, rng, sec, str));
     return;
   }
   
@@ -426,8 +472,13 @@ AliFMDReconstructor::ProcessSignal(UShort_t det,
   // Make rough multiplicity 
   Double_t mult     = Energy2Multiplicity(det, rng, sec, str, edep);
   // Get rid of nonsense mult
+  //if (mult > 20) { 
+  //  AliWarning(Form("The mutliplicity in FMD%d%c[%2d,%3d]=%f > 20 "
+  //               "(ADC: %d, Energy: %f)", det, rng, sec, str, mult, 
+  //               counts, edep));
+  // }
   if (mult < 0)  return; 
-  AliFMDDebug(5, ("FMD%d%c[%2d,%3d]: "
+  AliFMDDebug(10, ("FMD%d%c[%2d,%3d]: "
                    "ADC: %d, Counts: %d, Energy: %f, Mult: %f",
                  det, rng, sec, str, adc, counts, edep, mult));
   
@@ -538,7 +589,8 @@ AliFMDReconstructor::SubtractPedestal(UShort_t det,
   if (counts < noise * nf) counts = 0;
   if (counts > 0) AliDebugClass(15, "Got a hit strip");
 
-  return counts;
+  UShort_t ret = counts < 0 ? 0 : counts;
+  return ret;
 }
 
 
@@ -752,14 +804,14 @@ 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));
+  z     -= fCurrentVertex;
+  AliFMDGeometry::XYZ2REtaPhiTheta(x, y, z, r, deta, dphi, theta);
+  eta = deta;
+  phi = dphi;
 }
 
       
@@ -776,6 +828,9 @@ AliFMDReconstructor::FillESD(TTree*  /* digitsTree */,
   AliFMDDebug(2, ("Calling FillESD with two trees and one ESD"));
   // 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) { 
@@ -797,7 +852,7 @@ AliFMDReconstructor::FillESD(TTree*  /* digitsTree */,
   // - That's OK.  We just use it for the name of the directory -
   // nothing else.  Christian
   Int_t evno = esd->GetEventNumberInFile(); 
-  AliFMDDebug(1, ("Writing diagnostics histograms to FMD.Diag.root/%03d",evno));
+  AliFMDDebug(3, ("Writing diagnostics histograms to FMD.Diag.root/%03d",evno));
   TFile f("FMD.Diag.root", (first ? "RECREATE" : "UPDATE"));
   first = false;
   f.cd(); 
@@ -828,7 +883,6 @@ AliFMDReconstructor::FillESD(AliRawReader*, TTree* clusterTree,
   TTree* dummy = 0;
   FillESD(dummy, clusterTree, esd);
 }
-
 //____________________________________________________________________
 //
 // EOF