]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - MUON/AliMUONResponseV0.cxx
Coding conventions
[u/mrichter/AliRoot.git] / MUON / AliMUONResponseV0.cxx
index a820fcff52410b5c44ae5e7caed4a1b549e704c0..1828361019df05e36c352dd70f4277dee072b513 100644 (file)
 
 /* $Id$ */
 
-#include <TMath.h>
-#include <TRandom.h>
+// --------------------------
+// Class AliMUONResponseV0
+// --------------------------
+// Implementation of Mathieson response
+// ...
 
 #include "AliMUONResponseV0.h"
-#include "AliMUONGeometrySegmentation.h"
+
 #include "AliLog.h"
+#include "AliMUON.h"
+#include "AliMUONConstants.h"
+#include "AliMUONDigit.h"
+#include "AliMUONGeometrySegmentation.h"
+#include "AliMUONGeometryTransformer.h"
+#include "AliMUONHit.h"
+#include "AliMUONSegmentation.h"
+#include "AliMpArea.h"
+#include "AliMpDEManager.h"
+#include "AliMpVPadIterator.h"
+#include "AliMpVSegmentation.h"
+#include "AliRun.h"
+#include "Riostream.h"
+#include "TVector2.h"
+#include <TMath.h>
+#include <TRandom.h>
 
 ClassImp(AliMUONResponseV0)
        
+AliMUON* muon()
+{
+    return static_cast<AliMUON*>(gAlice->GetModule("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.
+  
+  const AliMUONGeometryTransformer* transformer = muon()->GetGeometryTransformer();
+  transformer->Global2Local(detElemId,xg,yg,zg,xl,yl,zl);
+}
+
+AliMUONSegmentation* Segmentation()
+{
+  static AliMUONSegmentation* segmentation = muon()->GetSegmentation();
+  return segmentation;
+}
+
 //__________________________________________________________________________
 AliMUONResponseV0::AliMUONResponseV0()
-  : AliMUONResponse()
+  : AliMUONResponse(),
+  fChargeSlope(0.0),
+  fChargeSpreadX(0.0),
+  fChargeSpreadY(0.0),
+  fSigmaIntegration(0.0),
+  fMaxAdc(0),
+  fZeroSuppression(0),
+  fChargeCorrel(0.0),
+  fMathieson(new AliMUONMathieson),
+  fChargeThreshold(1e-4)
 {
-// Default constructor
-
-  fMathieson = new AliMUONMathieson();
-  fChargeCorrel = 0;
+    // Normal constructor
+    AliDebug(1,Form("Default ctor"));
 }
 
    //_________________________________________________________________________
@@ -46,6 +96,7 @@ AliMUONResponseV0::AliMUONResponseV0(const AliMUONResponseV0& rhs)
    //__________________________________________________________________________
 AliMUONResponseV0::~AliMUONResponseV0()
 {
+  AliDebug(1,"");
   delete fMathieson;
 }
 
@@ -59,6 +110,33 @@ AliMUONResponseV0& AliMUONResponseV0::operator = (const AliMUONResponseV0& rhs)
   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  
+//
+
 }
 
   //__________________________________________________________________________
@@ -103,28 +181,151 @@ Float_t AliMUONResponseV0::IntXY(Int_t idDE, AliMUONGeometrySegmentation* segmen
 
   return fMathieson->IntXY(idDE, segmentation);
 }
+
+
   //-------------------------------------------
 Int_t  AliMUONResponseV0::DigitResponse(Int_t digit, AliMUONTransientDigit* /*where*/)
 {
-    // add white noise and do zero-suppression and signal truncation
-//     Float_t meanNoise = gRandom->Gaus(1, 0.2);
+  // \deprecated method
+  // Now part of the digitizer (where it belongs really), e.g. DigitizerV3
+  //
+  // 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 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();
+    if ( digit >  Saturation())          
+    {
+      digit=Saturation();
+    }
 
     return digit;
 }
 
+//_____________________________________________________________________________
+Float_t
+AliMUONResponseV0::GetAnod(Float_t x) const
+{
+  //
+  // 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;
+}
 
+//______________________________________________________________________________
+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();
+  
+  // Use that (dx,dy) to specify the area upon which
+  // we will iterate to spread charge into.
+  Double_t x,y,z;
+  Global2Local(detElemId,hit.X(),hit.Y(),hit.Z(),x,y,z);
+  x = GetAnod(x);
+  TVector2 hitPosition(x,y);
+  AliMpArea area(hitPosition,TVector2(dx,dy));
+  
+  // 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));
 
-
-
-
+  for ( Int_t cath = 0; cath < 2; ++cath )
+  {
+    Float_t qcath = qtot * ( cath == 0 ? currentCorrel : 1.0/currentCorrel);
+    
+    // Get an iterator to loop over pads, within the given area.
+    const AliMpVSegmentation* seg = 
+        Segmentation()->GetMpSegmentation(detElemId,cath);
+      
+    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() )
+      {
+        AliWarning(Form("Got an invalid iterator bug (area.Position() is within "
+                      " DE but the iterator is void) for detElemId %d cath %d",
+                      detElemId,cath));        
+      }
+      else
+      {
+        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;
+  }
+}