Fixing bugs in FMD reconstruction. Everything should work now
authorhdalsgaa <hdalsgaa@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 21 Oct 2009 14:14:45 +0000 (14:14 +0000)
committerhdalsgaa <hdalsgaa@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 21 Oct 2009 14:14:45 +0000 (14:14 +0000)
FMD/AliFMDDigit.h
FMD/AliFMDESDRevertexer.cxx
FMD/AliFMDESDRevertexer.h
FMD/AliFMDRawReader.cxx
FMD/AliFMDReconstructor.cxx
FMD/FMDrecLinkDef.h
FMD/libFMDrec.pkg

index 53cec00..466f2e9 100644 (file)
@@ -105,6 +105,22 @@ public:
    * @param c Counts 
    */
   void SetCount(UShort_t s, Short_t c);
+  /**
+   * Initialize all counts to appropriate values for this oversampling
+   * rate.  That is 
+   *
+   * @verbatim
+   *     Rate | Sample 1 | Sample 2 | Sample 3 | Sample 4 
+   *     -----+----------+----------+----------+----------
+   *     1    | 0        | -1       | -1       | -1
+   *     2    | 0        | 0        | -1       | -1
+   *     3    | 0        | 0        | 0        | -1
+   *     4    | 0        | 0        | 0        | 0
+   * @endverbatim
+   *
+   * @param rate Oversampling rate 
+   */
+  void SetDefaultCounts(UShort_t rate);
 protected:
   UShort_t fCount1;     // Digital signal 
   Short_t  fCount2;     // Digital signal (-1 if not used)
@@ -113,6 +129,20 @@ protected:
   ClassDef(AliFMDDigit,2)     // Normal FMD digit
 };
 
+inline void
+AliFMDDigit::SetDefaultCounts(UShort_t rate)
+{
+  switch (rate) { 
+  case 4: fCount4 = 0; // Fall through 
+  case 3: fCount3 = 0; // Fall through 
+  case 2: fCount2 = 0; // Fall through 
+  case 1: fCount1 = 0;
+    break;
+  default: 
+    fCount4 = fCount3 = fCount2 = fCount1 = 0;
+    break;
+  }
+}
 inline UShort_t 
 AliFMDDigit::Counts() const 
 {
index a0624b6..e993cbd 100644 (file)
@@ -1,18 +1,25 @@
 #include <AliFMDESDRevertexer.h>
 #include <AliFMDGeometry.h>
 #include <AliESDFMD.h>
+#include <TMath.h>
+#include <AliLog.h>
+
+ClassImp(AliFMDESDRevertexer)
+#if 0 // for emacs 
+;
+#endif
 
 //____________________________________________________________________
 AliFMDESDRevertexer::AliFMDESDRevertexer()
 {
   AliFMDGeometry* geom = AliFMDGeometry::Instance();
   geom->Init();
-  geom->InitTransforms();
+  geom->InitTransformations();
 }
 
 //____________________________________________________________________
 Bool_t
-AliFMDESDRevertexer::Revertex(AliESDFMD* fmdEsd, Double_t vz)
+AliFMDESDRevertexer::Revertex(AliESDFMD* fmdEsd, Double_t vz) const
 {
   if (!fmdEsd) return kFALSE;
   
@@ -24,30 +31,33 @@ AliFMDESDRevertexer::Revertex(AliESDFMD* fmdEsd, Double_t vz)
       UShort_t nsec = (ir == 0 ?  20 :  40);
       UShort_t nstr = (ir == 0 ? 512 : 256);
       for (UShort_t str = 0; str < nstr; str++) { 
-       Double_t phi, r, theta, oldTheta;
+       Double_t phi, r, theta;
        Double_t eta      = AliESDFMD::kInvalidEta;
-       Double_t oldEta   = fmdEsd->Eta(det, rng, sec, str);
+       Double_t oldEta   = fmdEsd->Eta(det, rng, 0, str);
+       if (oldEta == AliESDFMD::kInvalidEta) continue;
+
        Double_t oldTheta = Eta2Theta(oldEta);
        Bool_t   ret1     = PhysicalCoordinates(det, rng, 0, str, vz, 
-                                               eta, phi, r, theta));
-       fmdEsd->SetEta(det, rng, sec, str, eta);
+                                               eta, phi, r, theta);
+       fmdEsd->SetEta(det, rng, 0, str, eta);
 
        if (!ret1) {
          // If the was an error, then there's no reason to go on with
          // this strip-ring.  Note, that the eta is correctly set to
          // AliESDFMD::kInvalidMult. 
          AliWarning(Form("Failed to calculate eta, phi for "
-                         "FMD%d%c[%02d,%03d] with v_z=%9.4f" 
+                         "FMD%d%c[%02d,%03d] with v_z=%9.4f",
                          det, rng, 0, str, vz));
          ret = kFALSE;
          continue;
        }
-       
+
        Double_t corr = TMath::Abs(TMath::Cos(theta));
        if (fmdEsd->IsAngleCorrected()) 
          corr /= TMath::Abs(TMath::Cos(oldTheta));
        for (UShort_t sec = 0; sec < nsec; sec++) { 
          Double_t mult = fmdEsd->Multiplicity(det, rng, sec, str);
+         if (mult == AliESDFMD::kInvalidMult) continue;
          fmdEsd->SetMultiplicity(det, rng, sec, str, corr * mult);
        }
       }
index 67e8a1f..5d048e7 100644 (file)
@@ -1,7 +1,7 @@
 #ifndef ALIFMDESDREVERTEXER_H
 #define ALIFMDESDREVERTEXER_H
-# include <Rtypes.h>
-class ::AliESDFMD;
+#include <TObject.h>
+class AliESDFMD;
 
 //
 // Class to recaculate quantities in an AliESDFMD object based on new
@@ -11,7 +11,7 @@ class ::AliESDFMD;
 // z-coordinate of the primary interaction, to recalibrate the signals
 // in the FMD ESD object, without having to redo the reconstruction. 
 //
-class AliFMDESDRevertexer 
+class AliFMDESDRevertexer : public TObject
 {
 public:
   /** 
@@ -69,137 +69,9 @@ public:
    */  
   Double_t Eta2Theta(Double_t eta) const;
 protected:
+  ClassDef(AliFMDESDRevertexer,0) // Revertex and FMD ESD Object.
 };
 
-# ifndef __CINT__
-#  ifndef ALIFMDGEOMETRY_H
-#   include <AliFMDGeometry.h>
-#  endif
-#  ifndef ALIESDFMD_H
-#   include <AliESDFMD.h>
-#  endif
-
-//____________________________________________________________________
-inline 
-AliFMDESDRevertexer::AliFMDESDRevertexer()
-{
-  // Constructor 
-  AliFMDGeometry* geom = AliFMDGeometry::Instance();
-  geom->Init();
-  geom->InitTransformations();
-}
-
-//____________________________________________________________________
-inline Bool_t
-AliFMDESDRevertexer::Revertex(AliESDFMD* fmdEsd, Double_t vz) const
-{
-  // Revertex the passed ESD.   The passed ESD object will be modified
-  // directly. 
-  // 
-  // Parameters:
-  //   fmdEsd ESD object to revertex.
-  //   vz     New Z coordinate of primary vertex. 
-  //
-  // Return: 
-  //    kTRUE on success, kFALSE failure.
-  if (!fmdEsd) return kFALSE;
-  
-  Bool_t         ret  = kTRUE;
-  const UShort_t sec0 = 0;
-  for (UShort_t det = 1; det <= 3; det++) { 
-    UShort_t nrng = (det == 1 ? 1 : 2);
-    for (UShort_t ir = 0; ir < nrng; ir++) {
-      Char_t   rng  = (ir == 0 ? 'I' : 'O');
-      UShort_t nsec = (ir == 0 ?  20 :  40);
-      UShort_t nstr = (ir == 0 ? 512 : 256);
-      for (UShort_t str = 0; str < nstr; str++) { 
-       Double_t phi, r, theta;
-       Double_t eta      = AliESDFMD::kInvalidEta;
-       Double_t oldEta   = fmdEsd->Eta(det, rng, 0u, str);
-       Double_t oldTheta = Eta2Theta(oldEta);
-       Bool_t   ret1     = PhysicalCoordinates(det, rng, sec0, str, vz, 
-                                               eta, phi, r, theta);
-       fmdEsd->SetEta(det, rng, sec0, str, eta);
-
-       if (!ret1) {
-         // If the was an error, then there's no reason to go on with
-         // this strip-ring.  Note, that the eta is correctly set to
-         // AliESDFMD::kInvalidMult. 
-         ret = kFALSE;
-         continue;
-       }
-       
-       Double_t corr = 1; 
-       if (fmdEsd->IsAngleCorrected()) 
-         corr = TMath::Abs(TMath::Cos(theta) / TMath::Cos(oldTheta));
-       for (UShort_t sec = 0; sec < nsec; sec++) { 
-         Double_t mult = fmdEsd->Multiplicity(det, rng, sec, str);
-         fmdEsd->SetMultiplicity(det, rng, sec, str, corr * mult);
-       }
-      }
-    }
-  }
-
-  return ret;
-}
-
-//____________________________________________________________________
-inline Double_t
-AliFMDESDRevertexer::Eta2Theta(Double_t eta) const
-{
-  // Calculate the polar angle @f$ \theta@f$ corresponding to the
-  // psuedo-rapidity @f$ \eta@f$ 
-  // 
-  // Parameters:
-  //   eta Psuedo rapidity @f$ \eta=-\log[\tan(\theta/2)]@f$ 
-  // Return:
-  //   Polar angle @f$ \theta=2\atan[\exp(-\eta)]@f$
-  return 2 * TMath::ATan(TMath::Exp(-eta));
-}
-
-
-//____________________________________________________________________
-inline Bool_t
-AliFMDESDRevertexer::PhysicalCoordinates(UShort_t  det, 
-                                        Char_t    rng, 
-                                        UShort_t  sec, 
-                                        UShort_t  str,
-                                        Double_t  vz,
-                                        Double_t& eta, 
-                                        Double_t& phi,
-                                        Double_t& r,
-                                        Double_t& theta) const
-{
-  // Calculate the physical coordinates (@a eta, @a phi) corresponding
-  // to the detector coordinates (@a det, @a rng, @a sec, @a str).
-  // 
-  // Parameters:
-  //   det  The detector identifier 
-  //   rng  The ring identifier 
-  //   sec  The sector identifier 
-  //   str  The strip identifier 
-  //   vz   The z coordinate of the current primary interation vertex
-  //   eta  On return, the psuedo-rapidity
-  //   phi  On return, the azimuthal angle
-  //   r    On return, the radius
-  //   phi  On return, the polar angle
-  AliFMDGeometry* geom = AliFMDGeometry::Instance();
-  Double_t x=0, y=0, z=0;
-  geom->Detector2XYZ(det, rng, sec, str, x, y, z);
-
-  // Check that the conversion succeeded
-  if (x == 0 && y == 0 && z == 0) return kFALSE;
-  
-  // Correct for vertex offset. 
-  z     -= vz;
-  phi   =  TMath::ATan2(y, x);
-  r     =  TMath::Sqrt(y * y + x * x);
-  theta =  TMath::ATan2(r, z);
-  eta   = -TMath::Log(TMath::Tan(theta / 2));
-
-  return kTRUE;
-}
-# endif // __CINT__
 #endif
 //
 // Local Variables:
index 0206986..08260c0 100644 (file)
@@ -251,9 +251,9 @@ AliFMDRawReader::NewSample(AliAltroRawStreamV3& input,
   map->Timebin2Strip(sec, t, fPreSamp, fSampleRate[ddl], stroff, samp);
   str             = strbase + stroff;
       
-  AliFMDDebug(20, ("0x%04x/0x%03x/%04d maps to strip %3d sample %d " 
+  AliFMDDebug(20, ("0x%04x/0x%03x/%04d=%4d maps to strip %3d sample %d " 
                   "(pre: %d, min: %d, max: %d, rate: %d)",
-                 ddl, hwa, t, str, samp, fPreSamp, 
+                 ddl, hwa, t, adc, str, samp, fPreSamp, 
                  fMinStrip, fMaxStrip, fSampleRate[ddl]));
   if (str < 0) { 
     AliFMDDebug(10, ("Got presamples at timebin %d", i));
@@ -266,9 +266,9 @@ AliFMDRawReader::NewSample(AliAltroRawStreamV3& input,
   AliFMDDebug(15, ("Checking if strip %d (%d) in range [%d,%d]", 
                   lstrip, str, fMinStrip, fMaxStrip));
   if (lstrip < fMinStrip || lstrip > fMaxStrip) {
-    AliFMDDebug(20, ("Strip %03d-%d (%d,%d) from t=%d out of range (%3d->%3d)", 
+    AliFMDDebug(10, ("Strip %03d-%d (%d,%d) from t=%d out of range (%3d->%3d)", 
                    str, samp, lstrip, stroff, t, fMinStrip, fMaxStrip));
-    return -1;
+    adc = -1;
   }
   // Possibly do pedestal subtraction of signal 
   if (adc > 1023) 
@@ -450,16 +450,18 @@ AliFMDRawReader::ReadAdcs(TClonesArray* array)
                           det, ring, sec, str, samp, counts));
          // Check the cache of indicies
          Int_t idx = fSeen(det, ring, sec, str);
+         AliFMDDigit* digit = 0;
          if (idx == kUShortMax) { 
            // We haven't seen this strip yet. 
            fSeen(det, ring, sec, str) = idx = array->GetEntriesFast();
            AliFMDDebug(7,("making digit for FMD%d%c[%2d,%3d]-%d "
                           "from timebin %4d", 
                           det, ring, sec, str, samp, t));
-           new ((*array)[idx]) AliFMDDigit(det, ring, sec, str);
+           digit = new ((*array)[idx]) AliFMDDigit(det, ring, sec, str);
+           digit->SetDefaultCounts(fSampleRate[ddl]);
          }
-         AliFMDBaseDigit* digit = 
-           static_cast<AliFMDBaseDigit*>(array->At(idx));
+         else 
+           digit = static_cast<AliFMDDigit*>(array->At(idx));
          AliFMDDebug(10, ("Setting FMD%d%c[%2d,%3d]-%d "
                           "from timebin %4d=%4d (%4d)", 
                           det, ring, sec, str, samp, t, counts, data[i]));
@@ -490,8 +492,14 @@ Bool_t AliFMDRawReader::ReadSODevent(AliFMDCalibSampleRate* sampleRate,
   AliFMDParameters*   param = AliFMDParameters::Instance();
   AliFMDAltroMapping* map   = param->GetAltroMap();
   
-  while(fReader->ReadNextData(fData)) {
-    
+  AliAltroRawStreamV3  streamer(fReader);
+  streamer.Reset();
+  streamer.SelectRawData("FMD");
+  //fReader->GetDDLID();
+  //Int_t detID = fReader->GetDetectorID();
+  
+  //  while(fReader->ReadNextData(fData)) {
+  /*  
     Int_t ddl   = fReader->GetDDLID();
     Int_t detID = fReader->GetDetectorID();
     if (detectors) detectors[map->DDL2Detector(ddl)-1] = kTRUE;
@@ -516,10 +524,28 @@ Bool_t AliFMDRawReader::ReadSODevent(AliFMDCalibSampleRate* sampleRate,
     AliFMDDebug(20, (" # trailer words: %d, # payload words: %d", 
                     nTrailerWords, nPayloadWords));
     
+    ULong_t nPayloadWords = streamer.GetSOD...();
+    UInt_t   payloadWords* = streamer.GetSOD...();
+  */
+
+  
+  
+  while (streamer.NextDDL()) {
+    Int_t ddl   = streamer.GetDDLNumber();
+    Int_t detID = fReader->GetDetectorID();
+    if (detectors) detectors[map->DDL2Detector(ddl)-1] = kTRUE;
+    AliFMDDebug(0, (" From reader: DDL number is %d , det ID is %d",ddl,detID));
+
+    ULong_t  nPayloadWords = streamer.GetRCUPayloadSizeInSOD();
+    UChar_t* payloadData   = streamer.GetRCUPayloadInSOD();
+    UInt_t*  payloadWords  = reinterpret_cast<UInt_t*>(payloadData);
+    //UInt_t*   payloadWords  = streamer.GetRCUPayloadInSOD();
+
+    std::cout<<nPayloadWords<<"    "<<ddl<<std::endl;
+    for (ULong_t i = 1; i <= nPayloadWords ; i++, payloadWords++) {
+      UInt_t payloadWord = *payloadWords; // Get32bitWord(i);
     
-    for (ULong_t i = 1; i <= nPayloadWords ; i++) {
-      UInt_t payloadWord = Get32bitWord(i);
-      
+      std::cout<<i<<Form("  word: 0x%x",payloadWord)<<std::endl;
       // address is only 24 bit
       UInt_t address       = (0xffffff & payloadWord);
       UInt_t type          = ((address >> 21) & 0xf);
@@ -544,14 +570,15 @@ Bool_t AliFMDRawReader::ReadSODevent(AliFMDCalibSampleRate* sampleRate,
        readDataWord = kTRUE;  
       case 0x1: // Fec cmd
       case 0x2: // Fec write
-       i++;  
+               i++;  
+       payloadWords++;
        break;
       case 0x4: // Loop
       case 0x5: // Wait
        break;
       case 0x6: // End sequence
       case 0x7: // End Mem
-       i = nPayloadWords + 1;
+               i = nPayloadWords + 1;
        break;
       default:    
        break;
@@ -560,7 +587,7 @@ Bool_t AliFMDRawReader::ReadSODevent(AliFMDCalibSampleRate* sampleRate,
       //Don't read unless we have a FEC_RD
       if(!readDataWord)  continue;
 
-      UInt_t dataWord      = Get32bitWord(i);
+      UInt_t dataWord      = *payloadWords;//Get32bitWord(i);
       UInt_t data          = (0xFFFFF & dataWord) ;
       //UInt_t data          = (0xFFFF & dataWord) ;
        
@@ -674,7 +701,7 @@ Bool_t AliFMDRawReader::ReadSODevent(AliFMDCalibSampleRate* sampleRate,
       }
       AliFMDDebug(50, ("instruction 0x%x, dataword 0x%x",
                       instruction,dataWord));
-    }
+    } // End of loop over Result memory event
     
     UShort_t det,sector;
     Short_t strip;
@@ -747,9 +774,9 @@ Bool_t AliFMDRawReader::ReadSODevent(AliFMDCalibSampleRate* sampleRate,
                       strip_low[boards[i]], strip_high[boards[i]],
                       shift_clk[boards[i]], sample_clk[boards[i]],
                       pulse_size[boards[i]],pulse_length[boards[i]]));
-    }
+      }
     
-  }
+    }
   
   AliFMDParameters::Instance()->SetSampleRate(sampleRate);
   AliFMDParameters::Instance()->SetStripRange(stripRange);
index f6493fd..ec353a0 100644 (file)
 #include <TH2.h>
 #include <TFile.h>
 #include <climits>
-// Import revertexer into a private namespace (to prevent conflicts) 
-namespace { 
-# include "AliFMDESDRevertexer.h"
-}
+#include "AliFMDESDRevertexer.h"
 
 
 class AliRawReader;
@@ -321,13 +318,14 @@ AliFMDReconstructor::Reconstruct(TTree* digitsTree,
   // 
   AliFMDDebug(1, ("Reconstructing from digits in a tree"));
   GetVertex(fESD);
-  
+
+  static TClonesArray* digits = new TClonesArray("AliFMDDigit");
   TBranch *digitBranch = digitsTree->GetBranch("FMD");
   if (!digitBranch) {
     Error("Exec", "No digit branch for the FMD found");
     return;
   }
-  TClonesArray* digits = new TClonesArray("AliFMDDigit");
+  // TClonesArray* digits = new TClonesArray("AliFMDDigit");
   digitBranch->SetAddress(&digits);
 
   if (fMult)   fMult->Clear();
@@ -348,7 +346,7 @@ AliFMDReconstructor::Reconstruct(TTree* digitsTree,
   Int_t written = clusterTree->Fill();
   AliFMDDebug(10, ("Filled %d bytes into cluster tree", written));
   digits->Delete();
-  delete digits;
+  // delete digits;
 }
  
 
@@ -383,7 +381,8 @@ AliFMDReconstructor::ProcessDigit(AliFMDDigit* digit) const
   UShort_t sec = digit->Sector();
   UShort_t str = digit->Strip();
   Short_t  adc = digit->Counts();
-  
+  if (adc < 0) 
+    digit->Print();
   ProcessSignal(det, rng, sec, str, adc);
 }
 
@@ -407,7 +406,7 @@ AliFMDReconstructor::ProcessSignal(UShort_t det,
   AliFMDParameters* param  = AliFMDParameters::Instance();
   // Check that the strip is not marked as dead 
   if (param->IsDead(det, rng, sec, str)) {
-    AliFMDDebug(3, ("FMD%d%c[%2d,%3d] is dead", det, rng, sec, str));
+    AliFMDDebug(1, ("FMD%d%c[%2d,%3d] is dead", det, rng, sec, str));
     return;
   }
   
@@ -428,6 +427,11 @@ AliFMDReconstructor::ProcessSignal(UShort_t det,
   // Make rough multiplicity 
   Double_t mult     = Energy2Multiplicity(det, rng, sec, str, edep);
   // Get rid of nonsense mult
+  if (mult > 20) { 
+    AliWarning(Form("The mutliplicity in FMD%d%c[%2d,%3d]=%f > 20 "
+                   "(ADC: %d, Energy: %f)", det, rng, sec, str, mult, 
+                   counts, edep));
+  }
   if (mult < 0)  return; 
   AliFMDDebug(10, ("FMD%d%c[%2d,%3d]: "
                    "ADC: %d, Counts: %d, Energy: %f, Mult: %f",
@@ -540,7 +544,8 @@ AliFMDReconstructor::SubtractPedestal(UShort_t det,
   if (counts < noise * nf) counts = 0;
   if (counts > 0) AliDebugClass(15, "Got a hit strip");
 
-  return counts;
+  UShort_t ret = counts < 0 ? 0 : counts;
+  return ret;
 }
 
 
index 00f56d0..bf6ee9d 100644 (file)
@@ -26,6 +26,7 @@
 #pragma link C++ class  AliFMDRawReader+;
 #pragma link C++ class  AliFMDQADataMakerRec+;
 #pragma link C++ class  AliFMDOfflineTrigger+;
+#pragma link C++ class  AliFMDESDRevertexer+;
 
 #else
 # error Not for compilation 
index bab5842..0505030 100644 (file)
@@ -8,7 +8,8 @@ SRCS            =  AliFMDReconstructor.cxx      \
                   AliFMDRawReader.cxx          \
                   AliFMDRecPoint.cxx           \
                   AliFMDQADataMakerRec.cxx     \
-                  AliFMDOfflineTrigger.cxx
+                  AliFMDOfflineTrigger.cxx     \
+                  AliFMDESDRevertexer.cxx      
 
 HDRS           =  $(SRCS:.cxx=.h)
 DHDR           := FMDrecLinkDef.h