Various small fixes. Make sure Emacs knows it's C++ mode, and the like.
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 8 Mar 2006 14:44:32 +0000 (14:44 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 8 Mar 2006 14:44:32 +0000 (14:44 +0000)
New recconstruction scheme.  Create ESD object (AliESDFMD, added by
Peter to STEER).   Poisson method is not made in the reconstruction
pass.  It should be done one ESD's (the code is still here - will be
and AOD class or similar).

Fixes to the drawers.

Fetch all parameters from AliFMDParameters.  Made ready for CDB access.

Added a TODO file.

44 files changed:
FMD/AliFMD.cxx
FMD/AliFMDAlla.cxx
FMD/AliFMDAlla.h
FMD/AliFMDCalibSampleRate.cxx [new file with mode: 0644]
FMD/AliFMDCalibSampleRate.h [new file with mode: 0644]
FMD/AliFMDDigit.cxx
FMD/AliFMDDigit.h
FMD/AliFMDDigitizer.cxx
FMD/AliFMDDigitizer.h
FMD/AliFMDDigitizerAlla.cxx [new file with mode: 0644]
FMD/AliFMDDigitizerAlla.h [new file with mode: 0644]
FMD/AliFMDFloatMap.cxx [deleted file]
FMD/AliFMDFloatMap.h [deleted file]
FMD/AliFMDInput.cxx
FMD/AliFMDInput.h
FMD/AliFMDMap.cxx [deleted file]
FMD/AliFMDMap.h [deleted file]
FMD/AliFMDMultPoisson.cxx
FMD/AliFMDParameters.cxx
FMD/AliFMDParameters.h
FMD/AliFMDRawReader.cxx
FMD/AliFMDRawWriter.cxx
FMD/AliFMDReconstructor.cxx
FMD/AliFMDReconstructor.h
FMD/AliFMDUShortMap.h
FMD/AliFMDv1.cxx
FMD/Config.C
FMD/FMDbaseLinkDef.h
FMD/FMDrecLinkDef.h
FMD/FMDsimLinkDef.h
FMD/Reconstruct.C
FMD/Simulate.C
FMD/TODO [new file with mode: 0644]
FMD/libFMDbase.pkg
FMD/libFMDrec.pkg
FMD/libFMDsim.pkg
FMD/scripts/CompareESD.C [new file with mode: 0644]
FMD/scripts/Compile.C
FMD/scripts/DrawDigitsRecs.C [new file with mode: 0644]
FMD/scripts/DrawHitPatterns.C [new file with mode: 0644]
FMD/scripts/DrawHits.C
FMD/scripts/DrawHitsDigits.C
FMD/scripts/DrawHitsRecs.C [new file with mode: 0644]
FMD/scripts/TestMapIO.C [new file with mode: 0644]

index 50d9ac8..25e5c7b 100644 (file)
@@ -872,10 +872,10 @@ AliFMD::AddHitByFields(Int_t    track,
        && hit->Sector() == sector 
        && hit->Strip() == strip
        && hit->Track() == track) {
-      Warning("AddHit", "already had a hit in FMD%d%c[%2d,%3d] for track # %d,"
-             " adding energy (%f) to that hit (%f) -> %f", 
-             detector, ring, sector, strip, track, edep, hit->Edep(),
-             hit->Edep() + edep);
+      AliDebug(1, Form("already had a hit in FMD%d%c[%2d,%3d] for track # %d,"
+                      " adding energy (%f) to that hit (%f) -> %f", 
+                      detector, ring, sector, strip, track, edep, hit->Edep(),
+                      hit->Edep() + edep));
       hit->SetEdep(hit->Edep() + edep);
       return hit;
     }
index 1a508e7..3c51d7b 100644 (file)
@@ -45,6 +45,7 @@
 #include <AliFMDHit.h>
 #include "AliFMDv0.h"
 #include "AliFMDAlla.h"
+#include "AliFMDDigitizerAlla.h"
 #include "AliMagF.h"
 #include "AliRun.h"
 #include "AliMC.h"
@@ -409,7 +410,7 @@ AliFMDAlla::StepManager()
   if(id==fIdSens1||id==fIdSens2||id==fIdSens3||id==fIdSens4||id==fIdSens5) {
 
     if(gMC->IsTrackEntering())  {
-      vol[2]=copy;
+      vol[2]=copy-1;
       gMC->CurrentVolOffID(1,copy1);
       vol[1]=copy1;
       gMC->CurrentVolOffID(2,copy2);
@@ -475,12 +476,18 @@ AliFMDAlla::StepManager()
       Char_t   ring     = (vol[0] % 2) == 0 ? 'I' : 'O';
       UShort_t sector   = vol[1];
       UShort_t strip    = vol[2];
+      gMC->TrackPosition(pos);
+      TVector3 cur(pos.Vect());
+      TVector3 old(hits[0], hits[1], hits[2]);
+      cur -= old;
+      Double_t len = cur.Mag();
       AliFMDHit* h = new(lhits[fNhits++]) 
        AliFMDHit(fIshunt,
                  gAlice->GetMCApp()->GetCurrentTrackNumber(),
                  detector, ring, sector, strip, 
                  hits[0], hits[1], hits[2], hits[3], hits[4], hits[5], 
-                 hits[6], iPart, hits[8]);
+                 hits[6], iPart, hits[8], len, 
+                 gMC->IsTrackDisappeared()||gMC->IsTrackStop());
       if (hits[6] > 1 && p/mass > 1) fBad->Add(h);
     } // IsTrackExiting()
   }
@@ -496,6 +503,15 @@ AliFMDAlla::Response(Float_t Edep)
   if (Edep>0)
     charge=Int_t(gRandom->Gaus(chargeOnly,500));       
 }   
+
+//____________________________________________________________________
+AliDigitizer* 
+AliFMDAlla::CreateDigitizer(AliRunDigitizer* manager) const
+{
+  // Create a digitizer object 
+  AliFMDDigitizerAlla* digitizer = new AliFMDDigitizerAlla(manager);
+  return digitizer;
+}
 //--------------------------------------------------------------------------
 //
 // EOF
index fef5458..6b90767 100644 (file)
@@ -24,10 +24,9 @@ public:
   //  virtual void Hit2Digits(Int_t bgrEvent, Option_t *opt1=" ",
   //   Option_t *opt2=" ",Text_t *name=" "); // hit to digit for v1 :test  
  virtual void  Response( Float_t Edep);
-//private:
+ //private:
  //Int_t fCharge; 
-
-
+ AliDigitizer* CreateDigitizer(AliRunDigitizer* manager) const;
 protected:
    Int_t fIdSens1; // Sensetive volume  in FMD
    Int_t fIdSens2; // Sensetive volume  in FMD
@@ -46,5 +45,9 @@ protected:
 };
 
 #endif
-
+//
+// Local Variables:
+//   mode: C++
+// End:
+//
 
diff --git a/FMD/AliFMDCalibSampleRate.cxx b/FMD/AliFMDCalibSampleRate.cxx
new file mode 100644 (file)
index 0000000..de7357f
--- /dev/null
@@ -0,0 +1,70 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+//____________________________________________________________________
+//                                                                          
+//
+//
+#include "AliFMDCalibSampleRate.h"     // ALIFMDCALIBGAIN_H
+#include "AliFMDParameters.h"           // ALIFMDPARAMETERS_H
+
+//____________________________________________________________________
+ClassImp(AliFMDCalibSampleRate)
+#if 0
+  ; // This is here to keep Emacs for indenting the next line
+#endif
+
+//____________________________________________________________________
+AliFMDCalibSampleRate::AliFMDCalibSampleRate()
+  : fRates(3)
+{
+  fRates.Reset(0);
+}
+
+//____________________________________________________________________
+AliFMDCalibSampleRate::AliFMDCalibSampleRate(const AliFMDCalibSampleRate& o)
+  : TObject(o), fRates(o.fRates)
+{}
+
+//____________________________________________________________________
+AliFMDCalibSampleRate&
+AliFMDCalibSampleRate::operator=(const AliFMDCalibSampleRate& o)
+{
+  fRates     = o.fRates;
+  return (*this);
+}
+
+//____________________________________________________________________
+void
+AliFMDCalibSampleRate::Set(UShort_t ddl, UShort_t rate)
+{
+  if (ddl - AliFMDParameters::kBaseDDL < 0) return;
+  fRates[ddl - AliFMDParameters::kBaseDDL] = rate;
+}
+
+//____________________________________________________________________
+UShort_t
+AliFMDCalibSampleRate::Rate(UShort_t ddl) const
+{
+  if (ddl - AliFMDParameters::kBaseDDL < 0) return 0;
+  return fRates[ddl - AliFMDParameters::kBaseDDL];
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
diff --git a/FMD/AliFMDCalibSampleRate.h b/FMD/AliFMDCalibSampleRate.h
new file mode 100644 (file)
index 0000000..0d12399
--- /dev/null
@@ -0,0 +1,39 @@
+#ifndef ALIFMDCALIBSAMPLERATE_H
+#define ALIFMDCALIBSAMPLERATE_H
+/* Copyright(c) 1998-2000, ALICE Experiment at CERN, All rights
+ * reserved. 
+ *
+ * See cxx source for full Copyright notice                               
+ */
+#ifndef ROOT_TObject
+# include <TObject.h>
+#endif
+#ifndef ROOT_TArrayI
+# include <TArrayI.h>
+#endif
+//____________________________________________________________________
+//
+// Gain value and width for each strip in the FMD
+//
+class AliFMDCalibSampleRate : public TObject
+{
+public:
+  AliFMDCalibSampleRate();
+  AliFMDCalibSampleRate(const AliFMDCalibSampleRate& o);
+  AliFMDCalibSampleRate& operator=(const AliFMDCalibSampleRate& o);
+  void Set(UShort_t ddl, UShort_t rate);
+  UShort_t Rate(UShort_t ddl) const;
+protected:
+  TArrayI fRates; // Sample rates 
+  ClassDef(AliFMDCalibSampleRate,1); // Sample rates 
+};
+
+#endif
+//____________________________________________________________________
+//
+// Local Variables:
+//   mode: C++
+// End:
+//
+
+
index e40bfde..c600e6a 100644 (file)
 
 #include "AliFMDDigit.h"       // ALIFMDDIGIT_H
 #include "Riostream.h"         // ROOT_Riostream
+#include <TString.h>
 
 //====================================================================
 ClassImp(AliFMDBaseDigit)
 #if 0
-  ; // This is here to keep Emacs for indenting the next line
+  ; // This is here to keep Emacs from indenting the next line
 #endif
 
 //____________________________________________________________________
@@ -108,6 +109,13 @@ AliFMDBaseDigit::Print(Option_t* /* option*/) const
        << flush;
 }
 
+//____________________________________________________________________
+const char*
+AliFMDBaseDigit::GetName() const 
+{ 
+  return Form("FMD%d%c[%2d,%3d]", fDetector, fRing, fSector, fStrip);
+}
+
 //====================================================================
 ClassImp(AliFMDDigit)
 
index ba2ef9d..bf392fb 100644 (file)
@@ -26,7 +26,7 @@ public:
   UShort_t     Sector()                   const { return fSector;   }
   UShort_t     Strip()            const { return fStrip;    }
   virtual void Print(Option_t* opt="") const;
-
+  const char*  GetName() const;
 protected:
   UShort_t fDetector;  // (Sub) Detector # (1,2, or 3)
   Char_t   fRing;      // Ring ID ('I' or 'O')
index e2ea608..f175c2d 100644 (file)
 #include "AliFMDRing.h"                // ALIFMDRING_H
 #include "AliFMDHit.h"         // ALIFMDHIT_H
 #include "AliFMDDigit.h"       // ALIFMDDIGIT_H
+#include "AliFMDParameters.h"   // ALIFMDPARAMETERS_H
 #include <AliRunDigitizer.h>   // ALIRUNDIGITIZER_H
 #include <AliRun.h>            // ALIRUN_H
 #include <AliLoader.h>         // ALILOADER_H
@@ -231,9 +232,6 @@ AliFMDBaseDigitizer::AliFMDBaseDigitizer(AliRunDigitizer* manager)
 {
   // Normal CTOR
   AliDebug(1," processed");
-  SetVA1MipRange();
-  SetAltroChannelSize();
-  SetSampleRate();
   SetShapingTime();
 }
 
@@ -249,9 +247,6 @@ AliFMDBaseDigitizer::AliFMDBaseDigitizer(const Char_t* name,
 {
   // Normal CTOR
   AliDebug(1," processed");
-  SetVA1MipRange();
-  SetAltroChannelSize();
-  SetSampleRate();
   SetShapingTime();
 }
 
@@ -271,6 +266,16 @@ AliFMDBaseDigitizer::Init()
  
 
 //____________________________________________________________________
+UShort_t
+AliFMDBaseDigitizer::MakePedestal(UShort_t, 
+                                 Char_t, 
+                                 UShort_t, 
+                                 UShort_t) const 
+{ 
+  return 0; 
+}
+
+//____________________________________________________________________
 void
 AliFMDBaseDigitizer::SumContributions(AliFMD* fmd) 
 {
@@ -306,6 +311,7 @@ AliFMDBaseDigitizer::SumContributions(AliFMD* fmd)
   // Get number of entries in the tree 
   Int_t ntracks  = Int_t(hitsTree->GetEntries());
   
+  AliFMDParameters* param = AliFMDParameters::Instance();
   Int_t read = 0;
   // Loop over the tracks in the 
   for (Int_t track = 0; track < ntracks; track++)  {
@@ -325,10 +331,16 @@ AliFMDBaseDigitizer::SumContributions(AliFMD* fmd)
       UShort_t sector   = fmdHit->Sector();
       UShort_t strip    = fmdHit->Strip();
       Float_t  edep     = fmdHit->Edep();
+
+      // Check if strip is `dead' 
+      if (param->IsDead(detector, ring, sector, strip)) continue;
+      
+      // Give warning in case of double hit 
       if (fEdep(detector, ring, sector, strip).fEdep != 0)
        AliDebug(5, Form("Double hit in %d%c(%d,%d)", 
                         detector, ring, sector, strip));
       
+      // Sum energy deposition
       fEdep(detector, ring, sector, strip).fEdep  += edep;
       fEdep(detector, ring, sector, strip).fN     += 1;
       // Add this to the energy deposited for this strip
@@ -354,8 +366,9 @@ AliFMDBaseDigitizer::DigitizeHits(AliFMD* fmd) const
     AliFMDDetector* det = geometry->GetDetector(detector);
     if (!det) continue;
     for (UShort_t ringi = 0; ringi <= 1; ringi++) {
+      Char_t ring = ringi == 0 ? 'I' : 'O';
       // Get pointer to Ring
-      AliFMDRing* r = det->GetRing((ringi == 0 ? 'I' : 'O'));
+      AliFMDRing* r = det->GetRing(ring);
       if (!r) continue;
       
       // Get number of sectors 
@@ -373,20 +386,18 @@ AliFMDBaseDigitizer::DigitizeHits(AliFMD* fmd) const
          // VA1_ALICE channel. 
          if (strip % 128 == 0) last = 0;
          
-         Float_t edep = fEdep(detector, r->GetId(), sector, strip).fEdep;
-         ConvertToCount(edep, last, r->GetSiThickness(), 
-                        geometry->GetSiDensity(), counts);
+         Float_t edep = fEdep(detector, ring, sector, strip).fEdep;
+         ConvertToCount(edep, last, detector, ring, sector, strip, counts);
          last = edep;
-         AddDigit(fmd, detector, r->GetId(), sector, strip, 
-                  edep, UShort_t(counts[0]), 
-                  Short_t(counts[1]), Short_t(counts[2]));
+         AddDigit(fmd, detector, ring, sector, strip, edep, 
+                  UShort_t(counts[0]), Short_t(counts[1]), 
+                  Short_t(counts[2]));
 #if 0
          // This checks if the digit created will give the `right'
          // number of particles when reconstructed, using a naiive
          // approach.  It's here only as a quality check - nothing
          // else. 
-         CheckDigit(fEdep(detector, r->GetId(), sector, strip).fEdep,
-                    fEdep(detector, r->GetId(), sector, strip).fN,
+         CheckDigit(digit, fEdep(detector, ring, sector, strip).fN,
                     counts);
 #endif
        } // Strip
@@ -399,8 +410,10 @@ AliFMDBaseDigitizer::DigitizeHits(AliFMD* fmd) const
 void
 AliFMDBaseDigitizer::ConvertToCount(Float_t   edep, 
                                    Float_t   last,
-                                   Float_t   siThickness, 
-                                   Float_t   siDensity, 
+                                   UShort_t  detector, 
+                                   Char_t    ring, 
+                                   UShort_t  sector, 
+                                   UShort_t  strip,
                                    TArrayI&  counts) const
 {
   // Convert the total energy deposited to a (set of) ADC count(s). 
@@ -440,25 +453,25 @@ AliFMDBaseDigitizer::ConvertToCount(Float_t   edep,
   // 
   //                  = E + (l - E) * ext(-B * t)
   // 
-  Float_t  mipEnergy  = 1.664 * siThickness * siDensity;
-  Float_t  convF      = (1/mipEnergy*Float_t(fAltroChannelSize)/fVA1MipRange);
-  UShort_t ped       = MakePedestal();
-
+  AliFMDParameters* param = AliFMDParameters::Instance();
+  Float_t  convF          = 1/param->GetPulseGain(detector,ring,sector,strip);
+  UShort_t ped            = MakePedestal(detector,ring,sector,strip);
+  UInt_t   maxAdc         = param->GetAltroChannelSize();
+  UShort_t rate           = param->GetSampleRate(AliFMDParameters::kBaseDDL);
+  UShort_t size           = param->GetAltroChannelSize();
+  
   // In case we don't oversample, just return the end value. 
-  if (fSampleRate == 1) {
-    counts[0] = UShort_t(TMath::Min(edep * convF + ped, 
-                                   Float_t(fAltroChannelSize)));
+  if (rate == 1) {
+    counts[0] = UShort_t(TMath::Min(edep * convF + ped, Float_t(size)));
     return;
   }
   
   // Create a pedestal 
-  Int_t   n = fSampleRate;
   Float_t b = fShapingTime;
-  for (Ssiz_t i = 0; i < n;  i++) {
-    Float_t t = Float_t(i) / n;
+  for (Ssiz_t i = 0; i < rate;  i++) {
+    Float_t t = Float_t(i) / rate;
     Float_t s = edep + (last - edep) * TMath::Exp(-b * t);
-    counts[i] = UShort_t(TMath::Min(s * convF + ped, 
-                                   Float_t(fAltroChannelSize))); 
+    counts[i] = UShort_t(TMath::Min(s * convF + ped, Float_t(maxAdc)));
   }
 }
 
@@ -479,7 +492,6 @@ AliFMDDigitizer::AliFMDDigitizer(AliRunDigitizer* manager)
 {
   // Normal CTOR
   AliDebug(1," processed");
-  SetPedestal();
 }
 
 //____________________________________________________________________
@@ -505,7 +517,8 @@ AliFMDDigitizer::Exec(Option_t*)
   if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
 
   // Get the AliFMD object 
-  AliFMD* fmd = static_cast<AliFMD*>(fRunLoader->GetAliRun()->GetDetector("FMD"));
+  AliFMD* fmd = 
+    static_cast<AliFMD*>(fRunLoader->GetAliRun()->GetDetector("FMD"));
   if (!fmd) {
     AliError("Can not get FMD from gAlice");
     return;
@@ -555,10 +568,16 @@ AliFMDDigitizer::Exec(Option_t*)
 
 //____________________________________________________________________
 UShort_t
-AliFMDDigitizer::MakePedestal() const 
+AliFMDDigitizer::MakePedestal(UShort_t  detector, 
+                             Char_t    ring, 
+                             UShort_t  sector, 
+                             UShort_t  strip) const 
 {
   // Make a pedestal 
-  return UShort_t(TMath::Max(gRandom->Gaus(fPedestal, fPedestalWidth), 0.));
+  AliFMDParameters* param =AliFMDParameters::Instance();
+  Float_t           mean  =param->GetPedestal(detector,ring,sector,strip);
+  Float_t           width =param->GetPedestalWidth(detector,ring,sector,strip);
+  return UShort_t(TMath::Max(gRandom->Gaus(mean, width), 0.));
 }
 
 //____________________________________________________________________
@@ -579,18 +598,27 @@ AliFMDDigitizer::AddDigit(AliFMD*  fmd,
 
 //____________________________________________________________________
 void
-AliFMDDigitizer::CheckDigit(Float_t         /* edep */, 
+AliFMDDigitizer::CheckDigit(AliFMDDigit*    digit,
                            UShort_t        nhits,
                            const TArrayI&  counts) 
 {
   // Check that digit is consistent
-  Int_t integral = counts[0];
+  AliFMDParameters* param = AliFMDParameters::Instance();
+  UShort_t          det   = digit->Detector();
+  Char_t            ring  = digit->Ring();
+  UShort_t          sec   = digit->Sector();
+  UShort_t          str   = digit->Strip();
+  Float_t           mean  = param->GetPedestal(det,ring,sec,str);
+  Float_t           width = param->GetPedestalWidth(det,ring,sec,str);
+  UShort_t          range = param->GetVA1MipRange();
+  UShort_t          size  = param->GetAltroChannelSize();
+  Int_t             integral = counts[0];
   if (counts[1] >= 0) integral += counts[1];
   if (counts[2] >= 0) integral += counts[2];
-  integral -= Int_t(fPedestal + 2 * fPedestalWidth);
+  integral -= Int_t(mean + 2 * width);
   if (integral < 0) integral = 0;
   
-  Float_t convF = Float_t(fVA1MipRange) / fAltroChannelSize;
+  Float_t convF = Float_t(range) / size;
   Float_t mips  = integral * convF;
   if (mips > Float_t(nhits) + .5 || mips < Float_t(nhits) - .5) 
     Warning("CheckDigit", "Digit -> %4.2f MIPS != %d +/- .5 hits", 
index 8105066..60f0c92 100644 (file)
@@ -26,6 +26,7 @@ class TClonesArray;
 class AliFMD;
 class AliLoader;
 class AliRunLoader;
+class AliFMDDigit;
 
 
 //====================================================================
@@ -41,24 +42,22 @@ public:
   virtual Bool_t Init();
 
   // Extra member functions 
-  void     SetVA1MipRange(UShort_t r=20) { fVA1MipRange = r; }
-  void     SetAltroChannelSize(UShort_t s=1024) { fAltroChannelSize = s;}
-  void     SetSampleRate(UShort_t r=1) { fSampleRate = (r>2 ? 2 : r); }
-  void     SetShapingTime(Float_t t=10) { fShapingTime = t;  }
-  
-  UShort_t GetVA1MipRange()      const { return fVA1MipRange; }
-  UShort_t GetAltroChannelSize() const { return fAltroChannelSize; }
-  UShort_t GetSampleRate()       const { return fSampleRate; }
+  void     SetShapingTime(Float_t t=10) { fShapingTime = t;  }  
   Float_t  GetShapingTime()      const { return fShapingTime; }
 protected:
   virtual void     SumContributions(AliFMD* fmd);
   virtual void     DigitizeHits(AliFMD* fmd) const;
   virtual void     ConvertToCount(Float_t   edep, 
                                  Float_t   last,
-                                 Float_t   siThickness, 
-                                 Float_t   siDensity, 
+                                 UShort_t  detector, 
+                                 Char_t    ring, 
+                                 UShort_t  sector, 
+                                 UShort_t  strip,
                                  TArrayI&  counts) const;
-  virtual UShort_t MakePedestal() const { return 0; }
+  virtual UShort_t MakePedestal(UShort_t  detector, 
+                               Char_t    ring, 
+                               UShort_t  sector, 
+                               UShort_t  strip) const;
   virtual void     AddNoise(TArrayI&) const {}
   virtual void     AddDigit(AliFMD*  /* fmd      */,
                            UShort_t /* detector */, 
@@ -72,15 +71,12 @@ protected:
 
   AliRunLoader* fRunLoader;       // Run loader
   AliFMDEdepMap fEdep;             // Cache of Energy from hits 
-  UShort_t      fVA1MipRange;      // How many MIPs the pre-amp can do    
-  UShort_t      fAltroChannelSize; // Largest # to store in 1 ADC chan.
-  UShort_t      fSampleRate;       // Times the ALTRO samples pre-amp.
   Float_t       fShapingTime;      // Shaping profile parameter
   
   AliFMDBaseDigitizer(const AliFMDBaseDigitizer& o) 
     : AliDigitizer(o) {}
   AliFMDBaseDigitizer& operator=(const AliFMDBaseDigitizer&) { return *this; }
-  ClassDef(AliFMDBaseDigitizer,0) // Base class for FMD digitizers
+  ClassDef(AliFMDBaseDigitizer,1) // Base class for FMD digitizers
 };
 
 //====================================================================
@@ -91,11 +87,6 @@ public:
   AliFMDDigitizer(AliRunDigitizer * manager);
   virtual ~AliFMDDigitizer() {}
   virtual void  Exec(Option_t* option=0);
-  
-   
-  // Extra member functions 
-  void     SetPedestal(Float_t mean=10, Float_t width=.5);
-  void     GetPedestal(Float_t& mean,   Float_t& width) const;
 protected:
   virtual void     AddDigit(AliFMD*  fmd,
                            UShort_t detector, 
@@ -106,30 +97,15 @@ protected:
                            UShort_t count1, 
                            Short_t  count2, 
                            Short_t  count3) const;
-  virtual UShort_t MakePedestal() const;
-  virtual void     CheckDigit(Float_t         edep, 
+  virtual UShort_t MakePedestal(UShort_t  detector, 
+                               Char_t    ring, 
+                               UShort_t  sector, 
+                               UShort_t  strip) const;
+  virtual void     CheckDigit(AliFMDDigit*    digit,
                              UShort_t        nhits,
                              const TArrayI&  counts);
-  Float_t       fPedestal;         // Mean of pedestal 
-  Float_t       fPedestalWidth;    // Width of pedestal 
-  ClassDef(AliFMDDigitizer,0) // Make Digits from Hits
+  ClassDef(AliFMDDigitizer,1) // Make Digits from Hits
 };
-//____________________________________________________________________
-inline void 
-AliFMDDigitizer::SetPedestal(Float_t mean, Float_t width) 
-{
-  fPedestal      = mean;
-  fPedestalWidth = width;
-}
-
-//____________________________________________________________________
-inline void 
-AliFMDDigitizer::GetPedestal(Float_t& mean, Float_t& width)  const
-{
-  mean  = fPedestal;
-  width = fPedestalWidth;
-}
-
 
 //====================================================================
 class AliFMDSDigitizer : public AliFMDBaseDigitizer 
diff --git a/FMD/AliFMDDigitizerAlla.cxx b/FMD/AliFMDDigitizerAlla.cxx
new file mode 100644 (file)
index 0000000..32745ac
--- /dev/null
@@ -0,0 +1,211 @@
+ /**************************************************************************
+ * Copyright(c) 1998-2000, 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.                  *
+ **************************************************************************/
+
+ //////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  Forward Multiplicity Detector based on Silicon plates                    //
+//  This class contains the procedures simulation ADC  signal for            //
+//  the Forward Multiplicity detector  : hits -> digits                      //
+//  ADC signal consists                                                      //
+//   - number of detector;                                                   //
+//   - number of ring;                                                       //
+//   - number of sector;                                                     //
+//   - ADC signal in this channel                                            //
+//                                                                           //
+ //////////////////////////////////////////////////////////////////////////////
+
+#include <TTree.h> 
+#include <TVector.h>
+#include <TObjArray.h>
+#include <TFile.h>
+#include <TDirectory.h>
+#include <TRandom.h>
+#include "AliLog.h"
+#include "AliFMDDigitizerAlla.h"
+#include "AliFMD.h"
+#include "AliFMDHit.h"
+#include "AliFMDDigit.h"
+#include "AliRunDigitizer.h"
+#include "AliRun.h"
+#include "AliLoader.h"
+#include "AliRunLoader.h"
+#include <stdlib.h>
+#include <Riostream.h>
+
+ClassImp(AliFMDDigitizerAlla)
+#if 0
+  ;
+#endif
+
+//___________________________________________
+AliFMDDigitizerAlla::AliFMDDigitizerAlla()  
+  : AliDigitizer()
+{
+  // Default ctor - don't use it
+}
+
+//___________________________________________
+AliFMDDigitizerAlla::AliFMDDigitizerAlla(AliRunDigitizer* manager) 
+  : AliDigitizer(manager) 
+{
+  // ctor which should be used
+  //  fDebug =0;
+  AliDebug(1," processed");
+}
+
+//------------------------------------------------------------------------
+AliFMDDigitizerAlla::~AliFMDDigitizerAlla()
+{
+  // Destructor
+}
+
+ //------------------------------------------------------------------------
+Bool_t AliFMDDigitizerAlla::Init()
+{
+  // Initialization
+  // cout<<"AliFMDDigitizerAlla::Init"<<endl;
+  return kTRUE;
+}
+
+//---------------------------------------------------------------------
+
+void AliFMDDigitizerAlla::Exec(Option_t * /*option*/)
+{
+
+  // Conver hits to digits:
+  //  - number of detector;
+  //  - number of ring;
+  //  - number of sector;
+  //  - ADC signal in this channel
+  AliDebug(1," start...");
+  Float_t de[3][2][50][520];
+  for (Int_t i=0; i<3; i++)
+    for(Int_t j=0; j<2; j++)
+      for(Int_t k=0; k < 50; k++)
+       for (Int_t l=0; l < 520; l++)
+         de[i][j][k][l]=0;
+  Int_t nstrips[]  = { 512, 256 };
+  Int_t nsectors[] = {  20,  40 };
+  
+  AliRunLoader* outRL = 
+    AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
+  AliLoader* pOutFMD = outRL->GetLoader("FMDLoader");
+
+  AliRunLoader* inRL = 
+    AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0));
+
+  if (!inRL) {
+    AliError("Can not find Run Loader for input stream 0");
+    return;
+  }
+  if (!inRL->GetAliRun()) inRL->LoadgAlice();
+
+  AliFMD * fFMD = static_cast<AliFMD*>(inRL->GetAliRun()->GetDetector("FMD"));
+  if (!fFMD) {
+    AliError("Can not get FMD from gAlice");
+    return;
+  }
+
+  // Loop over files to digitize
+  Int_t nFiles = GetManager()->GetNinputs();
+  for (Int_t inputFile=0; inputFile < nFiles; inputFile++)  {
+    AliDebug(1,Form("Digitizing event number %d",
+                   fManager->GetOutputEventNr()));
+    if (fFMD) {
+       AliRunLoader* inRL = 
+        AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inputFile));
+       AliLoader* pInFMD = inRL->GetLoader("FMDLoader");
+       pInFMD->LoadHits("READ");
+      
+      
+      TTree* tH = pInFMD->TreeH();
+      if (!tH) {
+       pInFMD->LoadHits("read");
+       tH = pInFMD->TreeH();
+      }
+      TBranch* brHits = tH->GetBranch("FMD");
+      if (brHits)  fFMD->SetHitsAddressBranch(brHits);
+      else         AliFatal("Branch FMD hit not found");
+      TClonesArray *fFMDhits = fFMD->Hits ();
+      
+      Int_t ntracks    = tH->GetEntries();
+      for (Int_t track = 0; track < ntracks; track++) {
+       brHits->GetEntry(track);
+       Int_t nhits = fFMDhits->GetEntries ();
+       for (Int_t hit = 0; hit < nhits; hit++) {
+         AliFMDHit* fmdHit = 
+           static_cast<AliFMDHit*>(fFMDhits->UncheckedAt(hit));
+         Int_t detector = fmdHit->Detector();
+         Int_t iring    = fmdHit->Ring() == 'I' ? 0 : 1;
+         Int_t sector   = fmdHit->Sector();
+         Int_t strip    = fmdHit->Strip();
+         de[detector-1][iring][sector][strip] += fmdHit->Edep();
+       }          //hit loop
+      } //track loop
+    } //if FMD
+  }
+
+  // Put noise and make ADC signal
+   Float_t mipI = 1.664 * 0.04 * 2.33 / 22400;     // = 6.923e-6;
+   for (Int_t detector = 1; detector <= 3; detector++){
+     for (Int_t iring = 0; iring < 2; iring++) {
+       if (detector == 1 && iring == 1) continue;
+       char ring = (iring == 0 ? 'I' : 'O');
+       for (Int_t sector = 0; sector < nsectors[iring]; sector++) {
+        for (Int_t strip = 0; strip < nstrips[iring]; strip++) {
+          Int_t signal   = Int_t(de[detector-1][iring][sector][strip] / mipI);
+          Int_t pedestal = Int_t(gRandom->Gaus(500,250));
+          Int_t charge   = signal + pedestal;
+          if(charge <= 500) charge = 500; 
+          //dynamic range from MIP(0.155MeV) to 30MIP(4.65MeV)
+          //1024 ADC channels 
+          Float_t channelWidth = (22400 * 50) / 1024;
+          Int_t   adc          = Int_t(charge / channelWidth);
+          if (adc > 1023)  adc = 1023; 
+          fFMD->AddDigitByFields(detector, ring, sector, strip, adc);
+        } //strip
+       } //sector
+     } //iring
+   } // detector
+   
+   pOutFMD->LoadDigits("update");
+   TTree* treeD = pOutFMD->TreeD();
+   if (!treeD) {
+     pOutFMD->MakeTree("D");
+     treeD = pOutFMD->TreeD();
+   }
+
+   treeD->Reset();
+   TClonesArray* digits = fFMD->Digits();
+   fFMD->MakeBranchInTree(treeD, fFMD->GetName(), &(digits), 4000, 0);
+   if (!treeD->GetBranch("FMD")) AliFatal("No branch for FMD digits");
+   treeD->Fill();
+   pOutFMD->WriteDigits("OVERWRITE");
+   pOutFMD->UnloadHits();
+   pOutFMD->UnloadDigits();
+   fFMD->ResetDigits();
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
+
+
+
+
diff --git a/FMD/AliFMDDigitizerAlla.h b/FMD/AliFMDDigitizerAlla.h
new file mode 100644 (file)
index 0000000..4b50f6c
--- /dev/null
@@ -0,0 +1,42 @@
+#ifndef ALIFMDDIGITIZERALLA_H
+#define ALIFMDDIGITIZERALLA_H
+/* Copyright(c) 1998-2000, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+//
+// Resurection of Alla's old digitizer
+// This is to investigate the changes we've seen in the reconstructed 
+// particle multiplicity 
+//
+#include <AliDigitizer.h>
+#include <AliRunDigitizer.h>
+class TClonesArray;
+
+class AliFMDDigitizerAlla : public AliDigitizer 
+{
+public:
+  AliFMDDigitizerAlla();
+  AliFMDDigitizerAlla(AliRunDigitizer * manager);
+  virtual ~AliFMDDigitizerAlla();
+  virtual Bool_t Init();
+  // Do the main work
+  void  Exec(Option_t* option=0) ;
+  Int_t PutNoise(Int_t charge) {return (Int_t)(gRandom->Gaus(charge,500));}
+  TClonesArray *Digits() const {return fDigits;}
+  TClonesArray *Hits()   const {return fHits;}
+  enum {kBgTag = -1};
+private:
+  TClonesArray *fDigits;               // ! array with digits
+  TClonesArray *fHits;                 // List of hits
+  AliRunDigitizer* GetManager(){return fManager;}
+  ClassDef(AliFMDDigitizerAlla,0)
+};    
+#endif
+//
+// Local Variables:
+//  mode: C++
+// End:
+//
+
+
+
+
diff --git a/FMD/AliFMDFloatMap.cxx b/FMD/AliFMDFloatMap.cxx
deleted file mode 100644 (file)
index 5251a30..0000000
+++ /dev/null
@@ -1,128 +0,0 @@
-/**************************************************************
- * 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.   *
- **************************************************************/
-/* $Id$ */
-//__________________________________________________________
-// 
-// Map of per strip Float_t information
-// 
-// Created Mon Nov  8 12:51:51 2004 by Christian Holm Christensen
-// 
-#include "AliFMDFloatMap.h"    //ALIFMDFLOATMAP_H
-//__________________________________________________________
-ClassImp(AliFMDFloatMap)
-#if 0
-  ; // This is here to keep Emacs for indenting the next line
-#endif
-//__________________________________________________________
-AliFMDFloatMap::AliFMDFloatMap(const AliFMDFloatMap& other)
-  : AliFMDMap(other.fMaxDetectors,
-              other.fMaxRings,
-              other.fMaxSectors,
-              other.fMaxStrips),
-    fData(0)
-{
-  // Copy constructor
-  fData = new Float_t[fMaxDetectors * fMaxRings * fMaxSectors * fMaxStrips];
-  for (size_t i = 0; i < fMaxDetectors * fMaxRings 
-        * fMaxSectors * fMaxStrips; i++)
-    fData[i] = other.fData[i];
-}
-
-//__________________________________________________________
-AliFMDFloatMap::AliFMDFloatMap(size_t maxDet,
-                              size_t maxRing,
-                              size_t maxSec,
-                              size_t maxStr)
-  : AliFMDMap(maxDet, maxRing, maxSec, maxStr),
-    fData(0)
-{
-  // Constructor.
-  // Parameters:
-  //   maxDet  Maximum number of detectors
-  //   maxRing Maximum number of rings per detector
-  //   maxSec  Maximum number of sectors per ring
-  //   maxStr  Maximum number of strips per sector
-  fData = new Float_t[fMaxDetectors * fMaxRings * fMaxSectors * fMaxStrips];
-  Reset();
-}
-
-//__________________________________________________________
-AliFMDFloatMap&
-AliFMDFloatMap::operator=(const AliFMDFloatMap& other)
-{
-  // Assignment operator 
-  fMaxDetectors = other.fMaxDetectors;
-  fMaxRings     = other.fMaxRings;
-  fMaxSectors   = other.fMaxSectors;
-  fMaxStrips    = other.fMaxStrips;
-  if (fData) delete [] fData;
-  fData = new Float_t[fMaxDetectors * fMaxRings * fMaxSectors * fMaxStrips];
-  for (size_t i = 0; i < fMaxDetectors * fMaxRings 
-        * fMaxSectors * fMaxStrips; i++)
-    fData[i] = other.fData[i];
-  return *this;
-}
-
-//__________________________________________________________
-void
-AliFMDFloatMap::Reset(const Float_t& val)
-{
-  // Reset map to val
-  for (size_t i = 0; i < fMaxDetectors * fMaxRings 
-        * fMaxSectors * fMaxStrips; i++)
-    fData[i] = val;
-}
-
-//__________________________________________________________
-Float_t&
-AliFMDFloatMap::operator()(UShort_t det, 
-                          Char_t   ring, 
-                          UShort_t sec, 
-                          UShort_t str)
-{
-  // Get data
-  // Parameters:
-  //   det     Detector #
-  //   ring    Ring ID
-  //   sec     Sector #
-  //   str     Strip #
-  // Returns appropriate data
-  return fData[CalcIndex(det, ring, sec, str)];
-}
-
-//__________________________________________________________
-const Float_t&
-AliFMDFloatMap::operator()(UShort_t det, 
-                          Char_t   ring, 
-                          UShort_t sec, 
-                          UShort_t str) const
-{
-  // Get data
-  // Parameters:
-  //   det     Detector #
-  //   ring    Ring ID
-  //   sec     Sector #
-  //   str     Strip #
-  // Returns appropriate data
-  return fData[CalcIndex(det, ring, sec, str)];
-}
-
-//__________________________________________________________
-// 
-// EOF
-// 
-
diff --git a/FMD/AliFMDFloatMap.h b/FMD/AliFMDFloatMap.h
deleted file mode 100644 (file)
index 3f26353..0000000
+++ /dev/null
@@ -1,46 +0,0 @@
-#ifndef ALIFMDFLOATMAP_H
-#define ALIFMDFLOATMAP_H
-/* Copyright(c) 1998-2000, ALICE Experiment at CERN, All rights
- * reserved. 
- *
- * See cxx source for full Copyright notice                               
- */
-#ifndef ALIFMDMAP_H
-# include <AliFMDMap.h>
-#endif
-//____________________________________________________________________
-//
-// Array of floats indexed by strip identifier.
-//
-class AliFMDFloatMap : public AliFMDMap
-{
-public:
-  AliFMDFloatMap(size_t  maxDet = kMaxDetectors, 
-                size_t  maxRing= kMaxRings, 
-                size_t  maxSec = kMaxSectors, 
-                size_t  maxStr = kMaxStrips);
-  AliFMDFloatMap(const AliFMDFloatMap& o);
-  virtual ~AliFMDFloatMap() { delete [] fData; }
-  AliFMDFloatMap& operator=(const AliFMDFloatMap& o);
-  virtual void Reset(const Float_t& v=Float_t());
-  virtual Float_t& operator()(UShort_t det,
-                             Char_t   ring,
-                             UShort_t sec,
-                             UShort_t str);
-  virtual const Float_t& operator()(UShort_t det,
-                                   Char_t   ring,
-                                   UShort_t sec,
-                                   UShort_t str) const;
-protected:
-  Float_t* fData;
-  ClassDef(AliFMDFloatMap,1) // Map of floats
-};
-
-#endif
-//____________________________________________________________________
-//
-// Local Variables:
-//   mode: C++
-// End:
-//
-
index 3b30fb3..49f6ff4 100644 (file)
 //
 #include "AliFMDInput.h"       // ALIFMDHIT_H
 #include "AliLog.h"            // ALILOG_H
-#include "Riostream.h"         // ROOT_Riostream
-#include <TDatabasePDG.h>
-#include <TMath.h>
-#include <TString.h>
-#include <TParticle.h>
+#include "AliLoader.h"          // ALILOADER_H
+#include "AliRunLoader.h"       // ALIRUNLOADER_H
+#include "AliRun.h"             // ALIRUN_H
+#include "AliStack.h"           // ALISTACK_H
+#include "AliFMD.h"             // ALIFMD_H
 #include "AliFMDHit.h"         // ALIFMDHIT_H
-#include "AliFMDDigit.h"               // ALIFMDDigit_H
-#include "AliFMDMultStrip.h"           // ALIFMDMultStrip_H
-#include "AliFMDMultRegion.h"          // ALIFMDMultRegion_H
+#include "AliFMDDigit.h"       // ALIFMDDigit_H
+#include "AliFMDMultStrip.h"   // ALIFMDMultStrip_H
+#include "AliFMDMultRegion.h"  // ALIFMDMultRegion_H
+#include <TTree.h>              // ROOT_TTree
+#include <TParticle.h>          // ROOT_TParticle
+#include <TString.h>            // ROOT_TString
+#include <TDatabasePDG.h>       // ROOT_TDatabasePDG
+#include <TMath.h>              // ROOT_TMath
+#include <TGeoManager.h>        // ROOT_TGeoManager 
+#include <Riostream.h>         // ROOT_Riostream
 
 //____________________________________________________________________
 ClassImp(AliFMDInput)
@@ -152,6 +159,21 @@ AliFMDInput::Init()
     return kFALSE;
   }
   fTreeE = fLoader->TreeE();
+
+  // Optionally, get the geometry 
+  if (TESTBIT(fTreeMask, kGeometry)) {
+    TString fname(fRun->GetGeometryFileName());
+    if (fname.IsNull()) {
+      Warning("Init", "No file name for the geometry from AliRun");
+      fname = gSystem->DirName(fGAliceFile);
+      fname.Append("/geometry.root");
+    }
+    fGeoManager = TGeoManager::Import(fname.Data());
+    if (!fGeoManager) {
+      Fatal("Init", "No geometry manager found");
+      return kFALSE;
+    }
+  }
   
   fIsInit = kTRUE;
   return fIsInit;
@@ -208,9 +230,9 @@ AliFMDInput::Begin(Int_t event)
     if (fFMDLoader->LoadRecPoints()) return kFALSE;
     fTreeR = fFMDLoader->TreeR();
     if (!fArrayN) fArrayN = new TClonesArray("AliFMDMultStrip");
-    if (!fArrayP) fArrayP = new TClonesArray("AliFMDMultRegion");
-    fTreeR->SetBranchAddress("FMDNaiive",  &fArrayN);
-    fTreeR->SetBranchAddress("FMDPoisson", &fArrayP);
+    // if (!fArrayP) fArrayP = new TClonesArray("AliFMDMultRegion");
+    fTreeR->SetBranchAddress("FMD",  &fArrayN);
+    // fTreeR->SetBranchAddress("FMDPoisson", &fArrayP);
   }
   return kTRUE;
 }
index 95681bf..718cb1c 100644 (file)
 // classes for customized class that do some sort of analysis on the
 // various types of data produced by the FMD. 
 //
-#ifndef ALILOADER_H
-# include <AliLoader.h>
-#endif
-#ifndef ALIRUNLOADER_H
-# include <AliRunLoader.h>
-#endif
-#ifndef ALIRUN_H
-# include <AliRun.h>
-#endif
-#ifndef ALISTACK_H
-# include <AliStack.h>
-#endif
-#ifndef ALIFMD_H
-# include <AliFMD.h>
-#endif
-#ifndef ROOT_TTree
-# include <TTree.h>
-#endif
-#ifndef ROOT_TParticle
-# include <TParticle.h>
-#endif
+#include <TObject.h>
 #ifndef ROOT_TString
 # include <TString.h>
 #endif
+class AliRunLoader;
+class AliLoader;
+class AliStack;
+class AliRun;
+class AliFMD;
+class AliFMDHit;
+class TString;
+class TClonesArray;
+class TTree;
+class TGeoManager;
+class TParticle;
+
 
 //___________________________________________________________________
 class AliFMDInput : public TObject
 {
 public:
   enum ETrees {
-    kHits       = 1, 
-    kKinematics, 
-    kDigits, 
-    kSDigits, 
-    kHeader, 
-    kRecPoints
+    kHits       = 1,  // Hits
+    kKinematics,      // Kinematics (from sim)
+    kDigits,          // Digits
+    kSDigits,         // Summable digits 
+    kHeader,          // Header information 
+    kRecPoints,       // Reconstructed points
+    kGeometry         // Not really a tree 
   };
   AliFMDInput();
   AliFMDInput(const char* gAliceFile);
@@ -82,6 +75,7 @@ protected:
   TClonesArray* fArrayS;     // SDigit info array
   TClonesArray* fArrayN;     // Mult (single) info array
   TClonesArray* fArrayP;     // Mult (region) info array
+  TGeoManager*  fGeoManager; // Geometry manager
   Int_t         fTreeMask;   // Which tree's to load
   Bool_t        fIsInit;
   ClassDef(AliFMDInput,0)  //Hits for detector FMD
diff --git a/FMD/AliFMDMap.cxx b/FMD/AliFMDMap.cxx
deleted file mode 100644 (file)
index 245088e..0000000
+++ /dev/null
@@ -1,101 +0,0 @@
-/**************************************************************************
- * 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.                  *
- **************************************************************************/
-
-/* $Id$ */
-
-//____________________________________________________________________
-//                                                                          
-// Base class for caches of per-strip information.
-// This is used to index a strip. 
-// Data stored depends on derived class. 
-// This class provides some common infra-structure.
-// Derived classes sould define Reset, and operator(). 
-//
-#include "AliFMDMap.h"         // ALIFMDMAP_H
-#include "AliLog.h"
-
-//____________________________________________________________________
-ClassImp(AliFMDMap)
-#if 0
-  ; // This is here to keep Emacs for indenting the next line
-#endif
-
-//____________________________________________________________________
-AliFMDMap::AliFMDMap(size_t maxDet, 
-                    size_t maxRing, 
-                    size_t maxSec, 
-                    size_t maxStr)
-  : fMaxDetectors(maxDet), 
-    fMaxRings(maxRing), 
-    fMaxSectors(maxSec), 
-    fMaxStrips(maxStr)
-{
-  // Construct a map
-  //
-  // Parameters:
-  //     maxDet       Maximum # of detectors
-  //     maxRinf      Maximum # of rings
-  //     maxSec       Maximum # of sectors
-  //     maxStr       Maximum # of strips
-}
-
-//____________________________________________________________________
-Int_t 
-AliFMDMap::CheckIndex(size_t det, Char_t ring, size_t sec, size_t str) const
-{
-  // Check that the index supplied is OK.   Returns true index, or -1
-  // on error. 
-  size_t ringi = (ring == 'I' ||  ring == 'i' ? 0 : 1);
-  size_t idx = 
-    (det + fMaxDetectors * (ringi + fMaxRings * (sec + fMaxSectors * str)));
-  if (idx >= fMaxDetectors * fMaxRings * fMaxSectors * fMaxStrips) 
-    return -1;
-  return idx;
-}
-
-    
-//____________________________________________________________________
-size_t 
-AliFMDMap::CalcIndex(size_t det, Char_t ring, size_t sec, size_t str) const
-{
-  // Calculate index into storage from arguments. 
-  // 
-  // Parameters: 
-  //     det       Detector #
-  //     ring      Ring ID
-  //     sec       Sector # 
-  //     str       Strip # 
-  //
-  // Returns appropriate index into storage 
-  //
-  Int_t idx = CheckIndex(det, ring, sec, str);
-  if (idx < 0) {
-    size_t ringi = (ring == 'I' ||  ring == 'i' ? 0 : 1);
-    AliFatal(Form("CalcIndex", "Index (%d,'%c',%d,%d) out of bounds, "
-                 "in particular the %s index ", 
-                 det, ring, sec, str, 
-                 (det >= fMaxDetectors ? "Detector" : 
-                  (ringi >= fMaxRings ? "Ring" : 
-                   (sec >= fMaxSectors ? "Sector" : "Strip")))));
-    return 0;
-  }
-  return size_t(idx);
-}
-
-
-//___________________________________________________________________
-//
-// EOF
-//
diff --git a/FMD/AliFMDMap.h b/FMD/AliFMDMap.h
deleted file mode 100644 (file)
index 2dba2c5..0000000
+++ /dev/null
@@ -1,196 +0,0 @@
-#ifndef ALIFMDMAP_H
-#define ALIFMDMAP_H
-/* Copyright(c) 1998-2000, ALICE Experiment at CERN, All rights
- * reserved. 
- *
- * See cxx source for full Copyright notice                               
- */
-#ifndef ROOT_TObject
-# include <TObject.h>
-#endif 
-//____________________________________________________________________
-//
-// Base class for caches of per-strip information.
-// This is used to index a strip. 
-// Data stored depends on derived class. 
-//
-class AliFMDMap : public TObject 
-{
-public:
-  enum { 
-    kMaxDetectors = 3, 
-    kMaxRings     = 2, 
-    kMaxSectors   = 40, 
-    kMaxStrips    = 512
-  };
-  AliFMDMap(size_t maxDet = kMaxDetectors, 
-           size_t maxRing= kMaxRings, 
-           size_t maxSec = kMaxSectors, 
-           size_t maxStr = kMaxStrips);
-  virtual ~AliFMDMap() {}
-  Int_t CheckIndex(size_t det, Char_t ring, size_t sec, size_t str) const;
-protected:
-  size_t CalcIndex(size_t det, Char_t ring, size_t sec, size_t str) const;
-  size_t fMaxDetectors;             // Maximum # of detectors
-  size_t fMaxRings;                 // Maximum # of rings
-  size_t fMaxSectors;               // Maximum # of sectors
-  size_t fMaxStrips;                // Maximum # of strips
-  ClassDef(AliFMDMap, 1) // Cache of per strip information
-};
-
-#ifdef  MAY_USE_TEMPLATES
-#ifndef __VECTOR__
-# include <vector>
-#endif 
-//____________________________________________________________________
-//
-// Class template for classes that cache per strip information.  
-// Access to the data goes via 
-//
-//   Type& AliFMDMap<Type>::operator()(size_t detector,
-//                                     Char_t ring, 
-//                                     size_t sector,
-//                                     size_t strip);
-// 
-// (as well as a const version of this member function). 
-// The elements can be reset to the default value by calling 
-// AliFMDMap<Type>::Clear().  This resets the values to `Type()'. 
-//
-template <typename Type> 
-class AliFMDMap : public TObject 
-{
-public:
-  AliFMDMap(size_t maxDet=3, size_t maxRing=2, size_t maxSec=40, 
-           size_t maxStr=512);
-  virtual ~AliFMDMap() {}
-  void Clear(const Type& val=Type());
-  Type& operator()(size_t det, Char_t ring, size_t sec, size_t str);
-  const Type& operator()(size_t det, Char_t ring, size_t sec, size_t str)const;
-private:
-  typedef std::vector<Type> ValueVector; // Type of container
-  ValueVector fValues;                   // Contained values
-  size_t      fMaxDetectors;             // Maximum # of detectors
-  size_t      fMaxRings;                 // Maximum # of rings
-  size_t      fMaxSectors;               // Maximum # of sectors
-  size_t      fMaxStrips;                // Maximum # of strips
-  
-  size_t CalcIndex(size_t det, Char_t ring, size_t sec, size_t str) const;
-  ClassDef(AliFMDMap, 0); // Map of FMD index's to values 
-};
-
-
-//____________________________________________________________________
-template <typename Type>
-inline 
-AliFMDMap<Type>::AliFMDMap(size_t maxDet, 
-                          size_t maxRing, 
-                          size_t maxSec, 
-                          size_t maxStr)
-  : fValues(maxDet * maxRing * maxSec * maxStr), 
-    fMaxDetectors(maxDet), 
-    fMaxRings(maxRing), 
-    fMaxSectors(maxSec), 
-    fMaxStrips(maxStr)
-{
-  // Construct a map
-  //
-  // Parameters:
-  //     maxDet       Maximum # of detectors
-  //     maxRinf      Maximum # of rings
-  //     maxSec       Maximum # of sectors
-  //     maxStr       Maximum # of strips
-}
-
-
-//____________________________________________________________________
-template <typename Type>
-inline size_t 
-AliFMDMap<Type>::CalcIndex(size_t det, Char_t ring, size_t sec, size_t str) const
-{
-  // Calculate index into storage from arguments. 
-  // 
-  // Parameters: 
-  //     det       Detector #
-  //     ring      Ring ID
-  //     sec       Sector # 
-  //     str       Strip # 
-  //
-  // Returns appropriate index into storage 
-  //
-  size_t ringi = (ring == 'I' ||  ring == 'i' ? 0 : 1);
-  size_t idx = 
-    (det + fMaxDetectors * (ringi + fMaxRings * (sec + fMaxSectors * str)));
-  if (idx >= fMaxDetectors * fMaxRings * fMaxSectors * fMaxStrips) {
-    Fatal("CalcIndex", "Index (%d,'%c',%d,%d) out of bounds, "
-         "in particular the %s index", 
-         det, ring, sec, str, 
-         (det >= fMaxDetectors ? "Detector" : 
-          (ringi >= fMaxRings ? "Ring" : 
-           (sec >= fMaxSectors ? "Sector" : "Strip"))));
-    return 0;
-  }
-  return idx;
-}
-
-//____________________________________________________________________
-template <typename Type>
-inline void
-AliFMDMap<Type>::Clear(const Type& val) 
-{
-  // Resets stored values to the default value for that type 
-  for (size_t i = 0; i < fValues.size(); ++i) fValues[i] = val;
-}
-
-//____________________________________________________________________
-template <typename Type>
-inline Type& 
-AliFMDMap<Type>::operator()(size_t det, Char_t ring, size_t sec, size_t str)
-{
-  // Parameters: 
-  //     det       Detector #
-  //     ring      Ring ID
-  //     sec       Sector # 
-  //     str       Strip # 
-  //
-  // Returns data[det][ring][sec][str]
-  return fValues[CalcIndex(det, ring, sec, str)];
-}
-
-//____________________________________________________________________
-template <typename Type>
-inline const Type& 
-AliFMDMap<Type>::operator()(size_t det, 
-                           Char_t ring, 
-                           size_t sec, 
-                           size_t str) const
-{
-  // Parameters: 
-  //     det       Detector #
-  //     ring      Ring ID
-  //     sec       Sector # 
-  //     str       Strip # 
-  //
-  // Returns data[det][ring][sec][str]
-  return fValues[CalcIndex(det, ring, sec, str)];
-}
-
-
-//____________________________________________________________________
-// 
-// Some specialisations 
-//
-typedef AliFMDMap<UShort_t> AliFMDAdcMap;
-typedef AliFMDMap<std::pair<Float_t, UShort_t> > AliFMDEdepMap;
-
-#endif
-#endif 
-//____________________________________________________________________
-//
-// Local Variables:
-//   mode: C++
-// End:
-//
-// EOF
-//
-
-
index ef973d4..d1e50b1 100644 (file)
@@ -61,6 +61,7 @@ void
 AliFMDMultPoisson::PreEvent(TTree* tree, Float_t ipZ) 
 {
   // Reset internal data
+  AliDebug(30, Form("before event with vertex %f", ipZ));
   AliFMDMultAlgorithm::PreEvent(tree, ipZ);
   fCurrentVertexZ = ipZ;
   fEmpty.Reset(kFALSE);
@@ -87,7 +88,9 @@ AliFMDMultPoisson::ProcessDigit(AliFMDDigit*  digit,
   //                    vertex of this event 
   //
   if (!digit) return;
-  if (count < fThreshold) fEmpty(digit->Detector() - 1, 
+  AliDebug(30, Form("Processing digit %s (%s)", digit->GetName(), 
+                   count < fThreshold ? "empty" : "hit"));
+  if (count < fThreshold) fEmpty(digit->Detector(),
                                 digit->Ring(), 
                                 digit->Sector(), 
                                 digit->Strip()) = kTRUE;
@@ -182,7 +185,7 @@ AliFMDMultPoisson::PostEvent()
          Int_t   emptyStrips = 0;
          for (Int_t sector = minSector; sector < maxSector; sector++) 
            for (Int_t strip = minStrip; strip < maxStrip; strip++) 
-             if (fEmpty(sub->GetId() - 1, r->GetId(), sector, strip)) 
+             if (fEmpty(sub->GetId(), r->GetId(), sector, strip)) 
                emptyStrips++;
          
          // The total number of strips 
index b152a3a..40b84c0 100644 (file)
 // Eventually, this class will use the Conditions DB to get the
 // various parameters, which code can then request from here.
 //                                                       
-#include "AliFMDParameters.h"  // ALIFMDPARAMETERS_H
-#include "AliFMDGeometry.h"    // ALIFMDGEOMETRY_H
-#include "AliFMDRing.h"                // ALIFMDRING_H
-#include "AliLog.h"            // ALILOG_H
+#include "AliLog.h"               // ALILOG_H
+#include "AliFMDParameters.h"     // ALIFMDPARAMETERS_H
+#include "AliFMDGeometry.h"       // ALIFMDGEOMETRY_H
+#include "AliFMDRing.h"                   // ALIFMDRING_H
+#include "AliFMDCalibGain.h"       // ALIFMDCALIBGAIN_H
+#include "AliFMDCalibPedestal.h"   // ALIFMDCALIBPEDESTAL_H
+#include "AliFMDCalibSampleRate.h" // ALIFMDCALIBPEDESTAL_H
 #include <Riostream.h>
 
 //====================================================================
@@ -50,7 +53,14 @@ AliFMDParameters::Instance()
 
 //____________________________________________________________________
 AliFMDParameters::AliFMDParameters() 
-  : fSiDeDxMip(1.664)
+  : fSiDeDxMip(1.664), 
+    fFixedPulseGain(0), 
+    fEdepMip(0),
+    fZeroSuppression(0), 
+    fSampleRate(0), 
+    fPedestal(0), 
+    fPulseGain(0), 
+    fDeadMap(0)
 {
   // Default constructor 
   SetVA1MipRange();
@@ -61,19 +71,99 @@ AliFMDParameters::AliFMDParameters()
   SetPedestal();
   SetPedestalWidth();
   SetPedestalFactor();
+  SetThreshold();
 }
 
+//__________________________________________________________________
+Float_t
+AliFMDParameters::GetThreshold() const
+{
+  if (!fPulseGain) return fFixedThreshold;
+  return fPulseGain->Threshold();
+}
+
+//__________________________________________________________________
+Float_t
+AliFMDParameters::GetPulseGain(UShort_t detector, Char_t ring, 
+                              UShort_t sector, UShort_t strip) const
+{
+  // Returns the pulser calibrated gain for strip # strip in sector #
+  // sector or ring id ring of detector # detector. 
+  // 
+  // For simulation, this is normally set to 
+  // 
+  //       VA1_MIP_Range 
+  //    ------------------ * MIP_Energy_Loss
+  //    ALTRO_channel_size
+  // 
+  if (!fPulseGain) { 
+    if (fFixedPulseGain <= 0)
+      fFixedPulseGain = fVA1MipRange * GetEdepMip() / fAltroChannelSize;
+    return fFixedPulseGain;
+  }  
+  return fPulseGain->Value(detector, ring, sector, strip);
+}
+
+//__________________________________________________________________
+Bool_t
+AliFMDParameters::IsDead(UShort_t detector, Char_t ring, 
+                        UShort_t sector, UShort_t strip) const
+{
+  if (!fDeadMap) return kFALSE;
+  return fDeadMap->operator()(detector, ring, sector, strip);
+}
+
+//__________________________________________________________________
+UShort_t
+AliFMDParameters::GetZeroSuppression(UShort_t detector, Char_t ring, 
+                                    UShort_t sector, UShort_t strip) const
+{
+  if (!fZeroSuppression) return fFixedZeroSuppression;
+  // Need to map strip to ALTRO chip. 
+  return fZeroSuppression->operator()(detector, ring, sector, strip/128);
+}
+
+//__________________________________________________________________
+UShort_t
+AliFMDParameters::GetSampleRate(UShort_t ddl) const
+{
+  if (!fSampleRate) return fFixedSampleRate;
+  // Need to map sector to digitizier card. 
+  return fSampleRate->Rate(ddl);
+}
 
+//__________________________________________________________________
+Float_t
+AliFMDParameters::GetPedestal(UShort_t detector, Char_t ring, 
+                             UShort_t sector, UShort_t strip) const
+{
+  if (!fPedestal) return fFixedPedestal;
+  return fPedestal->Value(detector, ring, sector, strip);
+}
+
+//__________________________________________________________________
+Float_t
+AliFMDParameters::GetPedestalWidth(UShort_t detector, Char_t ring, 
+                                  UShort_t sector, UShort_t strip) const
+{
+  if (!fPedestal) return fFixedPedestalWidth;
+  return fPedestal->Width(detector, ring, sector, strip);
+}
+  
+                             
 
 //__________________________________________________________________
 Float_t
 AliFMDParameters::GetEdepMip() const 
 { 
   // Get energy deposited by a MIP in the silicon sensors
-  AliFMDGeometry* fmd = AliFMDGeometry::Instance();
-  return (fSiDeDxMip 
-         * fmd->GetRing('I')->GetSiThickness() 
-         * fmd->GetSiDensity());
+  if (fEdepMip <= 0){
+    AliFMDGeometry* fmd = AliFMDGeometry::Instance();
+    fEdepMip = (fSiDeDxMip 
+               * fmd->GetRing('I')->GetSiThickness() 
+               * fmd->GetSiDensity());
+  }
+  return fEdepMip;
 }
 
 //____________________________________________________________________
index 496d5af..754c194 100644 (file)
 #ifndef ROOT_TNamed
 # include <TNamed.h>
 #endif
+#ifndef ROOT_TArrayI
+# include <TArrayI.h>
+#endif
+#ifndef ALIFMDUSHORTMAP_H
+# include <AliFMDUShortMap.h>
+#endif
+#ifndef ALIFMDBOOLMAP_H
+# include <AliFMDBoolMap.h>
+#endif
+typedef AliFMDUShortMap AliFMDCalibZeroSuppression;
+typedef AliFMDBoolMap   AliFMDCalibDeadMap;
+class AliFMDCalibPedestal;
+class AliFMDCalibGain;
+class AliFMDCalibSampleRate;
 
 class AliFMDParameters : public TNamed
 {
 public:
   static AliFMDParameters* Instance();
   
-  // Set various parameters 
-  void     SetVA1MipRange(UShort_t r=20)          { fVA1MipRange = r; }
-  void     SetAltroChannelSize(UShort_t s=1024)   { fAltroChannelSize = s;}
-  void     SetChannelsPerAltro(UShort_t size=128) { fChannelsPerAltro = size; }
-  void     SetZeroSuppression(UShort_t s=0)       { fZeroSuppression = s; }
-  void     SetSampleRate(UShort_t r=1)            { fSampleRate = (r>2?2:r);}
-  void     SetPedestal(Float_t p=10)              { fPedestal = p; }
-  void     SetPedestalWidth(Float_t w=1)          { fPedestalWidth = w; }
-  void     SetPedestalFactor(Float_t f=3)         { fPedestalFactor = f; }
+  // Set various `Fixed' parameters 
+  void SetVA1MipRange(UShort_t r=20)          { fVA1MipRange = r; }
+  void SetAltroChannelSize(UShort_t s=1024)   { fAltroChannelSize = s;}
+  void SetChannelsPerAltro(UShort_t size=128) { fChannelsPerAltro = size; }
+  void SetPedestalFactor(Float_t f=3)         { fPedestalFactor = f; }
+
+  // Set various variable parameter defaults
+  void SetZeroSuppression(UShort_t s=0)       { fFixedZeroSuppression = s; }
+  void SetSampleRate(UShort_t r=1)            { fFixedSampleRate = (r>2?2:r);}
+  void SetPedestal(Float_t p=10)              { fFixedPedestal = p; }
+  void SetPedestalWidth(Float_t w=1)          { fFixedPedestalWidth = w; }
+  void SetThreshold(Float_t t=0)              { fFixedThreshold = t; }
 
-  // Get various parameters
+  // Get `Fixed' various parameters
   UShort_t GetVA1MipRange()          const { return fVA1MipRange; }
   UShort_t GetAltroChannelSize()     const { return fAltroChannelSize; }
   UShort_t GetChannelsPerAltro()     const { return fChannelsPerAltro; }
-  UShort_t GetZeroSuppression()      const { return fZeroSuppression; }
-  UShort_t GetSampleRate()           const { return fSampleRate; }
   Float_t  GetEdepMip()              const;
-  Float_t  GetPedestal()            const { return fPedestal; }
-  Float_t  GetPedestalWidth()       const { return fPedestalWidth; }
   Float_t  GetPedestalFactor()      const { return fPedestalFactor; }
 
+  // Get variable parameters 
+  Bool_t   IsDead(UShort_t detector, 
+                 Char_t ring, 
+                 UShort_t sector, 
+                 UShort_t strip) const;
+  Float_t  GetThreshold() const;
+  Float_t  GetPulseGain(UShort_t detector, 
+                       Char_t ring, 
+                       UShort_t sector, 
+                       UShort_t strip) const;
+  Float_t  GetPedestal(UShort_t detector, 
+                      Char_t ring, 
+                      UShort_t sector, 
+                      UShort_t strip) const;
+  Float_t  GetPedestalWidth(UShort_t detector, 
+                           Char_t ring, 
+                           UShort_t sector, 
+                           UShort_t strip) const;
+  UShort_t GetZeroSuppression(UShort_t detector, 
+                             Char_t ring, 
+                             UShort_t sector, 
+                             UShort_t strip) const;
+  UShort_t GetSampleRate(UShort_t ddl) const;
+
   enum { 
     kBaseDDL = 0x1000 // DDL offset for the FMD
   };
@@ -52,18 +89,28 @@ protected:
   virtual ~AliFMDParameters() {}
   static AliFMDParameters* fgInstance; // Static singleton instance
   
-  const Float_t fSiDeDxMip;        // MIP dE/dx in Silicon
-  UShort_t      fVA1MipRange;      // # MIPs the pre-amp can do    
-  UShort_t      fAltroChannelSize; // Largest # to store in 1 ADC ch.
-  UShort_t      fChannelsPerAltro; // Number of pre-amp. channels/adc chan.
-  UShort_t      fZeroSuppression;  // Threshold for zero-suppression
-  UShort_t      fSampleRate;       // Times the ALTRO samples pre-amp.
-  Float_t       fPedestal;         // Pedestal to subtract
-  Float_t       fPedestalWidth;    // Width of pedestal
-  Float_t       fPedestalFactor;   // Number of pedestal widths
+  const Float_t   fSiDeDxMip;            // MIP dE/dx in Silicon
+  UShort_t        fVA1MipRange;          // # MIPs the pre-amp can do    
+  UShort_t        fAltroChannelSize;     // Largest # to store in 1 ADC ch.
+  UShort_t        fChannelsPerAltro;     // Number of pre-amp. chan/adc chan.
+  Float_t         fPedestalFactor;       // Number of pedestal widths
 
+  Float_t         fFixedPedestal;        // Pedestal to subtract
+  Float_t         fFixedPedestalWidth;   // Width of pedestal
+  UShort_t        fFixedZeroSuppression; // Threshold for zero-suppression
+  UShort_t        fFixedSampleRate;      // Times the ALTRO samples pre-amp.
+  Float_t         fFixedThreshold;       //
+  mutable Float_t fFixedPulseGain;       //! Gain (cached)
+  mutable Float_t fEdepMip;              //! Cache of energy loss for a MIP
+  
+  AliFMDCalibZeroSuppression* fZeroSuppression; // Zero suppression from CDB
+  AliFMDCalibSampleRate*      fSampleRate;      // Sample rate from CDB 
+  AliFMDCalibPedestal*        fPedestal;        // Pedestals 
+  AliFMDCalibGain*            fPulseGain;       // Pulser gain
+  AliFMDCalibDeadMap*         fDeadMap;         // Pulser gain
+  
   
-  ClassDef(AliFMDParameters,1)
+  ClassDef(AliFMDParameters,2)
 };
 
 #endif
index 6f516a1..043f74d 100644 (file)
@@ -66,8 +66,6 @@ AliFMDRawReader::AliFMDRawReader(AliRawReader* reader, TTree* tree)
     fSampleRate(1)
 {
   // Default CTOR
-  AliFMDParameters* pars = AliFMDParameters::Instance();
-  fSampleRate = pars->GetSampleRate();
 }
 
 
@@ -85,11 +83,15 @@ AliFMDRawReader::Exec(Option_t*)
   TClonesArray* array = new TClonesArray("AliFMDDigit");
   fTree->Branch("FMD", &array);
 
+  // Get sample rate 
+  AliFMDParameters* pars = AliFMDParameters::Instance();
+  fSampleRate = pars->GetSampleRate(AliFMDParameters::kBaseDDL);
+
   // Use AliAltroRawStream to read the ALTRO format.  No need to
   // reinvent the wheel :-) 
   AliFMDRawStream input(fReader, fSampleRate);
   // Select FMD DDL's 
-  fReader->Select(AliFMDParameters::kBaseDDL >> 8);
+  fReader->Select(AliFMDParameters::kBaseDDL);
   
   Int_t    oldDDL      = -1;
   Int_t    count       = 0;
index f03e8df..649a17f 100644 (file)
@@ -55,8 +55,7 @@ AliFMDRawWriter::AliFMDRawWriter(AliFMD* fmd)
     fFMD(fmd)
 {
   AliFMDParameters* pars = AliFMDParameters::Instance();
-  fSampleRate            = pars->GetSampleRate();
-  fThreshold             = pars->GetZeroSuppression();
+  fSampleRate            = pars->GetSampleRate(AliFMDParameters::kBaseDDL);
   fChannelsPerAltro      = pars->GetChannelsPerAltro();
 }
 
@@ -137,6 +136,7 @@ AliFMDRawWriter::Exec(Option_t*)
   }
   digitBranch->SetAddress(&digits);
   
+  AliFMDParameters* pars = AliFMDParameters::Instance();
   Int_t nEvents = Int_t(digitTree->GetEntries());
   for (Int_t event = 0; event < nEvents; event++) {
     fFMD->ResetDigits();
@@ -179,6 +179,7 @@ AliFMDRawWriter::Exec(Option_t*)
       Char_t   ring   = digit->Ring();
       UShort_t sector = digit->Sector();
       UShort_t strip  = digit->Strip();
+      fThreshold      = pars->GetZeroSuppression(det, ring, sector, strip);
       if (det != prevDetector) {
        AliDebug(15, Form("FMD: New DDL, was %d, now %d",
                          AliFMDParameters::kBaseDDL + prevDetector - 1,
index b5e23f5..935a2f4 100644 (file)
 #include "AliFMDReconstructor.h"           // ALIFMDRECONSTRUCTOR_H
 #include "AliFMDRawStream.h"               // ALIFMDRAWSTREAM_H
 #include "AliFMDRawReader.h"               // ALIFMDRAWREADER_H
-#include "AliFMDMultAlgorithm.h"          // ALIFMDMULTALGORITHM_H
-#include "AliFMDMultPoisson.h"            // ALIFMDMULTPOISSON_H
-#include "AliFMDMultNaiive.h"             // ALIFMDMULTNAIIVE_H
+#include "AliFMDMultStrip.h"              // ALIFMDMULTNAIIVE_H
 #include "AliESD.h"                       // ALIESD_H
+#include <AliESDFMD.h>                    // ALIESDFMD_H
+#include <TFile.h>
 
 //____________________________________________________________________
 ClassImp(AliFMDReconstructor)
@@ -117,45 +117,20 @@ ClassImp(AliFMDReconstructor)
 //____________________________________________________________________
 AliFMDReconstructor::AliFMDReconstructor() 
   : AliReconstructor(),
-    fPedestal(0), 
-    fPedestalWidth(0),
-    fPedestalFactor(0)
+    fMult(0), 
+    fESDObj(0)
 {
-  // Make a new FMD reconstructor object - default CTOR.
-  AliFMDParameters* pars = AliFMDParameters::Instance();
-  fPedestal       = pars->GetPedestal();
-  fPedestalWidth  = pars->GetPedestalWidth();
-  fPedestalFactor = pars->GetPedestalFactor();
-  
-  fAlgorithms.Add(new AliFMDMultNaiive);
-  fAlgorithms.Add(new AliFMDMultPoisson);
-
-  TIter next(&fAlgorithms);
-  AliFMDMultAlgorithm* algorithm = 0;
-  while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
-    algorithm->PreRun(0);
+  // Make a new FMD reconstructor object - default CTOR.  
 }
   
 
 //____________________________________________________________________
 AliFMDReconstructor::AliFMDReconstructor(const AliFMDReconstructor& other) 
-  : AliReconstructor(),
-    fPedestal(0), 
-    fPedestalWidth(0),
-    fPedestalFactor(0)
+  : AliReconstructor(), 
+    fMult(other.fMult),
+    fESDObj(other.fESDObj)
 {
   // Copy constructor 
-  fPedestal       = other.fPedestal;
-  fPedestalWidth  = other.fPedestalWidth;
-  fPedestalFactor = other.fPedestalFactor;
-
-  
-  fAlgorithms.Delete();
-  TIter next(&(other.fAlgorithms));
-  AliFMDMultAlgorithm* algorithm = 0;
-  while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
-    fAlgorithms.Add(algorithm);
-  fAlgorithms.SetOwner(kFALSE);
 }
   
 
@@ -164,17 +139,8 @@ AliFMDReconstructor&
 AliFMDReconstructor::operator=(const AliFMDReconstructor& other) 
 {
   // Assignment operator
-  fPedestal       = other.fPedestal;
-  fPedestalWidth  = other.fPedestalWidth;
-  fPedestalFactor = other.fPedestalFactor;
-
-  fAlgorithms.Delete();
-  TIter next(&(other.fAlgorithms));
-  AliFMDMultAlgorithm* algorithm = 0;
-  while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
-    fAlgorithms.Add(algorithm);
-  fAlgorithms.SetOwner(kFALSE);
-
+  fMult   = other.fMult;
+  fESDObj = other.fESDObj;
   return *this;
 }
 
@@ -182,7 +148,9 @@ AliFMDReconstructor::operator=(const AliFMDReconstructor& other)
 AliFMDReconstructor::~AliFMDReconstructor() 
 {
   // Destructor 
-  fAlgorithms.Delete();
+  if (fMult)   fMult->Delete();
+  if (fMult)   delete fMult;
+  if (fESDObj) delete fESDObj;
 }
 
 //____________________________________________________________________
@@ -191,16 +159,27 @@ AliFMDReconstructor::Init(AliRunLoader* runLoader)
 {
   // Initialize the reconstructor 
   AliDebug(1, Form("Init called with runloader 0x%x", runLoader));
+  // Current vertex position
   fCurrentVertex = 0;
+  // Create array of reconstructed strip multiplicities 
+  fMult = new TClonesArray("AliFMDMultStrip", 51200);
+  // Create ESD output object 
+  fESDObj = new AliESDFMD;
+  
+  // Check that we have a run loader
   if (!runLoader) { 
     Warning("Init", "No run loader");
     return;
   }
+
+  // Check if we can get the header tree 
   AliHeader* header = runLoader->GetHeader();
   if (!header) {
     Warning("Init", "No header");
     return;
   }
+
+  // Check if we can get a simulation header 
   AliGenEventHeader* eventHeader = header->GenEventHeader();
   if (eventHeader) {
     TArrayF vtx;
@@ -224,6 +203,7 @@ AliFMDReconstructor::ConvertDigits(AliRawReader* reader,
                                   TTree* digitsTree) const
 {
   // Convert Raw digits to AliFMDDigit's in a tree 
+  AliDebug(1, "Reading raw data into digits tree");
   AliFMDRawReader rawRead(reader, digitsTree);
   // rawRead.SetSampleRate(fFMD->GetSampleRate());
   rawRead.Exec();
@@ -241,11 +221,11 @@ AliFMDReconstructor::Reconstruct(TTree* digitsTree,
   AliDebug(1, "Reconstructing from digits in a tree");
 #if 0
   if (fESD) {
-    // const AliESDVertex* vertex = fESD->GetVertex();
-    // if (vertex) {
-    //   AliDebug(1, Form("Got vertex from ESD: %f", vertex->GetZv()));
-    //   fCurrentVertex = vertex->GetZv();
-    // }
+    const AliESDVertex* vertex = fESD->GetVertex();
+    if (vertex) {
+      AliDebug(1, Form("Got vertex from ESD: %f", vertex->GetZv()));
+      fCurrentVertex = vertex->GetZv();
+    }
   }
 #endif  
   TBranch *digitBranch = digitsTree->GetBranch("FMD");
@@ -256,23 +236,40 @@ AliFMDReconstructor::Reconstruct(TTree* digitsTree,
   TClonesArray* digits = new TClonesArray("AliFMDDigit");
   digitBranch->SetAddress(&digits);
 
-  TIter next(&fAlgorithms);
-  AliFMDMultAlgorithm* algorithm = 0;
-  while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
-    algorithm->PreEvent(clusterTree, fCurrentVertex);
+  // TIter next(&fAlgorithms);
+  // AliFMDMultAlgorithm* algorithm = 0;
+  // while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) {
+  //   AliDebug(10, Form("PreEvent called for algorithm %s", 
+  //                     algorithm->GetName()));    
+  //   algorithm->PreEvent(clusterTree, fCurrentVertex);
+  // }
+  if (fMult)   fMult->Clear();
+  if (fESDObj) fESDObj->Clear();
+  
+  fNMult = 0;
+  fTreeR = clusterTree;
+  fTreeR->Branch("FMD", &fMult);
+  
+  AliDebug(5, "Getting entry 0 from digit branch");
   digitBranch->GetEntry(0);
   
+  AliDebug(5, "Processing digits");
   ProcessDigits(digits);
 
-  next.Reset();
-  algorithm = 0;
-  while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
-    algorithm->PostEvent();
-  clusterTree->Fill();
+  // next.Reset();
+  // algorithm = 0;
+  // while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) {
+  //   AliDebug(10, Form("PostEvent called for algorithm %s", 
+  //                     algorithm->GetName()));
+  // algorithm->PostEvent();
+  // }
+  Int_t written = clusterTree->Fill();
+  AliDebug(10, Form("Filled %d bytes into cluster tree", written));
   digits->Delete();
   delete digits;
 }
  
+
 //____________________________________________________________________
 void
 AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
@@ -284,38 +281,49 @@ AliFMDReconstructor::ProcessDigits(TClonesArray* digits) const
   AliDebug(1, Form("Got %d digits", nDigits));
   for (Int_t i = 0; i < nDigits; i++) {
     AliFMDDigit* digit = static_cast<AliFMDDigit*>(digits->At(i));
-    AliFMDGeometry* fmd = AliFMDGeometry::Instance();
-    AliFMDDetector* subDetector = fmd->GetDetector(digit->Detector());
-    if (!subDetector) { 
-      Warning("ProcessDigits", "Unknown detector: FMD%d" , digit->Detector());
-      continue;
-    }
-    
-    AliFMDRing* ring  = subDetector->GetRing(digit->Ring());
-    Float_t     ringZ = subDetector->GetRingZ(digit->Ring());
-    if (!ring) {
-      Warning("ProcessDigits", "Unknown ring: FMD%d%c", digit->Detector(), 
-             digit->Ring());
-      break;
-    }
+    AliFMDParameters* param  = AliFMDParameters::Instance();
+    // Check that the strip is not marked as dead 
+    if (param->IsDead(digit->Detector(), digit->Ring(), 
+                     digit->Sector(), digit->Strip())) continue;
+
+    // digit->Print();
+    // Get eta and phi 
+    Float_t eta, phi;
+    PhysicalCoordinates(digit, eta, phi);
     
-    Float_t  realZ    = fCurrentVertex + ringZ;
-    Float_t  stripR   = ((ring->GetHighR() - ring->GetLowR()) 
-                        / ring->GetNStrips() * (digit->Strip() + .5) 
-                        + ring->GetLowR());
-    Float_t  theta    = TMath::ATan2(stripR, realZ);
-    Float_t  phi      = (2 * TMath::Pi() / ring->GetNSectors() 
-                        * (digit->Sector() + .5));
-    Float_t  eta      = -TMath::Log(TMath::Tan(theta / 2));
+    // Substract pedestal. 
     UShort_t counts   = SubtractPedestal(digit);
     
-    TIter next(&fAlgorithms);
-    AliFMDMultAlgorithm* algorithm = 0;
-    while ((algorithm = static_cast<AliFMDMultAlgorithm*>(next()))) 
-      algorithm->ProcessDigit(digit, eta, phi, counts);
+    // Gain match digits. 
+    Double_t edep     = Adc2Energy(digit, eta, counts);
+    
+    // Make rough multiplicity 
+    Double_t mult     = Energy2Multiplicity(digit, edep);
+    
+    AliDebug(10, Form("FMD%d%c[%2d,%3d]: "
+                     "ADC: %d, Counts: %d, Energy: %f, Mult: %f", 
+                     digit->Detector(), digit->Ring(), digit->Sector(),
+                     digit->Strip(), digit->Counts(), counts, edep, mult));
+    
+    // Create a `RecPoint' on the output branch. 
+    AliFMDMultStrip* m = 
+      new ((*fMult)[fNMult]) AliFMDMultStrip(digit->Detector(), 
+                                            digit->Ring(), 
+                                            digit->Sector(),
+                                            digit->Strip(),
+                                            eta, phi, 
+                                            edep, mult,
+                                            AliFMDMult::kNaiive);
+    (void)m; // Suppress warnings about unused variables. 
+    fNMult++;
+
+    fESDObj->SetMultiplicity(digit->Detector(), digit->Ring(), 
+                            digit->Sector(),  digit->Strip(), mult);
+    fESDObj->SetEta(digit->Detector(), digit->Ring(), 
+                   digit->Sector(),  digit->Strip(), eta);
   }
 }
-      
+
 //____________________________________________________________________
 UShort_t
 AliFMDReconstructor::SubtractPedestal(AliFMDDigit* digit) const
@@ -325,24 +333,185 @@ AliFMDReconstructor::SubtractPedestal(AliFMDDigit* digit) const
   // load this to subtract a pedestal that was given in a database or
   // something like that. 
 
-  Int_t counts = 0;
-  Float_t ped = fPedestal + fPedestalFactor * fPedestalWidth;
+  Int_t             counts = 0;
+  AliFMDParameters* param  = AliFMDParameters::Instance();
+  Float_t           pedM   = param->GetPedestal(digit->Detector(), 
+                                               digit->Ring(), 
+                                               digit->Sector(), 
+                                               digit->Strip());
+  AliDebug(10, Form("Subtracting pedestal %f from signal %d", 
+                  pedM, digit->Counts()));
   if (digit->Count3() > 0)      counts = digit->Count3();
   else if (digit->Count2() > 0) counts = digit->Count2();
   else                          counts = digit->Count1();
-  counts = TMath::Max(Int_t(counts - ped), 0);
+  counts = TMath::Max(Int_t(counts - pedM), 0);
+  if (counts > 0) AliDebug(10, "Got a hit strip");
+  
   return  UShort_t(counts);
 }
 
 //____________________________________________________________________
+Float_t
+AliFMDReconstructor::Adc2Energy(AliFMDDigit* digit, 
+                               Float_t      /* eta */, 
+                               UShort_t     count) const
+{
+  // Converts number of ADC counts to energy deposited. 
+  // Note, that this member function can be overloaded by derived
+  // classes to do strip-specific look-ups in databases or the like,
+  // to find the proper gain for a strip. 
+  // 
+  // In this simple version, we calculate the energy deposited as 
+  // 
+  //    EnergyDeposited = cos(theta) * gain * count
+  // 
+  // where 
+  // 
+  //           Pre_amp_MIP_Range
+  //    gain = ----------------- * Energy_deposited_per_MIP
+  //           ADC_channel_size    
+  // 
+  // is constant and the same for all strips.
+
+  // Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
+  // Double_t edep  = TMath::Abs(TMath::Cos(theta)) * fGain * count;
+  AliFMDParameters* param = AliFMDParameters::Instance();
+  Float_t           gain  = param->GetPulseGain(digit->Detector(), 
+                                               digit->Ring(), 
+                                               digit->Sector(), 
+                                               digit->Strip());
+  Double_t          edep  = count * gain;
+  AliDebug(10, Form("Converting counts %d to energy via factor %f", 
+                   count, gain));
+  return edep;
+}
+
+//____________________________________________________________________
+Float_t
+AliFMDReconstructor::Energy2Multiplicity(AliFMDDigit* /* digit */, 
+                                        Float_t      edep) const
+{
+  // Converts an energy signal to number of particles. 
+  // Note, that this member function can be overloaded by derived
+  // classes to do strip-specific look-ups in databases or the like,
+  // to find the proper gain for a strip. 
+  // 
+  // In this simple version, we calculate the multiplicity as 
+  // 
+  //   multiplicity = Energy_deposited / Energy_deposited_per_MIP
+  // 
+  // where 
+  //
+  //   Energy_deposited_per_MIP = 1.664 * SI_density * SI_thickness 
+  // 
+  // is constant and the same for all strips 
+  AliFMDParameters* param   = AliFMDParameters::Instance();
+  Double_t          edepMIP = param->GetEdepMip();
+  Float_t           mult    = edep / edepMIP;
+  if (edep > 0) 
+    AliDebug(10, Form("Translating energy %f to multiplicity via "
+                    "divider %f->%f", edep, edepMIP, mult));
+  return mult;
+}
+
+//____________________________________________________________________
+void
+AliFMDReconstructor::PhysicalCoordinates(AliFMDDigit* digit, 
+                                        Float_t& eta, 
+                                        Float_t& phi) const
+{
+  // Get the eta and phi of a digit 
+  // 
+  // Get geometry. 
+  AliFMDGeometry* fmd = AliFMDGeometry::Instance();
+  AliFMDDetector* subDetector = fmd->GetDetector(digit->Detector());
+  if (!subDetector) { 
+    Warning("ProcessDigits", "Unknown detector: FMD%d" , digit->Detector());
+    return;
+  }
+    
+  // Get the ring - we should really use geometry->Detector2XYZ(...)
+  // here. 
+  AliFMDRing* ring  = subDetector->GetRing(digit->Ring());
+  Float_t     ringZ = subDetector->GetRingZ(digit->Ring());
+  if (!ring) {
+    Warning("ProcessDigits", "Unknown ring: FMD%d%c", digit->Detector(), 
+           digit->Ring());
+    return;
+  }
+    
+  // Correct for vertex offset. 
+  Float_t  realZ    = fCurrentVertex + ringZ;
+  Float_t  stripR   = ((ring->GetHighR() - ring->GetLowR()) 
+                      / ring->GetNStrips() * (digit->Strip() + .5) 
+                      + ring->GetLowR());
+  Float_t  theta    = TMath::ATan2(stripR, realZ);
+  phi               = (2 * TMath::Pi() / ring->GetNSectors() 
+                      * (digit->Sector() + .5));
+  eta               = -TMath::Log(TMath::Tan(theta / 2));
+}
+
+      
+
+//____________________________________________________________________
 void 
 AliFMDReconstructor::FillESD(TTree*  /* digitsTree */, 
                             TTree*  /* clusterTree */,
-                            AliESD* /* esd*/) const
+                            AliESD* esd) const
 {
   // nothing to be done
   // FIXME: The vertex may not be known when Reconstruct is executed,
   // so we may have to move some of that member function here. 
+  AliDebug(1, Form("Calling FillESD with two trees and one ESD"));
+  // fESDObj->Print();
+
+  if (esd) { 
+    AliDebug(1, Form("Writing FMD data to ESD tree"));
+    esd->SetFMDData(fESDObj);
+    // Let's check the data in the ESD
+#if 0
+    AliESDFMD* fromEsd = esd->GetFMDData();
+    if (!fromEsd) {
+      AliWarning("No FMD object in ESD!");
+      return;
+    }
+    for (UShort_t det = 1; det <= fESDObj->MaxDetectors(); det++) {
+      for (UShort_t ir = 0; ir < fESDObj->MaxRings(); ir++) {
+       Char_t ring = (ir == 0 ? 'I' : 'O');
+       for (UShort_t sec = 0; sec < fESDObj->MaxSectors(); sec++) {
+         for (UShort_t str = 0; str < fESDObj->MaxStrips(); str++) {
+           if (fESDObj->Multiplicity(det, ring, sec, str) != 
+               fromEsd->Multiplicity(det, ring, sec, str))
+             AliWarning(Form("Mult for FMD%d%c[%2d,%3d]",det,ring,sec,str));
+           if (fESDObj->Eta(det, ring, sec, str) != 
+               fromEsd->Eta(det, ring, sec, str))
+             AliWarning(Form("Eta for FMD%d%c[%2d,%3d]", det,ring,sec,str));
+           if (fESDObj->Multiplicity(det, ring, sec, str) > 0 && 
+               fESDObj->Multiplicity(det, ring, sec, str) 
+               != AliESDFMD::kInvalidMult) 
+             AliInfo(Form("Mult in FMD%d%c[%2d,%3d] is %f", det,ring,sec,str,
+                          fESDObj->Multiplicity(det, ring, sec, str)));
+         }
+       }
+      }
+    }
+#endif
+  }
+
+#if 0  
+  static Int_t evNo = -1;
+  evNo++;
+  if (esd) evNo = esd->GetEventNumber();
+  TString fname(Form("FMD.ESD.%03d.root", evNo));
+  TFile* file = TFile::Open(fname.Data(), "RECREATE");
+  if (!file) {
+    AliError(Form("Failed to open file %s", fname.Data()));
+    return;
+  }
+  fESDObj->Write();
+  file->Close();
+#endif
+    
 #if 0
   TClonesArray* multStrips  = 0;
   TClonesArray* multRegions = 0;
@@ -381,6 +550,7 @@ AliFMDReconstructor::Reconstruct(AliRawReader*,TTree*) const
 {
   // Cannot be used.  See member function with same name but with 2
   // TTree arguments.   Make sure you do local reconstrucion 
+  AliDebug(1, Form("Calling FillESD with loader and tree"));
   AliError("MayNotUse");
 }
 //____________________________________________________________________
@@ -389,6 +559,7 @@ AliFMDReconstructor::Reconstruct(AliRunLoader*) const
 {
   // Cannot be used.  See member function with same name but with 2
   // TTree arguments.   Make sure you do local reconstrucion 
+  AliDebug(1, Form("Calling FillESD with loader"));
   AliError("MayNotUse");
 }
 //____________________________________________________________________
@@ -397,6 +568,7 @@ AliFMDReconstructor::Reconstruct(AliRunLoader*, AliRawReader*) const
 {
   // Cannot be used.  See member function with same name but with 2
   // TTree arguments.   Make sure you do local reconstrucion 
+  AliDebug(1, Form("Calling FillESD with loader and raw reader"));
   AliError("MayNotUse");
 }
 //____________________________________________________________________
@@ -405,6 +577,7 @@ AliFMDReconstructor::FillESD(AliRawReader*,TTree*,AliESD*) const
 {
   // Cannot be used.  See member function with same name but with 2
   // TTree arguments.   Make sure you do local reconstrucion 
+  AliDebug(1, Form("Calling FillESD with raw reader, tree, and ESD"));
   AliError("MayNotUse");
 }
 //____________________________________________________________________
@@ -413,6 +586,7 @@ AliFMDReconstructor::FillESD(AliRunLoader*,AliESD*) const
 {
   // Cannot be used.  See member function with same name but with 2
   // TTree arguments.   Make sure you do local reconstrucion 
+  AliDebug(1, Form("Calling FillESD with loader and ESD"));
   AliError("MayNotUse");
 }
 //____________________________________________________________________
@@ -421,6 +595,7 @@ AliFMDReconstructor::FillESD(AliRunLoader*,AliRawReader*,AliESD*) const
 {
   // Cannot be used.  See member function with same name but with 2
   // TTree arguments.   Make sure you do local reconstrucion 
+  AliDebug(1, Form("Calling FillESD with loader, raw reader, and ESD"));
   AliError("MayNotUse");
 }
 
index 85793bd..51b9f0b 100644 (file)
@@ -32,9 +32,6 @@
 #ifndef ALIRECONSTRUCTOR_H
 # include <AliReconstructor.h>
 #endif
-#ifndef ROOT_TObjArray
-# include <TObjArray.h>
-#endif
 
 //____________________________________________________________________
 class TTree;
@@ -43,6 +40,7 @@ class AliFMDDigit;
 class AliRawReader;
 class AliRunLoader;
 class AliESD;
+class AliESDFMD;
 
 //____________________________________________________________________
 class AliFMDReconstructor: public AliReconstructor 
@@ -73,12 +71,17 @@ private:
 protected:
   virtual void     ProcessDigits(TClonesArray* digits) const;
   virtual UShort_t SubtractPedestal(AliFMDDigit* digit) const;
+  virtual Float_t  Adc2Energy(AliFMDDigit* digit, Float_t eta, 
+                             UShort_t count) const;
+  virtual Float_t  Energy2Multiplicity(AliFMDDigit* digit, Float_t edep) const;
+  virtual void     PhysicalCoordinates(AliFMDDigit* digit, Float_t& eta, 
+                                      Float_t& phi) const;
   
-  TObjArray            fAlgorithms;    // Array of algorithms
-  Float_t               fPedestal;      // Pedestal to subtract
-  Float_t               fPedestalWidth; // Width of pedestal
-  Float_t               fPedestalFactor;// Number of pedestal widths 
+  mutable TClonesArray* fMult;          // Cache of RecPoints
+  mutable Int_t         fNMult;         // Number of entries in fMult 
+  mutable TTree*        fTreeR;         // Output tree 
   mutable Float_t       fCurrentVertex; // Z-coordinate of primary vertex
+  mutable AliESDFMD*    fESDObj;        // ESD output object
   AliESD*               fESD;
   
   ClassDef(AliFMDReconstructor, 0)  // class for the FMD reconstruction
index 2f8bbfa..e54ef17 100644 (file)
@@ -1,5 +1,5 @@
-#ifndef ALIFMDEDEPMAP_H
-#define ALIFMDEDEPMAP_H
+#ifndef ALIFMDUSHORTMAP_H
+#define ALIFMDUSHORTMAP_H
 /* Copyright(c) 1998-2000, ALICE Experiment at CERN, All rights
  * reserved. 
  *
index 3155ec7..5c29d9e 100644 (file)
@@ -137,7 +137,7 @@ AliFMDv1::VMC2FMD(Int_t copy, TLorentzVector& v,
   }
 #endif
   if (sector < 1 || sector > n) {
-    AliWarning(Form("Step", "sector # %d out of range (0-%d)", sector-1, n-1));
+    AliWarning(Form("sector # %d out of range (0-%d)", sector-1, n-1));
     return kFALSE;
   }
   sector--;
@@ -243,7 +243,7 @@ AliFMDv1::StepManager()
     if (mc->IsTrackOut())         what.Append("out ");
     
     Int_t mother = gAlice->GetMCApp()->GetPrimary(trackno);
-    AliWarning(Form("Step", "Track # %5d deposits a lot of energy\n" 
+    AliWarning(Form("Track # %5d deposits a lot of energy\n" 
                    "  Volume:    %s\n" 
                    "  Momentum:  (%7.4f,%7.4f,%7.4f)\n"
                    "  PDG:       %d (%s)\n" 
index c863044..b20dee2 100644 (file)
@@ -113,9 +113,10 @@ enum EG_t {
   kMuonCocktailCent1Single,    //
   kMuonCocktailPer1Single,     //
   kMuonCocktailPer4Single,
-  kFlatFMD1, 
-  kFlatFMD2, 
-  kFlatFMD3,
+  kFMD1Flat, 
+  kFMD2Flat, 
+  kFMD3Flat,
+  kFMDFlat,
   kEgMax
 };
 
@@ -177,7 +178,8 @@ const char* egName[kEgMax] = {
   "kMuonCocktailPer4Single",
   "kFMD1Flat",
   "kFMD2Flat",
-  "kFMD3Flat"
+  "kFMD3Flat",
+  "kFMDFlat"
 };
 
 //____________________________________________________________________
@@ -1574,10 +1576,10 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
       gGener=gener;
     }
     break;
-  case kFlatFMD1: 
+  case kFMD1Flat: 
     {
       comment = comment.Append(" Flat in FMD1 range");
-      AliGenBox* gener = AliGenBox(2000);
+      AliGenBox* gener = new AliGenBox(2000);
       gener->SetPart(211);
       gener->SetMomentumRange(3,4);
       gener->SetPhiRange(0, 360);
@@ -1585,10 +1587,10 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
       gGener = gener;
     }
     break;
-  case kFlatFMD2: 
+  case kFMD2Flat: 
     {
       comment = comment.Append(" Flat in FMD2 range");
-      AliGenBox* gener = AliGenBox(2000);
+      AliGenBox* gener = new AliGenBox(2000);
       gener->SetPart(211);
       gener->SetMomentumRange(3,4);
       gener->SetPhiRange(0, 360);
@@ -1596,10 +1598,10 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
       gGener = gener;
     }
     break;
-  case kFlatFMD3: 
+  case kFMD3Flat: 
     {
       comment = comment.Append(" Flat in FMD3 range");
-      AliGenBox* gener = AliGenBox(2000);
+      AliGenBox* gener = new AliGenBox(2000);
       gener->SetPart(211);
       gener->SetMomentumRange(3,4);
       gener->SetPhiRange(0, 360);
@@ -1607,6 +1609,25 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
       gGener = gener;
     }
     break;
+  case kFMDFlat:
+    {
+      comment = comment.Append(" Flat in FMD range");
+      AliGenCocktail* gener = AliGenCocktail("FMD cocktail");
+      gener->SetPart(211);
+      gener->SetMomentumRange(3,4);
+      gener->SetPhiRange(0, 360);
+      AliGenBox* gener3 = new AliGenBox(2000);
+      gener3->SetThetaRange(155.97, 176.73);
+      gener->AddGenerator(gener3, "FMD3", .33);
+      AliGenBox* gener2 = new AliGenBox(2000);
+      gener2->SetThetaRange(2.95, 20.42);
+      gener->AddGenerator(gener2, "FMD2", .33);
+      AliGenBox* gener1 = new AliGenBox(2000);
+      gener1->SetThetaRange(0.77, 3.08);
+      gener->AddGenerator(gener1, "FMD1", .34);
+      gGener = gener;
+    }
+    break;
     
   default: break;
   }
index 4450905..a2933f4 100644 (file)
@@ -12,8 +12,6 @@
 #pragma link C++ class  AliFMDBaseDigit+;
 #pragma link C++ class  AliFMDDigit+;
 #pragma link C++ class  AliFMDSDigit+;
-#pragma link C++ class  AliFMDMap+;
-#pragma link C++ class  AliFMDFloatMap+;
 #pragma link C++ class  AliFMDBoolMap+;
 #pragma link C++ class  AliFMDUShortMap+;
 #pragma link C++ class  AliFMD1+;
@@ -25,6 +23,7 @@
 #pragma link C++ class  AliFMDParameters+;
 #pragma link C++ class  AliFMDCalibPedestal+;
 #pragma link C++ class  AliFMDCalibGain+;
+#pragma link C++ class  AliFMDCalibSampleRate+;
 
 #else
 # error Not for compilation 
index b913499..7b9a253 100644 (file)
@@ -14,9 +14,9 @@
 // #pragma link C++ class  AliFMDMap<UShort_t>;
 // #pragma link C++ typedef AliFMDAdcMap;
 #pragma link C++ class  AliFMDReconstructor+;
-#pragma link C++ class  AliFMDMultAlgorithm+;
-#pragma link C++ class  AliFMDMultNaiive+;
-#pragma link C++ class  AliFMDMultPoisson+;
+// #pragma link C++ class  AliFMDMultAlgorithm+;
+// #pragma link C++ class  AliFMDMultNaiive+;
+// #pragma link C++ class  AliFMDMultPoisson+;
 #pragma link C++ class  AliFMDMult+;
 #pragma link C++ class  AliFMDMultRegion+;
 #pragma link C++ class  AliFMDMultStrip+;
index 1922a5f..9d64f20 100644 (file)
@@ -30,6 +30,7 @@
 #pragma link C++ class  AliFMDBaseDigitizer+;
 #pragma link C++ class  AliFMDDigitizer+;
 #pragma link C++ class  AliFMDSDigitizer+;
+#pragma link C++ class  AliFMDDigitizerAlla+;
 #pragma link C++ class  AliFMDRawWriter+;
 #pragma link C++ class  AliFMDAlla+;
 
index 20da3c4..4f6cfbb 100644 (file)
 void 
 Reconstruct()
 {
+  AliLog::SetModuleDebugLevel("FMD", 2);
   AliReconstruction rec;   
   rec.SetRunLocalReconstruction("FMD");
   rec.SetRunVertexFinder(kFALSE);
-  // rec.SetRunTracking(kFALSE); 
+  rec.SetRunTracking(""); 
   rec.SetFillESD("FMD"); 
-  rec.SetInput("./");
+  // rec.SetInput("./");
   rec.Run(); 
 }
 
index 0b5a0f1..bb66d05 100644 (file)
@@ -23,11 +23,11 @@ Simulate()
  sim.SetConfigFile("$(ALICE_ROOT)/FMD/Config.C");
  // sim.SetMakeSDigits("FMD");
  sim.SetMakeDigits("FMD");
- // sim.SetWriteRawData("FMD");
+ sim.SetWriteRawData("FMD");
  // sim.SetMakeDigitsFromHits("FMD");
  TStopwatch w;
  w.Start();
- sim.Run(10); 
+ sim.Run(1); 
  w.Stop();
  w.Print();
 }
diff --git a/FMD/TODO b/FMD/TODO
new file mode 100644 (file)
index 0000000..c9727e2
--- /dev/null
+++ b/FMD/TODO
@@ -0,0 +1,9 @@
+       Things to do in the FMD offline code
+       ------------------------------------
+
+* Split rings/cone into half rings/cones for alignment.
+* In AliFMDGeometry, get the global transformation matrices for each
+  sensitive volume (use TGeoIterator). 
+* Try implement alignment objects. 
+* Get calibration parameters (in AliFMDParameters) from CDB.
+* Implement Detector2XYZ and XYZ2Detector via TGeoMatrix. 
index 099580e..3e0b9ad 100644 (file)
@@ -3,19 +3,19 @@
 # $Id$
 
 SRCS           =  AliFMDDigit.cxx              \
-                  AliFMDMap.cxx                \
-                  AliFMDFloatMap.cxx           \
                   AliFMDBoolMap.cxx            \
                   AliFMDUShortMap.cxx          \
                   AliFMDCalibPedestal.cxx      \
                   AliFMDCalibGain.cxx          \
+                  AliFMDCalibSampleRate.cxx    \
                   AliFMDParameters.cxx         \
                   AliFMDGeometry.cxx           \
                   AliFMDRing.cxx               \
                   AliFMDDetector.cxx           \
                   AliFMD1.cxx                  \
                   AliFMD2.cxx                  \
-                  AliFMD3.cxx
+                  AliFMD3.cxx                  
+
 HDRS           =  $(SRCS:.cxx=.h)
 DHDR           := FMDbaseLinkDef.h
 # Un comment to use old code with seperate simulator task
index 1d55f19..22f4d8f 100644 (file)
@@ -5,9 +5,6 @@
 SRCS           =  AliFMDReconstructor.cxx      \
                   AliFMDRawStream.cxx          \
                   AliFMDRawReader.cxx          \
-                  AliFMDMultAlgorithm.cxx      \
-                  AliFMDMultNaiive.cxx         \
-                  AliFMDMultPoisson.cxx        \
                   AliFMDMultRegion.cxx         \
                   AliFMDMult.cxx               \
                   AliFMDMultStrip.cxx          
@@ -17,6 +14,11 @@ EINCLUDE     := $(ALICE)/RAW
 # Un comment to use old code with seperate simulator task
 # PACKCXXFLAGS := $(CXXFLAGS) -DUSE_PRE_MOVE
 
+# Not used any more 
+#                 AliFMDMultAlgorithm.cxx      \
+#                 AliFMDMultNaiive.cxx         \
+#                 AliFMDMultPoisson.cxx        \
+
 #
 # EOF
 #
index bc8cc08..739dfb8 100644 (file)
@@ -14,6 +14,7 @@ SRCS          =  AliFMD.cxx                           \
                   AliFMDG3OldSimulator.cxx             \
                   AliFMDHit.cxx                        \
                   AliFMDDigitizer.cxx                  \
+                  AliFMDDigitizerAlla.cxx              \
                   AliFMDEdepMap.cxx                    \
                   AliFMDRawWriter.cxx                  \
                   AliFMDAlla.cxx               
diff --git a/FMD/scripts/CompareESD.C b/FMD/scripts/CompareESD.C
new file mode 100644 (file)
index 0000000..c1e6d45
--- /dev/null
@@ -0,0 +1,92 @@
+#ifndef __CINT__
+#include <AliESD.h>
+#include <TObjArray.h>
+#include <TFile.h>
+#include <TSystemDirectory.h>
+#include <TString.h>
+#include <TChain.h>
+#include <iostream>
+#endif
+
+void
+CompareESD()
+{
+  TChain*           chain = new TChain("esdTree");
+  TSystemDirectory* dir   = new TSystemDirectory(".", ".");
+  TList*            files = dir->GetListOfFiles();
+  files->Sort();
+  TIter next(files);
+  TSystemFile*      file  = 0;
+  TObjArray*        other = new TObjArray;
+  while ((file = static_cast<TSystemFile*>(next()))) {
+    TString fname(file->GetName());
+    if (fname.Contains("AliESDs")) 
+      chain->AddFile(fname.Data());
+    if (fname.Contains("FMD.ESD.")) {
+      TFile* o = TFile::Open(fname.Data());
+      if (!o) {
+       Warning("CompareESD", "Failed to open file %s", fname.Data());
+       other->Add(0);
+       continue;
+      }
+      AliESDFMD* fmdEsd = static_cast<AliESDFMD*>(o->Get("AliESDFMD"));
+      if (!fmdEsd) {
+       Warning("CompareESD", "Failed to get FMD object from reference %s", 
+               fname.Data());
+       other->Add(0);
+       continue;
+      }
+      other->Add(fmdEsd);
+    }
+  }
+  delete dir;
+  
+  AliESD* esd = 0;
+  chain->SetBranchAddress("ESD", &esd);
+  Int_t   n   = chain->GetEntries();
+  for (Int_t i = 0; i < n; i++) {
+    Int_t read = chain->GetEntry(i);
+    if (read <= 0) break;
+    std::cout << "Entry # " << i << std::endl;
+    if (!esd) break;
+    esd->Print();
+    AliESDFMD* esdObj = esd->GetFMDData();
+    AliESDFMD* refObj = static_cast<AliESDFMD*>(other->At(i));
+    if (!esdObj) {
+      Warning("CompareESD", "no FMD object in ESD");
+      continue;
+    }
+    if (!refObj) {
+      Warning("CompareESD", "no reference FMD object");
+      continue;
+    }
+    std::cout << " Comparing ... " << std::endl;
+    for (UShort_t det = 1; det <= esdObj->MaxDetectors(); det++) {
+      std::cout << "FMD" << det << std::endl;
+      for (UShort_t ir = 0; ir < esdObj->MaxRings(); ir++) {
+       Char_t ring = (ir == 0 ? 'I' : 'O');
+       std::cout << " Ring " << ring << std::endl;
+       for (UShort_t sec = 0; sec < esdObj->MaxSectors(); sec++) {
+         std::cout << "  Sector # " << sec << std::endl;
+         for (UShort_t str = 0; str < esdObj->MaxStrips(); str++) {
+           if (esdObj->Multiplicity(det, ring, sec, str) != 
+               refObj->Multiplicity(det, ring, sec, str))
+             Warning("CompareESD", 
+                     "Mult for FMD%d%c[%2d,%3d]",det,ring,sec,str);
+           if (esdObj->Eta(det, ring, sec, str) != 
+               refObj->Eta(det, ring, sec, str))
+             Warning("CompareESD", 
+                     "Eta for FMD%d%c[%2d,%3d]", det,ring,sec,str);
+           if (esdObj->Multiplicity(det, ring, sec, str) > 0) 
+             Info("CompareESD", "Mult in FMD%c%d[%2d,%3d] is %f",
+                  det,ring,sec,str,
+                  esdObj->Multiplicity(det, ring, sec, str));
+         }
+       }
+      }
+    }
+  }
+}
+
+      
+  
index 8e83cc7..db8ec10 100644 (file)
@@ -5,14 +5,14 @@
 // include path to contain the relevant directories. 
 //
 void
-Compile(const char* script)
+Compile(const char* script, Option_t* option="g")
 {
   gSystem->Load("libFMDutil.so");
   gSystem->SetIncludePath("-I`root-config --incdir` "
                          "-I${ALICE_ROOT}/include " 
                          "-I${ALICE_ROOT}/FMD "
                          "-I${ALICE_ROOT}/geant3/TGeant3");
-  gROOT->ProcessLine(Form(".L %s++g", script));
+  gROOT->ProcessLine(Form(".L %s+%s", script, option));
 }
 
 //____________________________________________________________________
diff --git a/FMD/scripts/DrawDigitsRecs.C b/FMD/scripts/DrawDigitsRecs.C
new file mode 100644 (file)
index 0000000..23c9336
--- /dev/null
@@ -0,0 +1,183 @@
+//
+// $Id$
+//
+// 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 <TH2D.h>
+#include <AliFMDHit.h>
+#include <AliFMDDigit.h>
+#include <AliFMDInput.h>
+#include <AliFMDUShortMap.h>
+#include <AliFMDFloatMap.h>
+#include <AliFMDMultStrip.h>
+#include <AliFMDMultRegion.h>
+#include <iostream>
+#include <TStyle.h>
+#include <TArrayF.h>
+#include <TCanvas.h>
+
+class DrawDigitsRecs : public AliFMDInputDigits
+{
+private:
+  TH2D* fAdcVsSingle; // Histogram 
+  TH2D* fAdcVsRegion; // Histogram 
+  TH2D* fSingleVsRegion; // Histogram 
+  AliFMDUShortMap fMap;
+  AliFMDFloatMap fEta;
+  AliFMDFloatMap fPhi;
+  AliFMDFloatMap fMult;
+public:
+  TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max) 
+  {
+    TArrayF bins(n+1);
+    Float_t dp   = n / TMath::Log10(max / min);
+    Float_t pmin = TMath::Log10(min);
+    bins[0]      = min;
+    for (Int_t i = 1; i < n+1; i++) {
+      Float_t p = pmin + i / dp;
+      bins[i]   = TMath::Power(10, p);
+    }
+    return bins;
+  }
+  DrawDigitsRecs(Int_t m=1100, Double_t amin=-0.5, Double_t amax=1099.5,
+                Int_t n=105, Double_t mmin=-0.5, Double_t mmax=20.5) 
+  { 
+    AddLoad(kRecPoints);
+    fAdcVsSingle = new TH2D("adcVsSingle", "ADC vs. Multiplicity (strip)", 
+                           m, amin, amax, n, mmin, mmax);
+    fAdcVsSingle->SetXTitle("ADC value");
+    fAdcVsSingle->SetYTitle("Strip Multiplicity");
+
+    mmin *= 20;
+    mmax *= 20;
+    amin *= 50;
+    amax *= 50;
+    fAdcVsRegion = new TH2D("adcVsRegion", "ADC vs Muliplcity (region)", 
+                           m, amin, amax, n, mmin, mmax);
+    fAdcVsRegion->SetXTitle("ADC value");
+    fAdcVsRegion->SetYTitle("Region Multiplicity");
+
+    fSingleVsRegion = new TH2D("singleVsRegion", "Single vs Region", 
+                              n, mmin, mmax, n, mmin, mmax);
+    fSingleVsRegion->SetXTitle("Strip Multiplicity");
+    fSingleVsRegion->SetYTitle("Region Multiplicity");
+  }
+  Bool_t Begin(Int_t ev) 
+  {
+    fMap.Reset();
+    fEta.Reset();
+    fPhi.Reset();
+    return AliFMDInputDigits::Begin(ev);
+  }
+  Bool_t ProcessDigit(AliFMDDigit* digit) 
+  {
+    // Cache the energy loss 
+    if (!digit) return kFALSE;
+    UShort_t det = digit->Detector();
+    Char_t   rng = digit->Ring();
+    UShort_t sec = digit->Sector();
+    UShort_t str = digit->Strip();
+    if (str > 511) {
+      AliWarning(Form("Bad strip number %d in digit", str));
+      return kTRUE;
+    }
+    fMap(det, rng, sec, str) = digit->Counts();
+    return kTRUE;
+  }
+  Bool_t Event() 
+  {
+    if (!AliFMDInputDigits::Event()) return kFALSE;
+    Int_t nEv = fTreeR->GetEntries();
+    for (Int_t i = 0; i < nEv; i++) {
+      Int_t recoRead  = fTreeR->GetEntry(i);
+      if (recoRead <= 0) continue;
+      Int_t nSingle = fArrayN->GetEntries();
+      if (nSingle <= 0) continue;
+      for (Int_t j = 0; j < nSingle; j++) {
+       AliFMDMultStrip* single = 
+         static_cast<AliFMDMultStrip*>(fArrayN->At(j));
+       if (!single) continue;
+       UShort_t det = single->Detector();
+       Char_t   rng = single->Ring();
+       UShort_t sec = single->Sector();
+       UShort_t str = single->Strip();
+       if (str > 511) {
+         AliWarning(Form("Bad strip number %d in single", str));
+         continue;
+       }
+       fAdcVsSingle->Fill(fMap(det, rng, sec, str), single->Particles());
+       fEta(det, rng, sec, str)  = single->Eta();
+       fPhi(det, rng, sec, str)  = single->Phi() * 180 / TMath::Pi();
+       fMult(det, rng, sec, str) = single->Particles();
+      }    
+
+      Int_t nRegion = fArrayP->GetEntries();
+      if (nRegion <= 0) continue;
+      for (Int_t j = 0; j < nRegion; j++) {
+       AliFMDMultRegion* region = 
+         static_cast<AliFMDMultRegion*>(fArrayP->At(j));
+       if (!region) continue;
+       UShort_t det  = region->Detector();
+       Char_t   rng  = region->Ring();
+       UInt_t   suma = 0;
+       Float_t  summ = 0;
+       
+       for (size_t sec = 0; sec < AliFMDMap::kMaxSectors; sec++) {
+         Float_t phi = fPhi(det, rng, sec, 0);     
+         Bool_t ok   = (phi <= region->MaxPhi() && phi >= region->MinPhi());
+         if (!ok) continue;
+         for (size_t str = 0; str < AliFMDMap::kMaxStrips; str++) {
+           Float_t eta  = fEta(det, rng, sec, str);
+           Float_t sign = eta < 0 ? -1. : 1.;
+           ok           = (sign * eta <= sign * region->MaxEta() 
+                           && sign * eta >= sign * region->MinEta());
+           if (!ok) continue;
+           suma += fMap(det, rng, sec, str);
+           summ += fMult(det, rng, sec, str);
+         }
+       }
+       fAdcVsRegion->Fill(suma, region->Particles()); 
+       fSingleVsRegion->Fill(summ, region->Particles());
+     }    
+    }
+    return kTRUE;
+  }
+  Bool_t Finish()
+  {
+    gStyle->SetPalette(1);
+    gStyle->SetOptTitle(0);
+    gStyle->SetCanvasColor(0);
+    gStyle->SetCanvasBorderSize(0);
+    gStyle->SetPadColor(0);
+    gStyle->SetPadBorderSize(0);
+
+    new TCanvas("c1", "ADC vs. Strip Multiplicity");
+    fAdcVsSingle->SetStats(kFALSE);
+    fAdcVsSingle->Draw("COLZ");
+
+    new TCanvas("c2", "ADC vs. Region Multiplicity");
+    fAdcVsRegion->SetStats(kFALSE);
+    fAdcVsRegion->Draw("COLZ");
+
+    new TCanvas("c3", "Single vs. Region Multiplicity");
+    fSingleVsRegion->SetStats(kFALSE);
+    fSingleVsRegion->Draw("COLZ");
+
+    return kTRUE;
+  }
+
+  ClassDef(DrawDigitsRecs,0);
+  
+};
+
+//____________________________________________________________________
+//
+// EOF
+//
diff --git a/FMD/scripts/DrawHitPatterns.C b/FMD/scripts/DrawHitPatterns.C
new file mode 100644 (file)
index 0000000..1c3bca2
--- /dev/null
@@ -0,0 +1,221 @@
+//
+// $Id$
+//
+// Script that contains a class to draw hits, 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 <TH2D.h>
+#include <AliFMDHit.h>
+#include <AliFMDInput.h>
+#include <iostream>
+#include <TStyle.h>
+#include <TArrayF.h>
+#include <TMarker3DBox.h>
+#include <TGeoManager.h>
+#include <TMath.h>
+#include <TApplication.h>
+#include <TButton.h>
+#include <TParticle.h>
+#include <TCanvas.h>
+#include <TView.h>
+#include <TVirtualX.h>
+
+class DrawHitPatterns : public AliFMDInputHits
+{
+private:
+  Bool_t     fWait;
+  Bool_t     fPrimary;
+  Int_t      fPdg;
+  TObjArray* fMarkers;
+  TObjArray* fHits;
+  TCanvas*   fCanvas;
+  TPad*      fPad;
+  TButton*   fButton;
+  TButton*   fZoom;
+  TButton*   fPick;
+  Bool_t     fZoomMode;
+public:
+  static DrawHitPatterns* fgInstance;
+  DrawHitPatterns(Bool_t primary=kTRUE, Int_t pdg=211) 
+    : fPrimary(primary), fPdg(pdg), 
+      fCanvas(0), fPad(0), fButton(0)
+  { 
+    AddLoad(kKinematics);
+    AddLoad(kGeometry);
+    fMarkers = new TObjArray;
+    fHits    = new TObjArray;
+    fMarkers->SetOwner(kTRUE);
+    fHits->SetOwner(kFALSE);
+    fgInstance = this;
+  }
+  void Continue() { fWait = kFALSE; }
+  void Zoom() { fZoomMode = kTRUE; }
+  void Pick() { fZoomMode = kFALSE; }
+  void ExecuteEvent(Int_t event, Int_t px, Int_t py) 
+  {
+    Info("ExecuteEvent", "Event %d, at (%d,%d)", px, py);
+    static Float_t x0, y0, x1, y1;
+    static Int_t   px0, py0, pxOld,pyOld;
+    static Bool_t  lineDraw;
+    if (px == 0 && py == 0) return;
+    if (!fZoomMode && fPad->GetView()) {
+      fPad->GetView()->ExecuteRotateView(event, px, py);
+      return;
+    }
+    gPad->SetCursor(kCross);
+    switch (event) {
+    case kButton1Down: 
+      gPad->TAttLine::Modify();
+      x0 = fPad->AbsPixeltoX(px);
+      y0 = fPad->AbsPixeltoY(py);
+      px0 = pxOld = px;
+      py0 = pyOld = py;
+      lineDraw = kFALSE;
+      return;
+    case kButton1Motion:
+      if (lineDraw) 
+       gVirtualX->DrawBox(px0, py0, pxOld, pyOld, TVirtualX::kHollow);
+      pxOld = px;
+      pyOld = py;
+      lineDraw = kTRUE;
+      gVirtualX->DrawBox(px0, py0, pxOld, pyOld, TVirtualX::kHollow);
+      return;
+    case kButton1Up:
+      fPad->GetCanvas()->FeedbackMode(kFALSE);
+      if (px == px0 || py == py0) return;
+      x1 = gPad->AbsPixeltoX(px);
+      y1 = gPad->AbsPixeltoY(py);
+      if (x1 < x0) std::swap(x0, x1); 
+      if (y1 < y0) std::swap(y0, y1); 
+      gPad->Range(x0, y0, x1, y1);
+      gPad->Modified();
+      return;
+    }
+  }
+  Int_t  DistanceToPrimitive(Int_t px, Int_t py) 
+  {
+    Info("DistanceToPrimitive", "@ (%d,%d)", px, py);
+    gPad->SetCursor(kCross);
+    Float_t xmin = gPad->GetX1();
+    Float_t xmax = gPad->GetX2();
+    Float_t dx   = .02 * (xmax - xmin);
+    Float_t x    = gPad->AbsPixeltoX(px);
+    if (x < xmin + dx || x > xmax - dx) return 9999;
+    return (fZoomMode ? 0 : 7);
+  }
+  Bool_t Begin(Int_t event) 
+  {
+    if (!fCanvas) {
+      fCanvas = new TCanvas("display", "Display", 700, 700);
+      fCanvas->SetFillColor(1);
+      fCanvas->ToggleEventStatus();
+      fPad = new TPad("view3D", "3DView", 0.0, 0.1, 1.0, 1.0, 1, 0, 0);
+      fCanvas->cd();
+      fPad->Draw();
+    }
+    if (!fButton) {
+      fCanvas->cd();
+      fButton = new TButton("Continue", 
+                           "DrawHitPatterns::fgInstance->Continue()",
+                           0, 0, .5, .05);
+      fButton->Draw();
+      fZoom = new TButton("Zoom", 
+                         "DrawHitPatterns::fgInstance->Zoom()",
+                         .5, 0, .75, .05);
+      fZoom->Draw();
+      fPick = new TButton("Pick", 
+                         "DrawHitPatterns::fgInstance->Pick()",
+                         .75, 0, 1, .05);
+      fPick->Draw();
+    }
+    Info("Begin", "Clearing canvas");
+    // fCanvas->Clear();
+    if (!fGeoManager) {
+      Warning("End", "No geometry manager");
+      return kFALSE;
+    }
+    Info("Begin", "Drawing geometry");
+    fPad->cd();
+    fGeoManager->GetTopVolume()->Draw();
+    Info("Begin", "Adjusting view");
+    Int_t irep;
+    if (fPad->GetView()) {
+      fPad->GetView()->SetView(-200, -40, 80, irep);
+      fPad->GetView()->Zoom();
+      fPad->Modified();
+      fPad->cd();
+    }
+    Info("Begin", "Drawing button");
+
+    return AliFMDInputHits::Begin(event);
+  }
+    
+  Bool_t ProcessHit(AliFMDHit* hit, TParticle* p) 
+  {
+    if (!hit) {
+      std::cout << "No hit" << std::endl;
+      return kFALSE;
+    }
+
+    if (!p) {
+      std::cout << "No track" << std::endl;
+      return kFALSE;
+    }
+    if (fPrimary && !p->IsPrimary()) return kTRUE;
+    if (hit->IsStop()) return kTRUE;
+    if (fPdg != 0) {
+      if (TMath::Abs(hit->Pdg()) != fPdg) return kTRUE;
+    }
+    fHits->Add(hit);
+    Float_t  size  = TMath::Min(TMath::Max(hit->Edep() * .1, .1), 1.);
+    Float_t  pt    = TMath::Sqrt(hit->Py()*hit->Py()+hit->Px()*hit->Px());
+    Float_t  theta = TMath::ATan2(pt, hit->Pz());
+    Float_t  phi   = TMath::ATan2(hit->Py(), hit->Px());
+    TMarker3DBox* marker = new  TMarker3DBox(hit->X(), hit->Y(), hit->Z(),
+                                            size, size, size,
+                                            theta, phi);
+    marker->SetLineColor((p->IsPrimary() ? 3 : 4));
+    marker->SetRefObject(hit);
+    fMarkers->Add(marker);
+    return kTRUE;
+  }
+  Bool_t End() 
+  {
+    // fGeoManager->GetTopVolume()->Draw();
+    fPad->cd();
+    fMarkers->Draw();
+    AppendPad();
+    fPad->Modified(kTRUE);
+    fPad->Update();
+    fPad->cd();
+    fCanvas->Modified(kTRUE);
+    fCanvas->Update();
+    fCanvas->cd();
+    
+    fWait = kTRUE;
+    while (fWait) {
+      gApplication->StartIdleing();
+      gSystem->InnerLoop();
+      gApplication->StopIdleing();
+    }
+    Info("End", "After idle loop");
+    fMarkers->Delete();
+    fHits->Clear();
+    Info("End", "After clearing caches");
+    return AliFMDInputHits::End();
+  }
+  ClassDef(DrawHitPatterns,0);
+};
+
+DrawHitPatterns* DrawHitPatterns::fgInstance = 0;
+
+//____________________________________________________________________
+//
+// EOF
+//
index eb82a1c..fdbb614 100644 (file)
@@ -16,6 +16,7 @@
 #include <iostream>
 #include <TStyle.h>
 #include <TArrayF.h>
+#include <TParticle.h>
 
 class DrawHits : public AliFMDInputHits
 {
@@ -26,9 +27,13 @@ public:
   TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max) 
   {
     TArrayF bins(n+1);
+    bins[0]      = min;
+    if (n <= 20) {
+      for (Int_t i = 1; i < n+1; i++) bins[i] = bins[i-1] + (max-min)/n;
+      return bins;
+    }
     Float_t dp   = n / TMath::Log10(max / min);
     Float_t pmin = TMath::Log10(min);
-    bins[0]      = min;
     for (Int_t i = 1; i < n+1; i++) {
       Float_t p = pmin + i / dp;
       bins[i]   = TMath::Power(10, p);
@@ -40,20 +45,30 @@ public:
           Int_t pdg=211) 
     : fPdg(pdg)
   { 
+    AddLoad(kKinematics);
     TArrayF tkine(MakeLogScale(n, tmin, tmax));
     TArrayF eloss(MakeLogScale(m, emin, emax));
-    fElossVsPMQ = new TH2D("bad", "#Delta E/#Delta x vs. p/(mq^{2})", 
-                          tkine.fN-1, tkine.fArray, eloss.fN-1, eloss.fArray);
+    TString name(Form("elossVsP%s", (fPdg == 0 ? "MQ" : "")));
+    TString title(Form("#Delta E/#Delta x vs. p%s", 
+                      (fPdg == 0 ? "/(mq^{2})" : "")));
+    fElossVsPMQ = new TH2D(name.Data(), title.Data(), 
+                          tkine.fN-1, tkine.fArray, 
+                          eloss.fN-1, eloss.fArray);
     fElossVsPMQ->SetXTitle(Form("p%s", (fPdg == 0 ? "/(mq^{2})" : " [GeV]")));
     fElossVsPMQ->SetYTitle("#Delta E/#Delta x [MeV/cm]");
   }
-  Bool_t ProcessHit(AliFMDHit* hit, TParticle*) 
+  Bool_t ProcessHit(AliFMDHit* hit, TParticle* p) 
   {
     if (!hit) {
       std::cout << "No hit" << std::endl;
       return kFALSE;
     }
 
+    if (!p) {
+      std::cout << "No track" << std::endl;
+      return kFALSE;
+    }
+    if (!p->IsPrimary()) return kTRUE;
     if (hit->IsStop()) return kTRUE;
     Float_t x = hit->P();
     Float_t y = hit->Edep()/hit->Length();
@@ -78,9 +93,10 @@ public:
     gStyle->SetPadColor(0);
     gStyle->SetPadBorderSize(0);
     fElossVsPMQ->SetStats(kFALSE);
-    fElossVsPMQ->Draw("COLZ");
+    fElossVsPMQ->Draw("COLZ box");
     return kTRUE;
   }
+  ClassDef(DrawHits,0);
 };
 
 //____________________________________________________________________
index 3e43c80..cf8471e 100644 (file)
@@ -63,6 +63,10 @@ public:
     Char_t   rng = hit->Ring();
     UShort_t sec = hit->Sector();
     UShort_t str = hit->Strip();
+    if (str > 511) {
+      AliWarning(Form("Bad strip number %d in hit", str));
+      return kTRUE;
+    }
     fMap(det, rng, sec, str).fEdep += hit->Edep();
     fMap(det, rng, sec, str).fN++;
     return kTRUE;
@@ -83,6 +87,10 @@ public:
        Char_t   rng = digit->Ring();
        UShort_t sec = digit->Sector();
        UShort_t str = digit->Strip();
+       if (str > 511) {
+         AliWarning(Form("Bad strip number %d in digit", str));
+         continue;
+       }
        fElossVsAdc->Fill(fMap(det, rng, sec, str).fEdep, digit->Counts());
       }    
     }
@@ -100,6 +108,8 @@ public:
     fElossVsAdc->Draw("COLZ");
     return kTRUE;
   }
+
+  ClassDef(DrawHitsDigits,0);
 };
 
 //____________________________________________________________________
diff --git a/FMD/scripts/DrawHitsRecs.C b/FMD/scripts/DrawHitsRecs.C
new file mode 100644 (file)
index 0000000..4f3ae0b
--- /dev/null
@@ -0,0 +1,167 @@
+//
+// $Id$
+//
+// 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 <TH2D.h>
+#include <AliFMDHit.h>
+#include <AliFMDMultStrip.h>
+#include <AliFMDMultRegion.h>
+#include <AliFMDDigit.h>
+#include <AliFMDInput.h>
+#include <AliFMDEdepMap.h>
+#include <AliFMDFloatMap.h>
+#include <iostream>
+#include <TStyle.h>
+#include <TArrayF.h>
+#include <TCanvas.h>
+#include <TParticle.h>
+#include <TClonesArray.h>
+#include <TTree.h>
+#include <AliStack.h>
+#include <AliLog.h>
+
+class DrawHitsRecs : public AliFMDInputHits
+{
+private:
+  TH2D* fElossVsMult; // Histogram 
+  TH2D* fHitsVsStrip;  // Histogram 
+  TH2D* fHitsVsRegion;  // Histogram 
+  AliFMDEdepMap fMap;
+  AliFMDFloatMap fEta;
+  AliFMDFloatMap fPhi;
+  AliFMDFloatMap fMult;
+  Bool_t fPrimary;
+public:
+  TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max) 
+  {
+    TArrayF bins(n+1);
+    Float_t dp   = n / TMath::Log10(max / min);
+    Float_t pmin = TMath::Log10(min);
+    bins[0]      = min;
+    for (Int_t i = 1; i < n+1; i++) {
+      Float_t p = pmin + i / dp;
+      bins[i]   = TMath::Power(10, p);
+    }
+    return bins;
+  }
+  DrawHitsRecs(Bool_t primary=kFALSE,
+              Int_t n=900, Double_t emin=1e-3, Double_t emax=10, 
+              Int_t m=21, Double_t mmin=-0.5, Double_t mmax=20.5) 
+  { 
+    fPrimary = primary;
+    AddLoad(kRecPoints);
+    if (fPrimary) AddLoad(kKinematics);
+    TArrayF eloss(MakeLogScale(n, emin, emax));
+    TArrayF mults(m+1);
+    mults[0] = mmin;
+    for (Int_t i = 1; i < m+1; i++) mults[i] = mults[i-1] + (mmax-mmin)/m;
+    fElossVsMult = new TH2D("elossVsMult", 
+                           "#Delta E vs. Multiplicity (single)", 
+                          eloss.fN-1, eloss.fArray, mults.fN-1, mults.fArray);
+    fElossVsMult->SetXTitle("#Delta E/#Delta x [MeV/cm]");
+    fElossVsMult->SetYTitle("Strip Multiplicity");
+
+    Double_t omin = -.5;
+    Double_t omax = 7.5;
+    Int_t    o    = 8;
+    fHitsVsStrip = new TH2D("hitsVsStrip", 
+                           "# of Hits vs. Multiplicity (strip)",
+                          o, omin, omax, m, mmin, mmax);
+    fHitsVsStrip->SetXTitle("# of Hits");
+    fHitsVsStrip->SetYTitle("Strip Multiplicity");
+  }
+  Bool_t Begin(Int_t ev) 
+  {
+    fMap.Reset();
+    return AliFMDInputHits::Begin(ev);
+  }
+  Bool_t ProcessHit(AliFMDHit* hit, TParticle*) 
+  {
+    // Cache the energy loss 
+    if (!hit) return kFALSE;
+    UShort_t det = hit->Detector();
+    Char_t   rng = hit->Ring();
+    UShort_t sec = hit->Sector();
+    UShort_t str = hit->Strip();
+    if (str > 511) {
+      AliWarning(Form("Bad strip number %d in hit", str));
+      return kTRUE;
+    }
+    if (fPrimary) {
+      TParticle* kine = fStack->Particle(hit->Track());
+      if (!kine) return kTRUE;
+      if (!kine->IsPrimary()) return kTRUE;
+    }
+    
+    fMap(det, rng, sec, str).fEdep += hit->Edep();
+    fMap(det, rng, sec, str).fN++;
+    return kTRUE;
+  }
+  Bool_t Event() 
+  {
+    if (!AliFMDInputHits::Event()) return kFALSE;
+    Int_t nEv = fTreeR->GetEntries();
+    for (Int_t i = 0; i < nEv; i++) {
+      Int_t singleRead  = fTreeR->GetEntry(i);
+      if (singleRead <= 0) continue;
+      Int_t nSingle = fArrayN->GetEntries();
+      if (nSingle <= 0) continue;
+      for (Int_t j = 0; j < nSingle; j++) {
+       AliFMDMultStrip* single = 
+         static_cast<AliFMDMultStrip*>(fArrayN->At(j));
+       if (!single) continue;
+       UShort_t det = single->Detector();
+       Char_t   rng = single->Ring();
+       UShort_t sec = single->Sector();
+       UShort_t str = single->Strip();
+       if (str > 511) {
+         AliWarning(Form("Bad strip number %d in single", str));
+         continue;
+       }
+       if (fMap(det, rng, sec, str).fEdep > 0) 
+         fElossVsMult->Fill(fMap(det, rng, sec, str).fEdep, 
+                            single->Particles());
+       if (fMap(det, rng, sec, str).fN > 0) 
+         fHitsVsStrip->Fill(fMap(det, rng, sec, str).fN, 
+                            single->Particles());
+       fEta(det, rng, sec, str)  = single->Eta();
+       fPhi(det, rng, sec, str)  = single->Phi() * 180 / TMath::Pi();
+      }      
+    }
+    return kTRUE;
+  }
+  Bool_t Finish()
+  {
+    gStyle->SetPalette(1);
+    gStyle->SetOptTitle(0);
+    gStyle->SetCanvasColor(0);
+    gStyle->SetCanvasBorderSize(0);
+    gStyle->SetPadColor(0);
+    gStyle->SetPadBorderSize(0);
+
+    new TCanvas("c1", "Energy loss vs. Strip Multiplicity");
+    fElossVsMult->SetStats(kFALSE);
+    fElossVsMult->Draw("COLZ");
+
+    new TCanvas("c2", "# of Hits vs. Strip Multiplicity");
+    fHitsVsStrip->SetStats(kFALSE);
+    fHitsVsStrip->Draw("COLZ");
+
+    return kTRUE;
+  }
+
+  ClassDef(DrawHitsRecs,0);
+};
+
+//____________________________________________________________________
+//
+// EOF
+//
diff --git a/FMD/scripts/TestMapIO.C b/FMD/scripts/TestMapIO.C
new file mode 100644 (file)
index 0000000..527d8af
--- /dev/null
@@ -0,0 +1,60 @@
+void 
+WriteTree()
+{
+  TFile* file = TFile::Open("map.root", "RECREATE");
+  TTree* tree = new TTree("T", "T");
+  AliFMDFloatMap* m = new AliFMDFloatMap(1, 1, 1, 3);
+  tree->Branch("map", "AliFMDFloatMap", &m);
+  for (int i = 0; i < 3; i++) m->operator()(1,'I',0,i) = i + 1;
+  tree->Fill();
+  file->Write();
+  file->Close();
+}
+
+void 
+ReadTree()
+{
+  TFile* file = TFile::Open("map.root", "READ");
+  TTree* tree = static_cast<TTree*>(file->Get("T"));
+  AliFMDFloatMap* m = 0;
+  tree->SetBranchAddress("map", &m);
+  tree->GetEntry(0);
+  for (int i = 0; i < 3; i++) {
+    std::cout << "Map(1,'I',0," << i << "): " << m->operator()(1,'I',0,i)
+             << std::endl;
+  }
+  file->Close();
+}
+
+  
+void
+WriteMap() 
+{
+  TFile* file = TFile::Open("map.root", "RECREATE");
+  AliFMDFloatMap* m = new AliFMDFloatMap(1, 1, 1, 3);
+  for (int i = 0; i < 3; i++) m->operator()(1,'I',0,i) = i + 1;
+  m.Write("map");
+  file->Close();
+}
+
+ReadMap() 
+{
+  TFile* file = TFile::Open("map.root", "READ");
+  AliFMDFloatMap* m = static_cast<AliFMDFloatMap*>(file->Get("map"));
+  std::cout << "Got map " << map << std::endl;
+  for (int i = 0; i < 3; i++) {
+    std::cout << "Map(1,'I',0," << i << "): " << m->operator()(1,'I',0,i)
+             << std::endl;
+  }
+  file->Close();
+}
+
+
+void
+TestMapIO()
+{
+  WriteMap();
+  ReadMap();
+  WriteTree();
+  ReadTree();
+}