]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - MUON/AliMUONResponseV0.cxx
cleaning (Christian)
[u/mrichter/AliRoot.git] / MUON / AliMUONResponseV0.cxx
index 7f0a255af21762b5b56845281ddaa6cbb4ed71a7..81747d63e50a2ea4ca2f0c35ab375792228de70a 100644 (file)
 
 /* $Id$ */
 
+// --------------------------
+// Class AliMUONResponseV0
+// --------------------------
+// Implementation of 
+// Mathieson response
 
 #include "AliMUONResponseV0.h"
-
-#include "AliLog.h"
 #include "AliMUON.h"
 #include "AliMUONConstants.h"
 #include "AliMUONDigit.h"
 #include "AliMUONGeometryTransformer.h"
 #include "AliMUONHit.h"
 #include "AliMUONSegmentation.h"
+
 #include "AliMpArea.h"
 #include "AliMpDEManager.h"
-#include "AliMpPlaneType.h"
-#include "AliMpStationType.h"
 #include "AliMpVPadIterator.h"
+#include "AliMpSegmentation.h"
 #include "AliMpVSegmentation.h"
+
 #include "AliRun.h"
+#include "AliLog.h"
+
 #include "Riostream.h"
 #include "TVector2.h"
 #include <TMath.h>
 #include <TRandom.h>
 
+/// \cond CLASSIMP
 ClassImp(AliMUONResponseV0)
+/// \endcond
        
 AliMUON* muon()
 {
@@ -48,11 +56,11 @@ AliMUON* muon()
 void Global2Local(Int_t detElemId, Double_t xg, Double_t yg, Double_t zg,
                   Double_t& xl, Double_t& yl, Double_t& zl)
 {  
-  // ideally should be : 
-  // Double_t x,y,z;
-  // AliMUONGeometry::Global2Local(detElemId,xg,yg,zg,x,y,z);
-  // but while waiting for this geometry singleton, let's go through
-  // AliMUON still.
+  /// ideally should be : 
+  /// Double_t x,y,z;
+  /// AliMUONGeometry::Global2Local(detElemId,xg,yg,zg,x,y,z);
+  /// but while waiting for this geometry singleton, let's go through
+  /// AliMUON still.
   
   const AliMUONGeometryTransformer* transformer = muon()->GetGeometryTransformer();
   transformer->Global2Local(detElemId,xg,yg,zg,xl,yl,zl);
@@ -72,91 +80,61 @@ AliMUONResponseV0::AliMUONResponseV0()
   fChargeSpreadY(0.0),
   fSigmaIntegration(0.0),
   fMaxAdc(0),
+  fSaturation(0),
   fZeroSuppression(0),
   fChargeCorrel(0.0),
   fMathieson(new AliMUONMathieson),
   fChargeThreshold(1e-4)
 {
-    // Normal constructor
+    /// Normal constructor
     AliDebug(1,Form("Default ctor"));
 }
 
-   //_________________________________________________________________________
-AliMUONResponseV0::AliMUONResponseV0(const AliMUONResponseV0& rhs)
- : AliMUONResponse(rhs)
-{
-// Protected copy constructor
-
-  AliFatal("Not implemented.");
-}
-
-   //__________________________________________________________________________
+//__________________________________________________________________________
 AliMUONResponseV0::~AliMUONResponseV0()
 {
+/// Destructor
+
   AliDebug(1,"");
   delete fMathieson;
 }
 
-   //________________________________________________________________________
-AliMUONResponseV0& AliMUONResponseV0::operator = (const AliMUONResponseV0& rhs)
-{
-// Protected assignement operator
-
-  if (this == &rhs) return *this;
-
-  AliFatal("Not implemented.");
-    
-  return *this;  
-}
-
 //______________________________________________________________________________
 void
 AliMUONResponseV0::Print(Option_t*) const
 {
+/// Printing
+
   cout << " ChargeSlope=" << fChargeSlope
     << " ChargeSpreadX,Y=" << fChargeSpreadX
     << fChargeSpreadY
     << " ChargeCorrelation=" << fChargeCorrel
     << endl;
-  
-//Float_t fChargeSlope;              // Slope of the charge distribution
-//Float_t fChargeSpreadX;            // Width of the charge distribution in x
-//Float_t fChargeSpreadY;            // Width of the charge distribution in y
-//Float_t fSigmaIntegration;         // Number of sigma's used for charge distribution
-//Int_t   fMaxAdc;                   // Maximum ADC channel
-//Int_t   fSaturation;               // Pad saturation in ADC channel
-//Int_t   fZeroSuppression;          // Zero suppression threshold
-//Float_t fChargeCorrel;             // amplitude of charge correlation on 2 cathods
-//                                   // is RMS of ln(q1/q2)
-//AliMUONMathieson* fMathieson;      // pointer to mathieson fct
-//Float_t fChargeThreshold;          // Charges below this threshold are = 0  
-//
-
 }
 
   //__________________________________________________________________________
 void AliMUONResponseV0::SetSqrtKx3AndDeriveKx2Kx4(Float_t SqrtKx3)
 {
-  // Set to "SqrtKx3" the Mathieson parameter K3 ("fSqrtKx3")
-  // in the X direction, perpendicular to the wires,
-  // and derive the Mathieson parameters K2 ("fKx2") and K4 ("fKx4")
-  // in the same direction
+  /// Set to "SqrtKx3" the Mathieson parameter K3 ("fSqrtKx3")
+  /// in the X direction, perpendicular to the wires,
+  /// and derive the Mathieson parameters K2 ("fKx2") and K4 ("fKx4")
+  /// in the same direction
   fMathieson->SetSqrtKx3AndDeriveKx2Kx4(SqrtKx3);
 }
        
   //__________________________________________________________________________
 void AliMUONResponseV0::SetSqrtKy3AndDeriveKy2Ky4(Float_t SqrtKy3)
 {
-  // Set to "SqrtKy3" the Mathieson parameter K3 ("fSqrtKy3")
-  // in the Y direction, along the wires,
-  // and derive the Mathieson parameters K2 ("fKy2") and K4 ("fKy4")
-  // in the same direction
+  /// Set to "SqrtKy3" the Mathieson parameter K3 ("fSqrtKy3")
+  /// in the Y direction, along the wires,
+  /// and derive the Mathieson parameters K2 ("fKy2") and K4 ("fKy4")
+  /// in the same direction
   fMathieson->SetSqrtKy3AndDeriveKy2Ky4(SqrtKy3);
 }
   //__________________________________________________________________________
-Float_t AliMUONResponseV0::IntPH(Float_t eloss)
+Float_t AliMUONResponseV0::IntPH(Float_t eloss) const
 {
-  // Calculate charge from given ionization energy loss
+  /// Calculate charge from given ionization energy loss
   Int_t nel;
   nel= Int_t(eloss*1.e9/27.4);
   Float_t charge=0;
@@ -170,43 +148,22 @@ Float_t AliMUONResponseV0::IntPH(Float_t eloss)
 }
 
   //-------------------------------------------
-Float_t AliMUONResponseV0::IntXY(Int_t idDE, AliMUONGeometrySegmentation* segmentation)
+Float_t AliMUONResponseV0::IntXY(Int_t idDE,
+                                AliMUONGeometrySegmentation* segmentation) 
+const
 {
- // Calculate charge on current pad according to Mathieson distribution
+ /// Calculate charge on current pad according to Mathieson distribution
 
   return fMathieson->IntXY(idDE, segmentation);
 }
 
 
-  //-------------------------------------------
-Int_t  AliMUONResponseV0::DigitResponse(Int_t digit, AliMUONTransientDigit* /*where*/)
-{
-//  FIXME : AliFatal("put the pedestal adding here!");
-    // add white noise and do zero-suppression and signal truncation
-//     Float_t meanNoise = gRandom->Gaus(1, 0.2);
-    // correct noise for slat chambers;
-    // one more field to add to AliMUONResponseV0 to allow different noises ????
-//    Float_t meanNoise = gRandom->Gaus(1., 0.2);
-//    Float_t noise     = gRandom->Gaus(0., meanNoise);
-    Float_t noise     = gRandom->Gaus(0., 1.0);
-    digit += TMath::Nint(noise); 
-    if ( digit <= ZeroSuppression()) digit = 0;
-    // if ( digit >  MaxAdc())          digit=MaxAdc();
-    if ( digit >  Saturation())          
-    {
-      digit=Saturation();
-    }
-
-    return digit;
-}
-
 //_____________________________________________________________________________
 Float_t
 AliMUONResponseV0::GetAnod(Float_t x) const
 {
-  //
-  // Return wire coordinate closest to x.
-  //
+  /// Return wire coordinate closest to x.
+
   Int_t n = Int_t(x/Pitch());
   Float_t wire = (x>0) ? n+0.5 : n-0.5;
   return Pitch()*wire;
@@ -216,17 +173,16 @@ AliMUONResponseV0::GetAnod(Float_t x) const
 void 
 AliMUONResponseV0::DisIntegrate(const AliMUONHit& hit, TList& digits)
 {
-  //
-  //
-  //
+  /// Go from 1 hit to a list of digits.
+  /// The energy deposition of that hit is first converted into charge
+  /// (in IntPH() method), and then this charge is dispatched on several
+  /// pads, according to the Mathieson distribution.
   
   digits.Clear();
   
   Int_t detElemId = hit.DetElemId();
   
-  //
   // Width of the integration area
-  //
   Double_t dx = SigmaIntegration()*ChargeSpreadX();
   Double_t dy = SigmaIntegration()*ChargeSpreadY();
   
@@ -238,80 +194,83 @@ AliMUONResponseV0::DisIntegrate(const AliMUONHit& hit, TList& digits)
   TVector2 hitPosition(x,y);
   AliMpArea area(hitPosition,TVector2(dx,dy));
   
-  //
-  // Get pulse height from energy loss
+  // Get pulse height from energy loss.
   Float_t qtot = IntPH(hit.Eloss());
   
+  // Get the charge correlation between cathodes.
   Float_t currentCorrel = TMath::Exp(gRandom->Gaus(0.0,ChargeCorrel()/2.0));
 
-  AliDebug(4,Form("DE=%d eloss=%e x,y,z=%e,%e,%e fCurrentCorrel=%e dx,dy=%e,%e",
-                  detElemId,hit.Eloss(),hit.X(),hit.Y(),hit.Z(),
-                  currentCorrel,dx,dy));
-      
-  AliMpStationType station = AliMpDEManager::GetStationType(detElemId);
-  
-  Int_t intOffset = 1;
-  if ( station == kStation1 || station == kStation2 )
-  {
-    intOffset = 0;
-  }
-  
   for ( Int_t cath = 0; cath < 2; ++cath )
   {
     Float_t qcath = qtot * ( cath == 0 ? currentCorrel : 1.0/currentCorrel);
     
-    AliDebug(4,Form("i=%d qtot=%e qcath=%e",cath+1,qtot,qcath));      
     // Get an iterator to loop over pads, within the given area.
-      const AliMpVSegmentation* seg = 
-        Segmentation()->GetMpSegmentation(detElemId,cath);
-      AliMpVPadIterator* it = seg->CreateIterator(area);
+    const AliMpVSegmentation* seg = 
+        AliMpSegmentation::Instance()->GetMpSegmentation(detElemId,cath);
       
-      if (!it)
+    AliMpVPadIterator* it = seg->CreateIterator(area);
+      
+    if (!it)
+    {
+      AliError(Form("Could not get iterator for detElemId %d",detElemId));
+      return;
+    }
+    
+    // Start loop over pads.
+    it->First();
+    
+    if ( it->IsDone() )
+    {
+      // Exceptional case : iterator is built, but is invalid from the start.
+      AliMpPad pad = seg->PadByPosition(area.Position(),kFALSE);
+      if ( pad.IsValid() )
       {
-        AliError(Form("Could not get iterator for detElemId %d",detElemId));
-        return;
+        AliWarning(Form("Got an invalid iterator bug (area.Position() is within "
+                      " DE but the iterator is void) for detElemId %d cath %d",
+                      detElemId,cath));        
       }
-      
-      // Start loop over pads.
-      it->First();
-      while ( !it->IsDone() )
+      else
       {
-        AliMpPad pad = it->CurrentItem();      
-        TVector2 lowerLeft(hitPosition-pad.Position()-pad.Dimensions());
-        TVector2 upperRight(lowerLeft + pad.Dimensions()*2.0);
-        Float_t qp = TMath::Abs(fMathieson->IntXY(lowerLeft.X(),lowerLeft.Y(),
-                                                  upperRight.X(),upperRight.Y()));
-
-
-        Int_t icharge = Int_t(qp*qcath);
-        
-        if ( qp > fChargeThreshold )
-        {
-          AliDebug(4,Form("ix,iy=%d,%d qp=%e",
-                          pad.GetIndices().GetFirst()+intOffset,
-                          pad.GetIndices().GetSecond()+intOffset,
-                          qp));
-          AliMUONDigit* d = new AliMUONDigit;
-          d->SetDetElemId(detElemId);
-          d->SetPadX(pad.GetIndices().GetFirst()+intOffset);
-          d->SetPadY(pad.GetIndices().GetSecond()+intOffset);
-          d->SetSignal(icharge);
-          d->AddPhysicsSignal(d->Signal());
-          d->SetCathode(cath);
-          Int_t manuId = pad.GetLocation().GetFirst();
-          Int_t manuChannel = pad.GetLocation().GetSecond();
-          AliMpPlaneType planeType = AliMpDEManager::GetPlaneType(detElemId,cath);
-          if ( planeType == kNonBendingPlane )
-          {
-            // FIXME: this should not be there, but integrated in the mapping files directly.
-            manuId |= (1<<11);
-          }
-          d->SetElectronics(manuId,manuChannel);
-          digits.Add(d);   
-        }       
-        it->Next();
+        AliError(Form("Got an invalid iterator bug for detElemId %d cath %d."
+                      "Might be a bad hit ? area.Position()=(%e,%e) "
+                      "Dimensions()=(%e,%e)",
+                      detElemId,cath,area.Position().X(),area.Position().Y(),
+                      area.Dimensions().X(),area.Dimensions().Y()));
       }
       delete it;
+      return;
+    }
+    
+    while ( !it->IsDone() )
+    {
+      // For each pad given by the iterator, compute the charge of that
+      // pad, according to the Mathieson distribution.
+      AliMpPad pad = it->CurrentItem();      
+      TVector2 lowerLeft(hitPosition-pad.Position()-pad.Dimensions());
+      TVector2 upperRight(lowerLeft + pad.Dimensions()*2.0);
+      Float_t qp = TMath::Abs(fMathieson->IntXY(lowerLeft.X(),lowerLeft.Y(),
+                                                upperRight.X(),upperRight.Y()));
+            
+      Int_t icharge = Int_t(qp*qcath);
+      
+      if ( qp > fChargeThreshold )
+      {
+        // If we're above threshold, then we create a digit,
+        // and fill it with relevant information, including electronics.
+        AliMUONDigit* d = new AliMUONDigit;
+        d->SetDetElemId(detElemId);
+        d->SetPadX(pad.GetIndices().GetFirst());
+        d->SetPadY(pad.GetIndices().GetSecond());
+        d->SetSignal(icharge);
+        d->AddPhysicsSignal(d->Signal());
+        d->SetCathode(cath);
+        d->SetElectronics(pad.GetLocation().GetFirst(),
+                          pad.GetLocation().GetSecond());
+        digits.Add(d);   
+      }       
+      it->Next();
+    }
+    delete it;
   }
 }