]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - ITS/AliITSsimulationSSD.cxx
add getter to cut on n cells
[u/mrichter/AliRoot.git] / ITS / AliITSsimulationSSD.cxx
index 4c01c9bdc717a1d37a2e7b3ec9e5557d0b80be3e..b28c28803109989f5f936ffa553b04f95ed25b7d 100644 (file)
@@ -21,6 +21,7 @@
 #include <TObjArray.h>
 #include <TRandom.h>
 
+#include <TGeoGlobalMagField.h>
 #include "AliITSmodule.h"
 #include "AliITSMapA2.h"
 #include "AliITSpList.h"
 #include "AliITShit.h"
 #include "AliITSdigitSSD.h"
 #include "AliRun.h"
+#include "AliMagF.h"
 #include "AliITSgeom.h"
 #include "AliITSsimulationSSD.h"
 #include "AliITSTableSSD.h"
-//#include "AliITSresponseSSD.h"
+#include <TF1.h>
+#include "AliMathBase.h"
 
+using std::endl;
+using std::cout;
 ClassImp(AliITSsimulationSSD)
 ////////////////////////////////////////////////////////////////////////
 //                                                                    //
@@ -52,7 +57,12 @@ AliITSsimulationSSD::AliITSsimulationSSD():AliITSsimulation(),
 fMapA2(0),
 fIonE(0.0),
 fDifConst(),
-fDriftVel(){
+fDriftVel(),
+fTimeResponse(NULL),
+fLorentz(kFALSE),
+fTanLorAngP(0),
+fTanLorAngN(0)
+{
     //default Constructor
     //Inputs:
     // none.
@@ -68,7 +78,12 @@ AliITSsimulation(dettyp),
 fMapA2(0),
 fIonE(0.0),
 fDifConst(),
-fDriftVel(){
+fDriftVel(),
+fTimeResponse(NULL),
+fLorentz(kFALSE),
+fTanLorAngP(0),
+fTanLorAngN(0)
+{
     // Constructor 
     // Input:
     //   AliITSDetTypeSim    Pointer to the SSD dettype to be used
@@ -77,6 +92,7 @@ fDriftVel(){
     // Return
     //   A standard constructed AliITSsimulationSSD class
 
+  fTimeResponse = new TF1("ftimeresponse",".5*x*exp(1.-.5*x)");
     Init();
 }
 //----------------------------------------------------------------------
@@ -90,13 +106,38 @@ void AliITSsimulationSSD::Init(){
   // Return
   //   none.
   AliITSsegmentationSSD* seg = (AliITSsegmentationSSD*)GetSegmentationModel(2);
+  AliITSSimuParam* simpar = fDetType->GetSimuParam();
   
   SetDriftVelocity(); // use default values in .h file
   SetIonizeE();       // use default values in .h file
   SetDiffConst();     // use default values in .h file
   fpList           = new AliITSpList(2,GetNStrips());
   fMapA2           = new AliITSMapA2(seg);
+  SetLorentzDrift(simpar->GetSSDLorentzDrift());
+  if (fLorentz) SetTanLorAngle();
 }
+
+//______________________________________________________________________
+Bool_t AliITSsimulationSSD::SetTanLorAngle() {
+    // This function set the Tangent of the Lorentz angles. 
+    // output: Bool_t : kTRUE in case of success
+    //
+
+    if(!fDetType) {
+      AliError("AliITSsimulationSPD::SetTanLorAngle: AliITSDetTypeSim* fDetType not set ");
+      return kFALSE;}
+
+    AliITSSimuParam* simpar = fDetType->GetSimuParam();
+    AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
+    if (!fld) AliFatal("The field is not initialized");
+    Double_t bz = fld->SolenoidField();
+
+    fTanLorAngN = TMath::Tan( simpar->LorentzAngleElectron(bz) );
+    fTanLorAngP = TMath::Tan( simpar->LorentzAngleHole(bz) );
+
+    return kTRUE;
+}
+
 //______________________________________________________________________
 AliITSsimulationSSD& AliITSsimulationSSD::operator=(
                                          const AliITSsimulationSSD &s){
@@ -111,6 +152,10 @@ AliITSsimulationSSD& AliITSsimulationSSD::operator=(
   this->fDifConst[1] = s.fDifConst[1];
   this->fDriftVel[0] = s.fDriftVel[0];
   this->fDriftVel[1] = s.fDriftVel[1];
+  this->fTimeResponse = s.fTimeResponse;
+  this->fLorentz   = s.fLorentz;
+  this->fTanLorAngP = s.fTanLorAngP;
+  this->fTanLorAngN = s.fTanLorAngN;
   return *this;
 }
 /*
@@ -132,7 +177,12 @@ AliITSsimulationSSD::AliITSsimulationSSD(const AliITSsimulationSSD &source):
 fMapA2(source.fMapA2),
 fIonE(source.fIonE),
 fDifConst(),
-fDriftVel(){
+fDriftVel(),
+fTimeResponse(source.fTimeResponse),
+fLorentz(source.fLorentz),
+fTanLorAngP(source.fTanLorAngP),
+fTanLorAngN(source.fTanLorAngN)
+{
   // copy constructor
   fDifConst[0] = source.fDifConst[0];
   fDifConst[1] = source.fDifConst[1];
@@ -143,6 +193,7 @@ fDriftVel(){
 AliITSsimulationSSD::~AliITSsimulationSSD() {
   // destructor
   delete fMapA2;
+  delete fTimeResponse;
   //delete fDCS;
 }
 //______________________________________________________________________
@@ -219,6 +270,8 @@ void AliITSsimulationSSD::HitsToAnalogDigits(AliITSmodule *mod,
   Double_t x1=0.0, y1=0.0, z1=0.0;
   Double_t de=0.0;
   Int_t module = mod->GetIndex();
+  Double_t tof = 0.;
+  
   
   AliITSsegmentationSSD* seg = (AliITSsegmentationSSD*)GetSegmentationModel(2);
   
@@ -229,6 +282,7 @@ void AliITSsimulationSSD::HitsToAnalogDigits(AliITSmodule *mod,
   module = mod->GetIndex();
   if ( mod->GetLayer() == 6 ) seg->SetLayer(6);
   if ( mod->GetLayer() == 5 ) seg->SetLayer(5);
+
   for(Int_t i=0; i<nhits; i++) {    
     // LineSegmentL returns 0 if the hit is entering
     // If hits is exiting returns positions of entering and exiting hits
@@ -240,6 +294,16 @@ void AliITSsimulationSSD::HitsToAnalogDigits(AliITSmodule *mod,
       cout << endl;
     } // end if
     if (mod->LineSegmentL(i, x0, x1, y0, y1, z0, z1, de, idtrack)) {
+
+      // Scale down dE/dx according to the hit's TOF wrt to the trigger
+      // Necessary for pileup simulation
+      // EF - 21/04/09
+      tof = mod->GetHit(i)->GetTOF();
+      tof *= 1.E+6; // convert time in microsecond
+      if(tof<2.) de = de * fTimeResponse->Eval(-1.*tof+2.);
+      else de = 0.;
+      //
+
       HitToDigit(module, x0, y0, z0, x1, y1, z1, de,tav);
       if (lasttrack != idtrack || i==(nhits-1)) {
        GetList(idtrack,i,module,pList,tav);
@@ -269,7 +333,15 @@ void AliITSsimulationSSD::HitToDigit(Int_t module, Double_t x0, Double_t y0,
   Double_t tdrift[2] = {0.,0.}; // time of drift
   Double_t w;
   Double_t inf[2], sup[2], par0[2];                 
-  
+  // Set up corrections for Lorentz drift (ExB)
+  Double_t tanLorAngP = fTanLorAngP;
+  Double_t tanLorAngN = fTanLorAngN;
+  if(seg->GetLayer()==6) {
+    tanLorAngP = -1.*fTanLorAngP;
+    tanLorAngN = -1.*fTanLorAngN;
+  }
+
   // Steps in the module are determined "manually" (i.e. No Geant)
   // NumOfSteps divide path between entering and exiting hits in steps 
   Int_t numOfSteps = NumOfSteps(x1, y1, z1, dex, dey, dez);
@@ -300,16 +372,23 @@ void AliITSsimulationSSD::HitToDigit(Int_t module, Double_t x0, Double_t y0,
       y=-y; // Lay6 module has sensor up-side-down!!!
     }
     
-    // w is the coord. perpendicular to the strips
-    //    Float_t xp=x*1.e+4,zp=z*1.e+4; // microns    
-    Float_t xp=x,zp=z; 
-    seg->GetPadTxz(xp,zp);
-
     Int_t k;
     //---------------------------------------------------------
     // Pside
     //------------------------------------------------------------
     k=0;
+
+    // w is the coord. perpendicular to the strips
+    //    Float_t xp=x*1.e+4,zp=z*1.e+4; // microns    
+    Float_t xp=x,zp=z; 
+
+    // correction for the Lorentz's angle
+    if(fLorentz) {
+      Float_t deltaxp = (y+(seg->Dy()*1.0E-4)/2)*tanLorAngP;
+      xp+=deltaxp;  
+    }
+
+    seg->GetPadTxz(xp,zp);
     
     // calculate drift time
     // y is the minimum path
@@ -323,8 +402,7 @@ void AliITSsimulationSSD::HitToDigit(Int_t module, Double_t x0, Double_t y0,
       if(GetDebug(4)) cout << "Dead SSD region, x,z="<<x<<","<<z<<endl;
       return; // There are dead region on the SSD sensitive volume!!!
     } // end if
-    
-      // sigma is the standard deviation of the diffusion gaussian
+    // sigma is the standard deviation of the diffusion gaussian
     if(tdrift[k]<0) return;
     
     sigma[k] = TMath::Sqrt(2*GetDiffConst(k)*tdrift[k]);
@@ -351,6 +429,18 @@ void AliITSsimulationSSD::HitToDigit(Int_t module, Double_t x0, Double_t y0,
     // Nside
     //-------------------------------------------------------
     k=1;
+
+    xp=x; zp=z; 
+
+    // correction for the Lorentz's angle
+    if(fLorentz) {
+      Float_t deltaxn = ((seg->Dy()*1.0E-4)/2-y)*tanLorAngN;
+      xp+=deltaxn;
+    }
+    
+
+    seg->GetPadTxz(xp,zp);
+
     tdrift[1] = ((seg->Dy()*1.0E-4)/2-y)/GetDriftVelocity(1);
     
     //tang[k]=TMath::Tan(tang[k]);
@@ -397,12 +487,12 @@ void AliITSsimulationSSD::ApplyNoise(AliITSpList *pList,Int_t module){
   Int_t ix;
   Double_t signal,noise;
   AliITSCalibrationSSD* res =(AliITSCalibrationSSD*)GetCalibrationModel(module);
-  
+   
   // Pside
   for(ix=0;ix<GetNStrips();ix++){      // loop over strips
     
     // noise is gaussian
-    noise  = (Double_t) gRandom->Gaus(0,res->GetNoiseP().At(ix));
+    noise  = (Double_t) gRandom->Gaus(0,res->GetNoiseP(ix));
     
     // need to calibrate noise 
     // NOTE. noise from the calibration database comes uncalibrated, 
@@ -412,7 +502,7 @@ void AliITSsimulationSSD::ApplyNoise(AliITSpList *pList,Int_t module){
     
     // noise comes in ADC channels from the calibration database
     // It needs to be converted back to electronVolts
-    noise /= res->GetDEvToADC(1.);
+    noise /= res->GetSSDDEvToADC(1.);
     
     // Finally, noise is added to the signal
     signal = noise + fMapA2->GetSignal(0,ix);//get signal from map
@@ -422,9 +512,9 @@ void AliITSsimulationSSD::ApplyNoise(AliITSpList *pList,Int_t module){
   
     // Nside
   for(ix=0;ix<GetNStrips();ix++){      // loop over strips
-    noise  = (Double_t) gRandom->Gaus(0,res->GetNoiseN().At(ix));// give noise to signal
+    noise  = (Double_t) gRandom->Gaus(0,res->GetNoiseN(ix));// give noise to signal
     noise *= (Double_t) res->GetGainN(ix); 
-    noise /= res->GetDEvToADC(1.);
+    noise /= res->GetSSDDEvToADC(1.);
     signal = noise + fMapA2->GetSignal(1,ix);//get signal from map
     fMapA2->SetHit(1,ix,signal); // give back signal to map
     if(signal>0.0) pList->AddNoise(1,ix,module,noise);
@@ -436,37 +526,38 @@ void AliITSsimulationSSD::ApplyCoupling(AliITSpList *pList,Int_t module) {
   // Apply the effect of electronic coupling between channels
   Int_t ix;
   Double_t signal=0;
-  AliITSCalibrationSSD* res =(AliITSCalibrationSSD*)GetCalibrationModel(module);
-  
+  //AliITSCalibrationSSD* res =(AliITSCalibrationSSD*)GetCalibrationModel(module);
+  AliITSSimuParam* res = fDetType->GetSimuParam();
+    
   Double_t *contrLeft  = new Double_t[GetNStrips()];
   Double_t *contrRight = new Double_t[GetNStrips()];
   
   // P side coupling
   for(ix=0;ix<GetNStrips();ix++){
-    if(ix>0) contrLeft[ix] = fMapA2->GetSignal(0,ix-1)*res->GetCouplingPL();
+    if(ix>0) contrLeft[ix] = fMapA2->GetSignal(0,ix-1)*res->GetSSDCouplingPL();
     else contrLeft[ix] = 0.0;
-    if(ix<(GetNStrips()-1)) contrRight[ix] = fMapA2->GetSignal(0,ix+1)*res->GetCouplingPR();
+    if(ix<(GetNStrips()-1)) contrRight[ix] = fMapA2->GetSignal(0,ix+1)*res->GetSSDCouplingPR();
     else contrRight[ix] = 0.0;
   } // loop over strips 
   
   for(ix=0;ix<GetNStrips();ix++){
-    signal = contrLeft[ix] + contrRight[ix] - res->GetCouplingPL() * fMapA2->GetSignal(0,ix)
-      - res->GetCouplingPR() * fMapA2->GetSignal(0,ix);
+    signal = contrLeft[ix] + contrRight[ix] - res->GetSSDCouplingPL() * fMapA2->GetSignal(0,ix)
+      - res->GetSSDCouplingPR() * fMapA2->GetSignal(0,ix);
     fMapA2->AddSignal(0,ix,signal);
     if(signal>0.0) pList->AddNoise(0,ix,module,signal);
   } // loop over strips 
   
   // N side coupling
   for(ix=0;ix<GetNStrips();ix++){
-    if(ix>0) contrLeft[ix] = fMapA2->GetSignal(1,ix-1)*res->GetCouplingNL();
+    if(ix>0) contrLeft[ix] = fMapA2->GetSignal(1,ix-1)*res->GetSSDCouplingNL();
     else contrLeft[ix] = 0.0;
-    if(ix<(GetNStrips()-1)) contrRight[ix] = fMapA2->GetSignal(1,ix+1)*res->GetCouplingNR();
+    if(ix<(GetNStrips()-1)) contrRight[ix] = fMapA2->GetSignal(1,ix+1)*res->GetSSDCouplingNR();
     else contrRight[ix] = 0.0;
   } // loop over strips 
   
   for(ix=0;ix<GetNStrips();ix++){
-    signal = contrLeft[ix] + contrRight[ix] - res->GetCouplingNL() * fMapA2->GetSignal(0,ix)
-      - res->GetCouplingNR() * fMapA2->GetSignal(0,ix);
+    signal = contrLeft[ix] + contrRight[ix] - res->GetSSDCouplingNL() * fMapA2->GetSignal(0,ix)
+      - res->GetSSDCouplingNR() * fMapA2->GetSignal(0,ix);
     fMapA2->AddSignal(1,ix,signal);
     if(signal>0.0) pList->AddNoise(1,ix,module,signal);
   } // loop over strips 
@@ -484,8 +575,8 @@ void AliITSsimulationSSD::ApplyDeadChannels(Int_t module) {
 
   for(Int_t i=0;i<GetNStrips();i++){
 
-    if(res->IsPChannelBad(i)) res->AddGainP(i,0.0);
-    if(res->IsNChannelBad(i)) res->AddGainN(i,0.0);
+    if(res->IsPChannelBad(i)) res->SetGainP(i,0.0);
+    if(res->IsNChannelBad(i)) res->SetGainN(i,0.0);
 
   } // loop over strips 
 
@@ -498,7 +589,7 @@ Float_t AliITSsimulationSSD::F(Float_t av, Float_t x, Float_t s) {
     Float_t sigm2 = sqrt2*s;
     Float_t integral;
 
-    integral = 0.5 * TMath::Erf( (x - av) / sigm2);
+    integral = 0.5 * AliMathBase::ErfFast( (x - av) / sigm2);
     return integral;
 }
 //______________________________________________________________________
@@ -628,7 +719,7 @@ void AliITSsimulationSSD::GetList(Int_t label,Int_t hit,Int_t mod,
     tav->Clear();
 }
 //----------------------------------------------------------------------
-void AliITSsimulationSSD::ChargeToSignal(Int_t module,AliITSpList *pList) {
+void AliITSsimulationSSD::ChargeToSignal(Int_t module,const AliITSpList *pList) {
     // charge to signal
     static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
     Float_t threshold = 0.;
@@ -640,6 +731,7 @@ void AliITSsimulationSSD::ChargeToSignal(Int_t module,AliITSpList *pList) {
     Float_t charges[3] = {0.0,0.0,0.0};
     Float_t signal;
     AliITSCalibrationSSD* res =(AliITSCalibrationSSD*)GetCalibrationModel(module);
+    AliITSSimuParam* simpar = fDetType->GetSimuParam();
 
     for(Int_t k=0;k<2;k++){         // both sides (0=Pside, 1=Nside)
       for(Int_t ix=0;ix<GetNStrips();ix++){     // loop over strips
@@ -656,15 +748,18 @@ void AliITSsimulationSSD::ChargeToSignal(Int_t module,AliITSpList *pList) {
        else signal /= res->GetGainN(ix);
 
        // signal is converted in unit of ADC
-       signal = res->GetDEvToADC(signal);
-       if(signal>4096.) signal = 4096.;//if exceeding, accumulate last one
+       signal = res->GetSSDDEvToADC(signal);
+       if(signal>4095.) signal = 4095.;//if exceeding, accumulate last one
 
        // threshold for zero suppression is set on the basis of the noise
        // A good value is 3*sigma_noise
-       if(k==0) threshold = res->GetNoiseP().At(ix);
-       else threshold = res->GetNoiseN().At(ix);
-       threshold *= res->GetZSThreshold(); // threshold at 3 sigma noise
+       if(k==0) threshold = res->GetNoiseP(ix);
+       else threshold = res->GetNoiseN(ix);
+
+       threshold *= simpar->GetSSDZSThreshold(); // threshold at 3 sigma noise
+
        if(signal < threshold) continue;
+       //cout<<signal<<" "<<threshold<<endl;
 
        digits[0] = k;
        digits[1] = ix;
@@ -750,4 +845,3 @@ istream &operator>>(istream &os,AliITSsimulationSSD &source){
 
 
 
-