/* $Id$ */
-#include <TMath.h>
-#include <TRandom.h>
+//-----------------------------------------------------------------------------
+// Class AliMUONResponseV0
+// --------------------------
+// Implementation of
+// Mathieson response
+//-----------------------------------------------------------------------------
#include "AliMUONResponseV0.h"
-#include "AliSegmentation.h"
+#include "AliMUON.h"
+#include "AliMUONConstants.h"
+#include "AliMUONDigit.h"
+#include "AliMUONGeometryTransformer.h"
+#include "AliMUONHit.h"
+#include "AliMUONConstants.h"
+
+#include "AliMpArea.h"
+#include "AliMpDEManager.h"
+#include "AliMpVPadIterator.h"
+#include "AliMpSegmentation.h"
+#include "AliMpVSegmentation.h"
+#include "AliMpCathodType.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()
+{
+ 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);
+}
+
//__________________________________________________________________________
AliMUONResponseV0::AliMUONResponseV0()
- : AliMUONResponse()
+ : AliMUONResponse(),
+ fChargeSlope(0.0),
+ fChargeSpreadX(0.0),
+ fChargeSpreadY(0.0),
+ fSigmaIntegration(0.0),
+ fMaxAdc(0),
+ fSaturation(0),
+ fZeroSuppression(0),
+ fChargeCorrel(0.0),
+ fMathieson(new AliMUONMathieson),
+ fChargeThreshold(1e-4),
+ fIsTailEffect(kFALSE)
{
-// Default constructor
+ /// Normal constructor
+ AliDebug(1,Form("Default ctor"));
+}
- fChargeCorrel = 0;
+//__________________________________________________________________________
+AliMUONResponseV0::AliMUONResponseV0(const AliMUONResponseV0& other)
+: AliMUONResponse(),
+fChargeSlope(0.0),
+fChargeSpreadX(0.0),
+fChargeSpreadY(0.0),
+fSigmaIntegration(0.0),
+fMaxAdc(0),
+fSaturation(0),
+fZeroSuppression(0),
+fChargeCorrel(0.0),
+fMathieson(0),
+fChargeThreshold(1e-4),
+fIsTailEffect(kFALSE)
+{
+ /// copy ctor
+ other.CopyTo(*this);
+}
+
+//__________________________________________________________________________
+AliMUONResponseV0&
+AliMUONResponseV0::operator=(const AliMUONResponseV0& other)
+{
+ /// Assignment operator
+ other.CopyTo(*this);
+ return *this;
+}
+
+//__________________________________________________________________________
+AliMUONResponseV0::~AliMUONResponseV0()
+{
+/// Destructor
+
+ AliDebug(1,"");
+ delete fMathieson;
+}
+
+//______________________________________________________________________________
+void
+AliMUONResponseV0::CopyTo(AliMUONResponseV0& other) const
+{
+ /// Copy *this to other
+ other.fChargeSlope=fChargeSlope;
+ other.fChargeSpreadX=fChargeSpreadX;
+ other.fChargeSpreadY=fChargeSpreadY;
+ other.fSigmaIntegration=fSigmaIntegration;
+ other.fMaxAdc=fMaxAdc;
+ other.fSaturation=fSaturation;
+ other.fZeroSuppression=fZeroSuppression;
+ other.fChargeCorrel=fChargeCorrel;
+ delete other.fMathieson;
+ other.fMathieson = new AliMUONMathieson(*fMathieson);
+ other.fChargeThreshold=fChargeThreshold;
+}
+
+//______________________________________________________________________________
+void
+AliMUONResponseV0::Print(Option_t*) const
+{
+/// Printing
+
+ cout << " ChargeSlope=" << fChargeSlope
+ << " ChargeSpreadX,Y=" << fChargeSpreadX
+ << fChargeSpreadY
+ << " ChargeCorrelation=" << fChargeCorrel
+ << endl;
}
//__________________________________________________________________________
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
- fSqrtKx3 = SqrtKx3;
- fKx2 = TMath::Pi() / 2. * (1. - 0.5 * fSqrtKx3);
- Float_t cx1 = fKx2 * fSqrtKx3 / 4. / TMath::ATan(Double_t(fSqrtKx3));
- fKx4 = cx1 / fKx2 / fSqrtKx3;
+ /// 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
- fSqrtKy3 = SqrtKy3;
- fKy2 = TMath::Pi() / 2. * (1. - 0.5 * fSqrtKy3);
- Float_t cy1 = fKy2 * fSqrtKy3 / 4. / TMath::ATan(Double_t(fSqrtKy3));
- fKy4 = cy1 / fKy2 / fSqrtKy3;
+ /// 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;
}
return charge;
}
-// -------------------------------------------
-Float_t AliMUONResponseV0::IntXY(AliSegmentation * segmentation)
+//_____________________________________________________________________________
+Float_t
+AliMUONResponseV0::GetAnod(Float_t x) const
{
-// Calculate charge on current pad according to Mathieson distribution
-//
- const Float_t kInversePitch = 1/fPitch;
-//
-// Integration limits defined by segmentation model
-//
- Float_t xi1, xi2, yi1, yi2;
- segmentation->IntegrationLimits(xi1,xi2,yi1,yi2);
- xi1=xi1*kInversePitch;
- xi2=xi2*kInversePitch;
- yi1=yi1*kInversePitch;
- yi2=yi2*kInversePitch;
-//
-// The Mathieson function
- Double_t ux1=fSqrtKx3*TMath::TanH(fKx2*xi1);
- Double_t ux2=fSqrtKx3*TMath::TanH(fKx2*xi2);
-
- Double_t uy1=fSqrtKy3*TMath::TanH(fKy2*yi1);
- Double_t uy2=fSqrtKy3*TMath::TanH(fKy2*yi2);
+ /// Return wire coordinate closest to x.
-
- return Float_t(4.*fKx4*(TMath::ATan(ux2)-TMath::ATan(ux1))*
- fKy4*(TMath::ATan(uy2)-TMath::ATan(uy1)));
+ Int_t n = Int_t(x/Pitch());
+ Float_t wire = (x>0) ? n+0.5 : n-0.5;
+ return Pitch()*wire;
}
-Int_t AliMUONResponseV0::DigitResponse(Int_t digit, AliMUONTransientDigit* /*where*/)
+//______________________________________________________________________________
+void
+AliMUONResponseV0::DisIntegrate(const AliMUONHit& hit, TList& digits, Float_t timeDif)
{
- // 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);
- digit+=(Int_t)noise;
- if ( digit <= ZeroSuppression()) digit = 0;
- // if ( digit > MaxAdc()) digit=MaxAdc();
- if ( digit > Saturation()) digit=Saturation();
-
- return digit;
-}
-
-
-
+ /// 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();
+ Double_t hitX = hit.X() ;
+ Double_t hitY = hit.Y() ;
+ Double_t hitZ = hit.Z() ;
+ // Width of the integration area
+ Double_t dx = SigmaIntegration()*ChargeSpreadX();
+ Double_t dy = SigmaIntegration()*ChargeSpreadY();
+
+ //Modify to take the tailing effect.
+ if(fIsTailEffect){
+ Double_t locX,locY,locZ,globXCenter,globYCenter,globZ;
+ Int_t para = 5; // This parameter is a natural number(excluding zero), higher the value less is the tailing effect
+ Double_t termA = 1.0;
+ Double_t termB = 1.0;
+ if(para>0){
+ for ( Int_t cath = AliMp::kCath0; cath <= AliMp::kCath1; ++cath )
+ {
+ // Get an iterator to loop over pads, within the given area.
+ const AliMpVSegmentation* seg =
+ AliMpSegmentation::Instance()
+ ->GetMpSegmentation(detElemId,AliMp::GetCathodType(cath));
+ AliMp::PlaneType plane = seg->PlaneType();
+
+ if(plane == AliMp::kBendingPlane) {
+ Global2Local(detElemId,hitX,hitY,hitZ,locX,locY,locZ);
+ AliMpPad pad = seg->PadByPosition(locX,locY,kFALSE);
+ if(pad.IsValid()){
+ Double_t locYCenter = pad.GetPositionY();
+ Double_t locXCenter = pad.GetPositionX();
+ const AliMUONGeometryTransformer* transformer = muon()->GetGeometryTransformer();
+ transformer->Local2Global(detElemId,locXCenter,locYCenter,locZ,globXCenter,globYCenter,globZ);
+ for(Int_t itime = 0; itime<para; itime++)
+ termA *= 10.0;
+
+ for(Int_t itime = 0; itime<Int_t((2*para) + 1); itime++)
+ termB *= (hitY - globYCenter) ;
+
+ hitY = hitY + termA*termB;
+ }// if the pad is a valid one
+ }// if bending plane
+ }// cathode loop
+ }// if para > 0 condn
+ }// if tail effect
+ // Use that (dx,dy) to specify the area upon which
+ // we will iterate to spread charge into.
+ Double_t x,y,z;
+ Global2Local(detElemId,hitX,hitY,hitZ,x,y,z);
+ x = GetAnod(x);
+ AliMpArea area(x,y,dx,dy);
+
+ // Get pulse height from energy loss.
+ Float_t qtot = IntPH(hit.Eloss());
+
+ // If from a pileup event we apply a reduction factor to the charge
+ if (timeDif!=0){
+ qtot = AliMUONConstants::ReducedQTot(qtot,timeDif);
+ }
+
+ // Get the charge correlation between cathodes.
+ Float_t currentCorrel = TMath::Exp(gRandom->Gaus(0.0,ChargeCorrel()/2.0));
+ for ( Int_t cath = AliMp::kCath0; cath <= AliMp::kCath1; ++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 =
+ AliMpSegmentation::Instance()
+ ->GetMpSegmentation(detElemId,AliMp::GetCathodType(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.GetPositionX(),area.GetPositionY(),
+ kFALSE);
+ if ( pad.IsValid() )
+ {
+ AliDebug(1, Form("Got an invalid iterator bug (area.Position() is within "
+ " DE but the iterator is void) for detElemId %d cath %d",
+ detElemId,cath));
+ }
+ else
+ {
+ AliDebug(1, 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.GetPositionX(),area.GetPositionY(),
+ area.GetDimensionX(),area.GetDimensionY()));
+ }
+ 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(TVector2(x,y)-TVector2(pad.GetPositionX(),pad.GetPositionY())-
+ TVector2(pad.GetDimensionX(),pad.GetDimensionY()));
+ TVector2 upperRight(lowerLeft + TVector2(pad.GetDimensionX(),pad.GetDimensionY())*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(detElemId,pad.GetManuId(),
+ pad.GetManuChannel(),cath);
+ d->SetPadXY(pad.GetIx(),pad.GetIy());
+ d->SetCharge(icharge);
+ digits.Add(d);
+ }
+ it->Next();
+ }
+ delete it;
+ }
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff