]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - ITS/AliITSsimulationSSD.cxx
Adding some further mother volumes to speed-up the overlap checking and particle...
[u/mrichter/AliRoot.git] / ITS / AliITSsimulationSSD.cxx
index f6a6ee6aa6bef66c808c75ccd8111db60bac1119..48de1a28e2b6fc49fcce9eaae4f3f9d3f7eb9367 100644 (file)
@@ -19,6 +19,8 @@
 #include <stdlib.h>
 #include <Riostream.h>
 #include <TObjArray.h>
+#include <TRandom.h>
+
 #include "AliITSmodule.h"
 #include "AliITSMapA2.h"
 #include "AliITSpList.h"
 #include "AliITSgeom.h"
 #include "AliITSsimulationSSD.h"
 #include "AliITSTableSSD.h"
-//#include "AliITSresponseSSD.h"
 
 ClassImp(AliITSsimulationSSD)
 ////////////////////////////////////////////////////////////////////////
 //                                                                    //
 // Author: Enrico Fragiacomo                                          //
 //         enrico.fragiacomo@ts.infn.it                               //
-// Last revised: march 2006                                           // 
+// Last revised: june 2008                                            // 
 //                                                                    //
 // AliITSsimulationSSD is the simulation of SSD.                     //
 ////////////////////////////////////////////////////////////////////////
@@ -111,6 +112,7 @@ AliITSsimulationSSD& AliITSsimulationSSD::operator=(
   this->fDriftVel[1] = s.fDriftVel[1];
   return *this;
 }
+/*
 //______________________________________________________________________
 AliITSsimulation& AliITSsimulationSSD::operator=(
                                          const AliITSsimulation &s){
@@ -122,6 +124,7 @@ AliITSsimulation& AliITSsimulationSSD::operator=(
   
   return *this;
 }
+*/
 //______________________________________________________________________
 AliITSsimulationSSD::AliITSsimulationSSD(const AliITSsimulationSSD &source):
     AliITSsimulation(source),
@@ -218,33 +221,33 @@ void AliITSsimulationSSD::HitsToAnalogDigits(AliITSmodule *mod,
   
   AliITSsegmentationSSD* seg = (AliITSsegmentationSSD*)GetSegmentationModel(2);
   
-    TObjArray *hits = mod->GetHits();
-    Int_t nhits     = hits->GetEntriesFast();
-    if (nhits<=0) return;
-    AliITSTableSSD * tav = new AliITSTableSSD(GetNStrips());
-    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
-      // Returns also energy loss
-      if(GetDebug(4)){
-       cout << i << " ";
-       cout << mod->GetHit(i)->GetXL() << " "<<mod->GetHit(i)->GetYL();
-       cout << " " << mod->GetHit(i)->GetZL();
-       cout << endl;
-      } // end if
-      if (mod->LineSegmentL(i, x0, x1, y0, y1, z0, z1, de, idtrack)) {
-       HitToDigit(module, x0, y0, z0, x1, y1, z1, de,tav);
-       if (lasttrack != idtrack || i==(nhits-1)) {
-         GetList(idtrack,i,module,pList,tav);
-       } // end if
-       lasttrack=idtrack;
+  TObjArray *hits = mod->GetHits();
+  Int_t nhits     = hits->GetEntriesFast();
+  if (nhits<=0) return;
+  AliITSTableSSD * tav = new AliITSTableSSD(GetNStrips());
+  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
+    // Returns also energy loss
+    if(GetDebug(4)){
+      cout << i << " ";
+      cout << mod->GetHit(i)->GetXL() << " "<<mod->GetHit(i)->GetYL();
+      cout << " " << mod->GetHit(i)->GetZL();
+      cout << endl;
+    } // end if
+    if (mod->LineSegmentL(i, x0, x1, y0, y1, z0, z1, de, idtrack)) {
+      HitToDigit(module, x0, y0, z0, x1, y1, z1, de,tav);
+      if (lasttrack != idtrack || i==(nhits-1)) {
+       GetList(idtrack,i,module,pList,tav);
       } // end if
-    }  // end loop over hits
-    delete tav; tav=0;
-    return;
+      lasttrack=idtrack;
+    } // end if
+  }  // end loop over hits
+  delete tav; tav=0;
+  return;
 }
 //----------------------------------------------------------------------
 void AliITSsimulationSSD::HitToDigit(Int_t module, Double_t x0, Double_t y0, 
@@ -256,8 +259,8 @@ void AliITSsimulationSSD::HitToDigit(Int_t module, Double_t x0, Double_t y0,
   
   AliITSsegmentationSSD* seg = (AliITSsegmentationSSD*)GetSegmentationModel(2);
   // Turns hits in SSD module into one or more digits.
-  Float_t tang[2] = {0.0,0.0};
-  seg->Angles(tang[0], tang[1]);//stereo<<->tan(stereo)~=stereo
+  //Float_t tang[2] = {0.0,0.0};
+  //seg->Angles(tang[0], tang[1]);//stereo<<->tan(stereo)~=stereo
   Double_t x, y, z;
   Double_t dex=0.0, dey=0.0, dez=0.0; 
   Double_t pairs; // pair generation energy per step.
@@ -272,59 +275,118 @@ void AliITSsimulationSSD::HitToDigit(Int_t module, Double_t x0, Double_t y0,
   // Enery loss is equally distributed among steps
   de    = de/numOfSteps;
   pairs = de/GetIonizeE(); // e-h pairs generated
+
+  //-----------------------------------------------------
+  // stepping
+  //-----------------------------------------------------
   for(Int_t j=0; j<numOfSteps; j++) {     // stepping
+
     x = x0 + (j+0.5)*dex;
     y = y0 + (j+0.5)*dey;
     if ( y > (seg->Dy()/2+10)*1.0E-4 ) {
       // check if particle is within the detector
       Warning("HitToDigit",
-             "hit out of detector y0=%e,y=%e,dey=%e,j =%e module=%d",
-             y0,y,dey,j,module);
+             "hit out of detector y0=%e,y=%e,dey=%e,j =%d module=%d,  exceed=%e",
+             y0,y,dey,j,module, y-(seg->Dy()/2+10)*1.0E-4);
       return;
     } // end if
     z = z0 + (j+0.5)*dez;
+
     if(GetDebug(4)) cout <<"HitToDigit "<<x<<" "<<y<<" "<<z<< " "
                         <<dex<<" "<<dey<<" "<<dez<<endl;
+
+    if(seg->GetLayer()==6) {
+      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;
+    
     // calculate drift time
     // y is the minimum path
     tdrift[0] = (y+(seg->Dy()*1.0E-4)/2)/GetDriftVelocity(0);
+    
+    w = xp; // P side strip number
+    
+    if((w<(-0.5)) || (w>(GetNStrips()-0.5))) {
+      // this check rejects hits in regions not covered by strips
+      // 0.5 takes into account boundaries 
+      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
+    if(tdrift[k]<0) return;
+    
+    sigma[k] = TMath::Sqrt(2*GetDiffConst(k)*tdrift[k]);
+    sigma[k] /= (GetStripPitch()*1.0E-4);  //units of Pitch
+    
+    if(sigma[k]==0.0) {        
+      Error("HitToDigit"," sigma[%d]=0",k);
+      exit(0);
+    } // end if
+    
+    par0[k] = pairs;
+    // we integrate the diffusion gaussian from -3sigma to 3sigma 
+    inf[k] = w - 3*sigma[k]; // 3 sigma from the gaussian average  
+    sup[k] = w + 3*sigma[k]; // 3 sigma from the gaussian average
+    // IntegrateGaussian does the actual
+    // integration of diffusion gaussian
+    IntegrateGaussian(k, par0[k], w, sigma[k], inf[k], sup[k],tav);
+    
+    //------------------------------------------------------
+    // end Pside
+    //-------------------------------------------------------
+    
+    //------------------------------------------------------
+    // Nside
+    //-------------------------------------------------------
+    k=1;
     tdrift[1] = ((seg->Dy()*1.0E-4)/2-y)/GetDriftVelocity(1);
     
-    for(Int_t k=0; k<2; k++) {   // both sides    remember: 0=Pside 1=Nside
-      
-      tang[k]=TMath::Tan(tang[k]);
-      
-      // w is the coord. perpendicular to the strips
-      Float_t xp=x*1.e+4,zp=z*1.e+4; // microns
-      seg->GetPadTxz(xp,zp);
-      if(k==0) w = xp; // P side strip number
-      else w = zp; // N side strip number
-      
-      if((w<(-0.5)) || (w>(GetNStrips()-0.5))) {
-       // this check rejects hits in regions not covered by strips
-       // 0.5 takes into account boundaries 
-       if(GetDebug(4)) cout << "x,z="<<x<<","<<z<<" w="<<w
-                            <<" Nstrips="<<GetNStrips()<<endl;
-       return; // There are dead region on the SSD sensitive volume.
-      } // end if
-      
+    //tang[k]=TMath::Tan(tang[k]);
+    
+    w = zp; // N side strip number
+    
+    if((w<(-0.5)) || (w>(GetNStrips()-0.5))) {
+      // this check rejects hits in regions not covered by strips
+      // 0.5 takes into account boundaries 
+      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
-      if(tdrift[k]<0) return;
-      sigma[k] = TMath::Sqrt(2*GetDiffConst(k)*tdrift[k]);
-      sigma[k] /= (GetStripPitch()*1.0E-4);  //units of Pitch
-      if(sigma[k]==0.0) {      
-       Error("HitToDigit"," sigma[%d]=0",k);
-       exit(0);
-      } // end if
-      
-      par0[k] = pairs;
-      // we integrate the diffusion gaussian from -3sigma to 3sigma 
-      inf[k] = w - 3*sigma[k]; // 3 sigma from the gaussian average  
-      sup[k] = w + 3*sigma[k]; // 3 sigma from the gaussian average
-      // IntegrateGaussian does the actual
-      // integration of diffusion gaussian
-      IntegrateGaussian(k, par0[k], w, sigma[k], inf[k], sup[k],tav);
-    }  // end for loop over side (0=Pside, 1=Nside)      
+    if(tdrift[k]<0) return;
+    
+    sigma[k] = TMath::Sqrt(2*GetDiffConst(k)*tdrift[k]);
+    sigma[k] /= (GetStripPitch()*1.0E-4);  //units of Pitch
+    
+    if(sigma[k]==0.0) {        
+      Error("HitToDigit"," sigma[%d]=0",k);
+      exit(0);
+    } // end if
+    
+    par0[k] = pairs;
+    // we integrate the diffusion gaussian from -3sigma to 3sigma 
+    inf[k] = w - 3*sigma[k]; // 3 sigma from the gaussian average  
+    sup[k] = w + 3*sigma[k]; // 3 sigma from the gaussian average
+    // IntegrateGaussian does the actual
+    // integration of diffusion gaussian
+    IntegrateGaussian(k, par0[k], w, sigma[k], inf[k], sup[k],tav);
+    
+    //-------------------------------------------------
+    // end Nside
+    //-------------------------------------------------
+    
+    
   } // end stepping
 }
 
@@ -334,12 +396,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, 
@@ -349,7 +411,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
@@ -359,9 +421,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);
@@ -373,37 +435,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 
@@ -417,20 +480,14 @@ void AliITSsimulationSSD::ApplyCoupling(AliITSpList *pList,Int_t module) {
 void AliITSsimulationSSD::ApplyDeadChannels(Int_t module) {
   // Kill dead channels setting gain to zero
 
-  Int_t deadentries;
-
   AliITSCalibrationSSD* res = (AliITSCalibrationSSD*)GetCalibrationModel(module);
 
-  deadentries = res->GetDeadPChannelsList().GetSize();
-  //cout<<module<<" "<<deadentries<<endl;
-  for(Int_t i=0; i<deadentries; i++) {
-    res->AddGainP(res->GetDeadPChannelsList().At(i),0.0);
-  }
+  for(Int_t i=0;i<GetNStrips();i++){
 
-  deadentries = res->GetDeadNChannelsList().GetSize();
-  for(Int_t i=0; i<deadentries; i++) {
-    res->AddGainN(res->GetDeadNChannelsList().At(i),0.0);
-  }
+    if(res->IsPChannelBad(i)) res->SetGainP(i,0.0);
+    if(res->IsNChannelBad(i)) res->SetGainN(i,0.0);
+
+  } // loop over strips 
 
 }
 
@@ -583,6 +640,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
@@ -590,6 +648,7 @@ void AliITSsimulationSSD::ChargeToSignal(Int_t module,AliITSpList *pList) {
        // if strip is dead -> gain=0
        if( ((k==0)&&(res->GetGainP(ix)==0)) || ((k==1)&&(res->GetGainN(ix)==0))) continue;
        
+       signal = fMapA2->GetSignal(k,ix);
        // signal has to be uncalibrated
        // In real life, gains are supposed to be calculated from calibration runs,
        // stored in the calibration DB and used in the reconstruction
@@ -598,15 +657,18 @@ void AliITSsimulationSSD::ChargeToSignal(Int_t module,AliITSpList *pList) {
        else signal /= res->GetGainN(ix);
 
        // signal is converted in unit of ADC
-       signal = res->GetDEvToADC(fMapA2->GetSignal(k,ix));
+       signal = res->GetSSDDEvToADC(signal);
        if(signal>4096.) signal = 4096.;//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;
@@ -693,4 +755,3 @@ istream &operator>>(istream &os,AliITSsimulationSSD &source){
 
 
 
-