Added script to draw calibrated raw data
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 15 Jul 2009 12:29:34 +0000 (12:29 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 15 Jul 2009 12:29:34 +0000 (12:29 +0000)
Added possible use of FMD RecoParam
Implemented FMD RecoParam

FMD/AliFMDRecoParam.cxx [new file with mode: 0644]
FMD/AliFMDRecoParam.h [new file with mode: 0644]
FMD/AliFMDReconstructor.cxx
FMD/AliFMDReconstructor.h
FMD/FMDrecLinkDef.h
FMD/libFMDrec.pkg
FMD/scripts/DrawCalibRaw.C [new file with mode: 0644]

diff --git a/FMD/AliFMDRecoParam.cxx b/FMD/AliFMDRecoParam.cxx
new file mode 100644 (file)
index 0000000..4b7ed95
--- /dev/null
@@ -0,0 +1,83 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+
+///////////////////////////////////////////////////////////////////////////////
+//
+// FMD Reconstruction Parameters 
+//
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
+#include "AliFMDRecoParam.h"
+
+ClassImp(AliDetectorRecoParam)
+#if 0 
+; // Don't delete - for Emacs
+#endif
+
+//____________________________________________________________________
+AliFMDRecoParam::AliFMDRecoParam(Float_t noiseFactor, 
+                                Bool_t angleCorrect,
+                                Bool_t sharingCorrect)
+  : AliDetectorRecoParam(),
+    fNoiseFactor(noiseFactor), 
+    fAngleCorrect(angleCorrect), 
+    fSharingCorrect(sharingCorrect)
+{
+  // Constructor
+  SetName("FMD");
+  SetTitle("FMD");
+}
+
+//____________________________________________________________________
+AliFMDRecoParam*
+AliFMDRecoParam::GetLowFluxParam()
+{
+  AliFMDRecoParam* p = new AliFMDRecoParam(10, kTRUE, kFALSE);
+  p->SetName("FMD_low_flux");
+  p->SetTitle("FMD low flux");
+  return p;
+}
+//____________________________________________________________________
+AliFMDRecoParam*
+AliFMDRecoParam::GetHighFluxParam()
+{
+  AliFMDRecoParam* p = new AliFMDRecoParam(10, kTRUE, kFALSE);
+  p->SetName("FMD_high_flux");
+  p->SetTitle("FMD high flux");
+  return p;
+}
+//____________________________________________________________________
+AliFMDRecoParam*
+AliFMDRecoParam::GetParam(AliRecoParam::EventSpecie_t specie)
+{
+  switch (specie) { 
+  case AliRecoParam::kDefault: 
+  case AliRecoParam::kCalib: 
+  case AliRecoParam::kHighMult: return GetHighFluxParam(); 
+  case AliRecoParam::kCosmic: 
+  case AliRecoParam::kLowMult:  return GetLowFluxParam(); 
+  }
+  return new AliFMDRecoParam();
+}
+
+
+//____________________________________________________________________
+//
+//
+// EOF
+//
diff --git a/FMD/AliFMDRecoParam.h b/FMD/AliFMDRecoParam.h
new file mode 100644 (file)
index 0000000..5a61487
--- /dev/null
@@ -0,0 +1,64 @@
+#ifndef ALIFMDRECOPARAM_H
+#define ALIFMDRECOPARAM_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+//
+//
+// FMD reconstruction parameters
+//
+//
+
+#include "AliDetectorRecoParam.h"
+
+class AliFMDRecoParam : public AliDetectorRecoParam
+{
+public: 
+  AliFMDRecoParam(Float_t noiseFactor=3, 
+                 Bool_t angleCorrect=kTRUE,
+                 Bool_t sharingCorrect=kFALSE);
+  virtual ~AliFMDRecoParam() {}
+  /** 
+   * Whether to do angle of passage correction 
+   * 
+   * @return @c true if we're to do angle of passage correction
+   */
+  Bool_t  AngleCorrect()   const { return fAngleCorrect; }
+  /** 
+   * Get the noise suppression factor
+   * 
+   * @return The number of noise levels away from the pedestal 
+   *         that are suppressed. 
+   */
+  Float_t NoiseFactor()    const { return fNoiseFactor; }
+  /** 
+   * Whether to do the sharing correction.  A single particle may
+   * traverse more than one strip due to it's incident angle.  In that
+   * case, part of it's signal is put in one strip, and another in
+   * it's adjacent strip.  The sharing correction corrects for this
+   * and adds the signal of the two strips into a single strip. 
+   * 
+   * @return @c true if the reconstruction should also do the sharing
+   * correction. 
+   */
+  Bool_t  SharingCorrect() const { return fSharingCorrect; }
+
+  void SetAngleCorrect(Bool_t doit) { fAngleCorrect = doit; }
+  void SetSharingCorrect(Bool_t doit) { fSharingCorrect = doit; }
+  void SetNoiseFactor(Float_t f) { fNoiseFactor = f; }
+
+  static AliFMDRecoParam* GetLowFluxParam();
+  static AliFMDRecoParam* GetHighFluxParam();
+  static AliFMDRecoParam* GetParam(AliRecoParam::EventSpecie_t specie);
+private:
+  Float_t fNoiseFactor;    // Noise suppression factor 
+  Bool_t  fAngleCorrect;   // Whether to do angle correction or not
+  Bool_t  fSharingCorrect; // Whether to do sharing correction or not
+  
+  ClassDef(AliFMDRecoParam, 1)
+};
+
+
+#endif
+// Local Variables:
+//  mode: C++ 
+// End:
index 19b4f12..17211be 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
@@ -253,6 +254,42 @@ 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 
@@ -270,7 +307,8 @@ AliFMDReconstructor::Reconstruct(AliRawReader* reader, TTree*) const
   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;
@@ -279,6 +317,8 @@ AliFMDReconstructor::Reconstruct(AliRawReader* reader, TTree*) const
     }
     ProcessSignal(det, rng, sec, str, adc);
   }
+  UseRecoParam(kFALSE);
+  
 }
 
 //____________________________________________________________________
@@ -298,6 +338,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) { 
@@ -307,6 +348,7 @@ AliFMDReconstructor::Digitize(AliRawReader* reader, TClonesArray* sdigits) const
     }
     DigitizeSignal(sdigits, det, rng, sec, str, sam, adc);
   }
+  UseRecoParam(kFALSE);
 }
 
 //____________________________________________________________________
@@ -345,7 +387,9 @@ AliFMDReconstructor::Reconstruct(TTree* digitsTree,
   digitBranch->GetEntry(0);
   
   AliFMDDebug(1, ("Processing digits"));
+  UseRecoParam(kTRUE);
   ProcessDigits(digits);
+  UseRecoParam(kFALSE);
 
   Int_t written = clusterTree->Fill();
   AliFMDDebug(10, ("Filled %d bytes into cluster tree", written));
index 5837fae..e7e4bf2 100644 (file)
@@ -33,6 +33,7 @@ class AliFMDDigit;
 class AliRawReader;
 class AliESDEvent;
 class AliESDFMD;
+class AliFMDRecoParam;
 class TH1;
 
 
@@ -362,7 +363,28 @@ protected:
                                       UShort_t str, 
                                       Float_t& eta, 
                                       Float_t& phi) 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. 
+   * 
+   * @param set 
+   */  
+  virtual void UseRecoParam(Bool_t set=kTRUE) const;
+  /** 
+   * Utility member function to get the reconstruction parameters for 
+   * this event
+   * 
+   * @return Pointer to AliFMDRecoParam object or null if not
+   * available. 
+   */
+  const AliFMDRecoParam* GetParameters() const;
+  /** 
+   * Get the numeric identifier of this detector
+   * 
+   * @return Should be 12
+   */  
+  Int_t GetIdentifier() const;
   enum Vertex_t {
     kNoVertex,   // Got no vertex
     kGenVertex,  // Got generator vertex 
@@ -373,8 +395,8 @@ protected:
   mutable TTree*        fTreeR;         // Output tree 
   mutable Float_t       fCurrentVertex; // Z-coordinate of primary vertex
   mutable AliESDFMD*    fESDObj;        // ESD output object
-  Float_t               fNoiseFactor;   // Factor of noise to check
-  Bool_t                fAngleCorrect;  // Whether to angle correct
+  mutable Float_t       fNoiseFactor;   // Factor of noise to check
+  mutable Bool_t        fAngleCorrect;  // Whether to angle correct
   mutable Vertex_t      fVertexType;    // What kind of vertex we got
   AliESDEvent*          fESD;           // ESD object(?)
   Bool_t                fDiagnostics;   // Wheter to do diagnostics
index 6613cd3..0ec214b 100644 (file)
@@ -20,6 +20,7 @@
 // #pragma link C++ class  AliFMDMap<UShort_t>;
 // #pragma link C++ typedef AliFMDAdcMap;
 #pragma link C++ class  AliFMDReconstructor+;
+#pragma link C++ class  AliFMDRecoParam+;
 #pragma link C++ class  AliFMDRecPoint+;
 #pragma link C++ class  AliFMDRawStream+;
 #pragma link C++ class  AliFMDRawReader+;
index 183861d..888d4af 100644 (file)
@@ -3,6 +3,7 @@
 # $Id$
 
 SRCS           =  AliFMDReconstructor.cxx      \
+                  AliFMDRecoParam.cxx          \
                   AliFMDRawStream.cxx          \
                   AliFMDRawReader.cxx          \
                   AliFMDRecPoint.cxx           \
diff --git a/FMD/scripts/DrawCalibRaw.C b/FMD/scripts/DrawCalibRaw.C
new file mode 100644 (file)
index 0000000..fc46cd0
--- /dev/null
@@ -0,0 +1,280 @@
+//____________________________________________________________________
+//
+// $Id: DrawDigits.C 30718 2009-01-22 16:07:40Z cholm $
+//
+// Script that contains a class to draw eloss from hits, versus ADC
+// counts from digits, using the AliFMDInputHits class in the util library. 
+//
+// It draws the energy loss versus the p/(mq^2).  It can be overlayed
+// with the Bethe-Bloc curve to show how the simulation behaves
+// relative to the expected. 
+//
+// Use the script `Compile.C' to compile this class using ACLic. 
+//
+#include <TH1D.h>
+#include <AliCDBManager.h>
+#include <AliFMDHit.h>
+#include <AliFMDDigit.h>
+#include <AliFMDInput.h>
+#include <AliFMDEdepMap.h>
+#include <AliFMDParameters.h>
+#include <AliFMDCalibStripRange.h>
+#include <AliFMDCalibSampleRate.h>
+#include <AliFMDCalibGain.h>
+#include <AliFMDCalibPedestal.h>
+#include <AliFMDRawReader.h>
+#include <iostream>
+#include <fstream>
+#include <TStyle.h>
+#include <TArrayF.h>
+#include <AliLog.h>
+#include <TSystem.h>
+#include <TFile.h>
+
+/** @class DrawDigits
+    @brief Draw hit energy loss versus digit ADC
+    @code 
+    Root> .L Compile.C
+    Root> Compile("DrawCalibRaw.C")
+    Root> DrawCalibRaw cr;
+    Root> cr.SetRawFile("file");
+    Root> cr.Run();
+    @endcode
+    @ingroup FMD_script
+ */
+class DrawCalibRaw : public AliFMDInput
+{
+private:
+  Double_t   fFactor;
+  TString    fCalibDir;
+  TH1D*      fELoss; // Histogram 
+  TObjArray* fDets;
+  TFile*     fOut;
+  Bool_t     fHasData;
+  Int_t      fGotNEvents;
+
+public:
+  //__________________________________________________________________
+  DrawCalibRaw(const       char*  file, 
+              const char* calibDir    = 0,
+              Double_t    noiseFactor = 5, 
+              Bool_t      save        = kTRUE,
+              Int_t       m           = 420, 
+              Double_t    mmin        = -0.5, 
+              Double_t    mmax        = 20.5) 
+    : AliFMDInput(), 
+      fFactor(noiseFactor),
+      fCalibDir(""), 
+      fELoss(0), 
+      fDets(0), 
+      fOut(0), 
+      fHasData(kFALSE),
+      fGotNEvents(0)
+  { 
+    if (calibDir) fCalibDir = calibDir;
+    AddLoad(kRawCalib);
+    fELoss = new TH1D("eLoss", "Scaled Energy loss", m, mmin, mmax);
+    fELoss->SetXTitle("#Delta E/#Delta E_{mip}");
+    fELoss->Sumw2();
+    fELoss->SetDirectory(0);
+    
+    SetRawFile(file);
+
+    if (!save) return;
+    fDets = new TObjArray(3);
+
+  }
+  //__________________________________________________________________
+  Bool_t CheckFile(const char* prefix, int number, TString& f)
+  {
+    f = (Form("%s%d.csv", prefix, number));
+    f = gSystem->Which(fCalibDir.Data(), f.Data());
+    return !f.IsNull();
+  }
+
+  //__________________________________________________________________
+  Bool_t Init()
+  {
+    AliCDBManager* cdb = AliCDBManager::Instance();
+    cdb->SetRun(0);
+    cdb->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+
+    AliFMDCalibStripRange* range = 0;
+    AliFMDCalibSampleRate* rate  = 0;
+    AliFMDCalibPedestal*   peds  = 0;
+    AliFMDCalibGain*       gains = 0;
+    for (Int_t i = 1; i <= 3; i++) { 
+      TString f;
+      if (CheckFile("conditions", i, f)) {
+       Info("Init", "Reading conditions for FMD%d from %s", i, f.Data());
+       if (!range) range = new AliFMDCalibStripRange;
+       if (!rate)  rate  = new AliFMDCalibSampleRate;
+       std::ifstream in(f.Data());
+       range->ReadFromFile(in);
+       rate->ReadFromFile(in);
+      }
+      if (CheckFile("peds", i, f)) {
+       Info("Init", "Reading pedestals for FMD%d from %s", i, f.Data());
+       if (!peds) peds = new AliFMDCalibPedestal;
+       std::ifstream in(f.Data());
+       peds->ReadFromFile(in);
+      }
+      if (CheckFile("gains", i, f)) {
+       Info("Init", "Reading gains for FMD%d from %s", i, f.Data());
+       if (!gains) gains = new AliFMDCalibGain;
+       std::ifstream in(f.Data());
+       gains->ReadFromFile(in);
+      }
+    }
+
+    Int_t mask = (AliFMDParameters::kDeadMap|
+                 AliFMDParameters::kZeroSuppression|
+                 AliFMDParameters::kAltroMap);
+
+    if (!range) mask |= AliFMDParameters::kStripRange;
+    if (!rate)  mask |= AliFMDParameters::kSampleRate;
+    if (!peds)  mask |= AliFMDParameters::kPedestal;
+    if (!gains) mask |= AliFMDParameters::kPulseGain;
+
+    AliFMDParameters* pars = AliFMDParameters::Instance();
+    pars->Init(kFALSE, mask);
+
+    if (range)  pars->SetStripRange(range);
+    if (rate)   pars->SetSampleRate(rate);
+    if (peds)   pars->SetPedestal(peds);
+    if (gains)  pars->SetGain(gains);
+    
+    Bool_t ret = AliFMDInput::Init();
+    
+    if (!fDets) return ret;
+
+    fOut  = TFile::Open(Form("histo_%s", fRawFile.Data()), "RECREATE");
+
+    Int_t m    = fELoss->GetXaxis()->GetNbins();
+    Int_t mmin = fELoss->GetXaxis()->GetXmin();
+    Int_t mmax = fELoss->GetXaxis()->GetXmax();
+    for (Int_t d = 1; d <= 3; d++) {
+      Int_t      nRng = (d == 1 ? 1 : 2);
+      TObjArray* det  = 0;
+      fDets->AddAt(det = new TObjArray(nRng), d-1);
+      det->SetName(Form("FMD%d", d));
+      TDirectory* detD = fOut->mkdir(det->GetName());
+      for (Int_t q = 0; q < nRng; q++) { 
+       Char_t r       = q == 0 ? 'I' : 'O';
+       Int_t  nSec    = q == 0 ?  20 :  40;
+       Int_t  nStr    = q == 0 ? 512 : 256;
+       TObjArray* rng = 0;
+       det->AddAt(rng = new TObjArray(nSec), q);
+       rng->SetName(Form("FMD%d%c", d, r));
+       TDirectory* rngD = detD->mkdir(rng->GetName());
+       for (Int_t s = 0; s < nSec; s++) { 
+         TObjArray* sec = 0;
+         rng->AddAt(sec = new TObjArray(nStr), s);
+         sec->SetName(Form("FMD%d%c_%02d", d, r, s));
+         TDirectory* secD = rngD->mkdir(sec->GetName());
+         for (Int_t t = 0; t < nStr; t++) { 
+           secD->cd();
+           TH1* str = new TH1D(Form("FMD%d%c_%02d_%03d", d, r, s, t), 
+                                Form("Scaled energy loss in FMD%d%c[%2d,%3d]",
+                                     d, r, s, t), m, mmin, mmax);
+           str->SetXTitle("#Delta E/#Delta E_{mip}");
+           // str->SetDirectory(secD);
+           sec->AddAt(str, t);
+         }
+       }
+      }
+    }
+
+    return ret;
+  }
+
+  //__________________________________________________________________
+  Bool_t Begin(Int_t e)
+  {
+    fHasData = kFALSE;
+    return AliFMDInput::Begin(e);
+  }
+
+  //__________________________________________________________________
+  Bool_t ProcessRawCalibDigit(AliFMDDigit* digit)
+  {
+    if (!digit) return kTRUE;
+
+    AliFMDParameters* parm = AliFMDParameters::Instance();
+    UShort_t d             =  digit->Detector();
+    Char_t   r             =  digit->Ring();
+    UShort_t s             =  digit->Sector();
+    UShort_t t             =  digit->Strip();
+    Double_t gain          =  parm->GetPulseGain(d, r, s, t);
+    Double_t ped           =  parm->GetPedestal(d, r, s, t);
+    Double_t pedW          =  parm->GetPedestalWidth(d, r, s, t);
+    Double_t adc           =  digit->Counts();
+    Double_t threshold     =  pedW * fFactor;
+    if (gain < 0.1 || gain > 10) return kTRUE;
+    if (pedW > 10) { 
+      Warning("ProcessRawCalibDigit", "FMD%d%c[%2d,%3d] is noisy: %f",
+             d, r, s, t, pedW);
+      return kTRUE;
+    }
+
+    if (fFMDReader && fFMDReader->IsZeroSuppressed(d-1))
+      adc += fFMDReader->NoiseFactor(d-1) * pedW;
+    else 
+      threshold            += ped;
+
+    if (adc < threshold) return kTRUE;
+    
+    Double_t mult = (adc-ped) / (gain * parm->GetDACPerMIP());
+
+    fHasData = kTRUE;
+    fELoss->Fill(mult);
+
+    // if (t >= 10) return kTRUE;
+    TObjArray* det = static_cast<TObjArray*>(fDets->At(d-1));
+    TObjArray* rng = static_cast<TObjArray*>(det->At(r == 'I' ? 0 : 1));
+    TObjArray* sec = static_cast<TObjArray*>(rng->At(s));
+    TH1*       str = static_cast<TH1*>(sec->At(t));
+    str->Fill(mult);
+
+    return kTRUE;
+  }
+  //__________________________________________________________________
+  Bool_t End()
+  {
+    if (fHasData) fGotNEvents++;
+    return AliFMDInput::End();
+  }
+  //__________________________________________________________________
+  Bool_t Finish()
+  {
+    std::cout << "A total of " << fGotNEvents << " with data" << std::endl;
+    gStyle->SetPalette(1);
+    gStyle->SetOptTitle(0);
+    gStyle->SetCanvasColor(0);
+    gStyle->SetCanvasBorderSize(0);
+    gStyle->SetPadColor(0);
+    gStyle->SetPadBorderSize(0);
+    fELoss->SetStats(kFALSE);
+    fELoss->SetFillColor(kRed);
+    fELoss->SetFillStyle(3001);
+    fELoss->Scale(1. / fELoss->GetEntries());
+    fELoss->DrawCopy("e1 bar");
+
+    if (fDets && fOut) { 
+      std::cout << "Flusing to disk ... " << std::flush;
+      fOut->cd();
+      fELoss->Write();
+      fOut->Write();
+      fOut->Close();
+      std::cout << "done" << std::endl;
+    }
+    return kTRUE;
+  }
+
+  ClassDef(DrawCalibRaw,0);
+};
+
+//____________________________________________________________________
+//
+// EOF
+//