BesTofDigitizerEcV4_dbs.cc

Go to the documentation of this file.

//---------------------------------------------------------------------------//
//      BOOST --- BESIII Object_Oriented Simulation Tool                     //
//---------------------------------------------------------------------------//
//Description: This is the Digitizer for the MRPC with doubled sided readout
//Author: Matthias Ullrich
//Created: May, 2012
//Modified:
//Comment:
//---------------------------------------------------------------------------//
//$Id: BesTofDigitizerEcV4.cc
 
#include "BesTofDigitizerEcV4_dbs.hh"
#include "BesTofDigi.hh"
#include "BesTofHit.hh"
#include "G4DigiManager.hh"
#include "BesTofGeoParameter.hh"
#include "Randomize.hh"
#include "TMath.h"
#include <math.h>
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "G4Svc/IG4Svc.h"
#include "G4Svc/G4Svc.h"
#include "TH1F.h"
#include "TNtuple.h"
#include "TFile.h"
 
 
BesTofDigitizerEcV4_dbs::BesTofDigitizerEcV4_dbs()
{
 
 
  m_storedata=false;
 
  if(m_storedata)
    {
  FileName = "output_mrpc_digitizerecv4.root";
  m_TreeFile=new TFile(TString(FileName),"RECREATE","TreeFile");
  m_ecv4tree = new TTree("ecv4tree", "");   //Create Tree
  m_ecv4tree->Branch("m_partId",&m_partID,"m_partId/D");
  m_ecv4tree->Branch("m_stripidentifier",&m_stripidentifier,"m_stripidentifier/D");
  m_ecv4tree->Branch("m_trackIndex",&m_trackIndex,"m_trackIndex/D");
  m_ecv4tree->Branch("m_signal_pc",&m_signal_pc,"m_signal_pc/D");
  m_ecv4tree->Branch("m_time_threshold",&m_time_threshold,"m_time_threshold/D");
  m_ecv4tree->Branch("m_time_1sthit",&m_time_1sthit,"m_time_1sthit/D");
  m_ecv4tree->Branch("m_time_1",&m_time_1,"m_time_1/D");
  m_ecv4tree->Branch("m_time_2",&m_time_2,"m_time_2/D");
  m_ecv4tree->Branch("m_trans_time_1",&m_trans_time_1,"m_trans_time_1/D");
  m_ecv4tree->Branch("m_trans_time_2",&m_trans_time_2,"m_trans_time_2/D");
  m_ecv4tree->Branch("m_firedstrip",&m_firedstrip,"m_firedstrip/D");
  m_ecv4tree->Branch("m_numberions",&m_numberions,"m_numberions/D");
  m_ecv4tree->Branch("m_numberofgaps_with_avalanches",&m_numberofgaps_with_avalanches,"m_numberofgaps_with_avalanches/D");
  m_ecv4tree->Branch("m_numberofstripsfired",&m_numberofstripsfired,"m_numberofstripsfired/D");
  m_ecv4tree->Branch("m_multihit",&m_multihit,"m_multihit/D");
  m_ecv4tree->Branch("m_inefficiency",&m_inefficiency,"m_inefficiency/D");
  m_ecv4tree->Branch("m_particle_true",&m_particle_true,"m_particle_true/D");
  m_ecv4tree->Branch("m_particle_has_signal",&m_particle_has_signal,"m_particle_has_signal/D");
 
    }
 
}
 
 
BesTofDigitizerEcV4_dbs::~BesTofDigitizerEcV4_dbs()
{
 
  if(m_storedata)
    {
  m_TreeFile->Write();
  m_TreeFile->Close();
    }
   
}
 
void BesTofDigitizerEcV4_dbs::Digitize(ScintSingle* single_module, BesTofDigitsCollection* DC)
{
 
     trackIndex=0;
     partId=-999;
     module_mrpc=-999;
     time=-999.;
 
 
    m_besTofDigitsCollection = DC;
    G4DigiManager* digiManager = G4DigiManager::GetDMpointer();
    G4int THCID = digiManager->GetHitsCollectionID("BesTofHitsCollection");
    m_THC = (BesTofHitsCollection*) (digiManager->GetHitsCollection(THCID));
    if (m_THC)
    {
 
 
 
      G4int n_hit=single_module->GetHitIndexes_mrpc()->size(); 
      
 
      vector<adc_info> avalanche_info;//Will store the basic information for signal creation etc.                                                         
          avalanche_info.clear();
 
      partId = single_module->GetPartId();
      module_mrpc = single_module->GetModule_mrpc();
         
      
         
      //Process the hits
      int deadarea=0;//DeadArea Counter: How often can't we associate an readoutstrip?
      
      for(G4int i=0;i<n_hit;i++)
        {  
 
 
          BesTofHit* hit = (*m_THC)[( *(single_module->GetHitIndexes_mrpc()) )[i]];
          
                 
          G4ThreeVector help = hit->GetPos();
          G4int readoutstripnumber = Calculate_Readoutstrip_number(help.x()/mm,help.y()/mm,partId,module_mrpc);
          
          if( readoutstripnumber==0 )   {deadarea++; continue;}
 
          
          adc_info struct_adc_info;
          struct_adc_info.gap =0; struct_adc_info.readoutstrip =0; struct_adc_info.ions =0; struct_adc_info.zdistance =0.; struct_adc_info.trackindex=0;struct_adc_info.time=0.; 
          struct_adc_info.pdg_code=0;//Intialize all values
 
          struct_adc_info.gap=Get_gapnumber(partId, help.z());
          struct_adc_info.readoutstrip=readoutstripnumber;
          struct_adc_info.ions=hit->GetIons();
          struct_adc_info.zdistance=Distance_for_avalanche(struct_adc_info.gap,help.z(), partId);
          struct_adc_info.time = hit->GetTime();
          struct_adc_info.trackindex=hit->GetG4Index();
          struct_adc_info.pdg_code=hit->GetPDGcode();
          struct_adc_info.transtime_1=Calculate_strip_transition_time_1(help.x()/mm,help.y()/mm,partId,module_mrpc);
          struct_adc_info.transtime_2=Calculate_strip_transition_time_2(help.x()/mm,help.y()/mm,partId,module_mrpc);
              struct_adc_info.x=help.x()/mm;
              struct_adc_info.y=help.y()/mm;
              struct_adc_info.partId=partId;
              struct_adc_info.module=module_mrpc;
 
          if(struct_adc_info.ions>0 && struct_adc_info.readoutstrip!=0 && struct_adc_info.time!=0.) avalanche_info.push_back(struct_adc_info);
 
        }//Close for-loop over hits
      
 
 
      G4int avalanche_info_size = avalanche_info.size();
 
 
      //Determine which readoutstrips has been fired and push the depeding on the strip number into one vector
      vector <int> fired_readoutstrips;
 
      fired_readoutstrips.clear();
      for(G4int i=0; i<avalanche_info_size; i++)
        {
 
          bool dlh=true;
          for(G4int j=0;j<fired_readoutstrips.size();j++)
        {
          if(fired_readoutstrips[j]==avalanche_info[i].readoutstrip){dlh=false; break;}
        }
          if(dlh==true)
        {
                  fired_readoutstrips.push_back(avalanche_info[i].readoutstrip);
        }
        }
 
 
 
      //Determine the adcinfo for each fired strips with the smallest time -> These information will be stored within the digis                         
      vector <int> mybest_adcinfo_inreadoutstrips(fired_readoutstrips.size());
      for(G4int i=0;i< fired_readoutstrips.size();i++)
        {
          double time_help=10000.;
          for(G4int j=0; j<avalanche_info_size; j++)
        {
          if((time_help>avalanche_info[j].time) && avalanche_info[j].readoutstrip==fired_readoutstrips[i])
            {
              time_help = avalanche_info[j].time;
              mybest_adcinfo_inreadoutstrips[i]=j;
            }//close if
        }//close for j
        }//close for i
 
 
 
 
      //Calculate the signal strength and determine the time information and store these information into the digi collection                           
      double stripidentifier_1=0; //Each strip will be read out by 2 NINO chips, each on one side!
      double stripidentifier_2=0;
      
      for(G4int i=0; i<fired_readoutstrips.size();i++)
        {
 
 
 
          //Determine the total number of ions!                                                                                                           
          G4int sumions=0;
              for(G4int j=0; j<avalanche_info_size; j++)
        {
          if(avalanche_info[j].readoutstrip==fired_readoutstrips[i]) 
            {
            sumions = sumions+ avalanche_info[j].ions;
            }
 
        }
 
 
              //Determine the total number of gaps with avalanches
              G4int numberofgaps_with_avalanches=0;
              for(G4int j=1;j<=12;j++)
                {
                  for(G4int k=0; k<avalanche_info_size; k++) { if(avalanche_info[k].gap ==j && avalanche_info[k].ions >0){numberofgaps_with_avalanches++; break;} }
                }
              
 
 
 
          G4double signal_pc=0; G4double time_of_threshold_exceedance=0;
          simulate_avalanche(fired_readoutstrips[i], &avalanche_info,&time_of_threshold_exceedance,&signal_pc);
          signal_pc=signal_pc;//1e-12*A*s;                                                                                                             
          time_of_threshold_exceedance=time_of_threshold_exceedance*s;//Stores the time when the threshold is reached.
 
 
          if(signal_pc<0.015) continue; //We only store signals wich are over the threshold
 
          
          trackIndex = avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].trackindex;
          
 
 
 
          // We add the threshold crossing time for the avalanche. This will be corrected later within the reconstruction.
          double resulting_phi = Calculate_resulting_phi(avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].x,avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].y
                      ,avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].partId,avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].module);
 
          if(resulting_phi<0)//The ninos one the left side of a module (odd  numbers) will see the signal earlier!
                         //Time 1 is the smaller transition time.
        {
          time = avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].time + time_of_threshold_exceedance + avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].transtime_1; 
          time_digi2=avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].time + time_of_threshold_exceedance + avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].transtime_2; 
        }
          else//The ninos one the right side of a module (even numbers) will see the signal earlier!
          //Time 1 is the smaller transition time.
        {
          time = avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].time + time_of_threshold_exceedance + avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].transtime_2; 
          time_digi2=avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].time + time_of_threshold_exceedance + avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].transtime_1; 
        }
 
          
          //time smear #Correct 
          //Uncertainities which are not considered in the digitization process. Following values are assumed:
          //Avalanche Process : Already considered with time of threshold!
          //FEE (NINO)        : Depending on the charge, will be calculated in the next step! = NINO_RES
          //Hit position      : 10 ps 
          //Readout           : 27 ps
          //Electronics/cable : 20 ps
          //-----------------------------------------
          // Sum =sqrt(10**2 + 27**2 +20**2 NINO_RES**2 ) time smear
 
          double  femtocolomb= 1000.;
 
          double NINO_RES=0;
          NINO_RES = Get_NINO_leadingedge_resolution(signal_pc*femtocolomb);
          double smear_time = sqrt(NINO_RES*NINO_RES + 27.0*27.0 + 10.0*10.0 + 20.0*20.0); //This is in ps
 
          time = Smear_gaussian(time,0,smear_time*0.001*ns);
          time_digi2=Smear_gaussian(time_digi2,0,smear_time*0.001*ns);
                
 
          stripidentifier_1= Produce_unique_identifier(module_mrpc,(avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].readoutstrip)*2-1);
          stripidentifier_2= Produce_unique_identifier(module_mrpc,(avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].readoutstrip)*2);
          
 
          //Resolution of Pulse length:
          // 27 ps TDC
          // Timeresolution of pulse length
          double NINO_pulselength_resolution = Get_NINO_pulselength_resolution( signal_pc*femtocolomb );
          double smear_time_2 = sqrt(NINO_pulselength_resolution*NINO_pulselength_resolution + 27.0*27.0 );
     
          double time2 =chargetotime(signal_pc*femtocolomb)+time; //The function returns the time in ns 
          time2 = Smear_gaussian(time2,0,smear_time_2*0.001*ns);
 
 
          //so far I do assume,that both Nino chips see the same signals!
          double time2_digi2= chargetotime(signal_pc*femtocolomb)+time_digi2;
          time2_digi2 = Smear_gaussian(time2_digi2,0,smear_time_2*0.001*ns);
 
 
          if(m_storedata)
        {
          m_partID=partId;
          m_stripidentifier=stripidentifier_1;
          m_trackIndex=trackIndex;
          m_signal_pc=signal_pc;
          m_time_threshold= time_of_threshold_exceedance;
          m_time_1sthit=avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].time;
          m_time_1=time;
          m_time_2=time2;
          m_trans_time_1=avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].transtime_1;
          m_trans_time_2=avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].transtime_2;
          m_firedstrip=fired_readoutstrips[i];
          m_numberions=sumions;
              m_numberofgaps_with_avalanches=numberofgaps_with_avalanches;
              m_numberofstripsfired = fired_readoutstrips.size();
          m_x=avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].x;
              m_y=avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].y;
          m_ecv4tree->Fill();
        }
 
          
 
          
 
          BesTofDigi *digi = new BesTofDigi;
          digi->SetTrackIndex(trackIndex);
          digi->SetPartId(partId);
          digi->SetModule_mrpc(stripidentifier_1);
          digi->SetTime1_mrpc(time);
          digi->SetTime2_mrpc(time2);
          m_besTofDigitsCollection->insert(digi);
 
 
 
 
          BesTofDigi *digi2 = new BesTofDigi;
          digi2->SetTrackIndex(trackIndex);
          digi2->SetPartId(partId);
          digi2->SetModule_mrpc(stripidentifier_2);
          digi2->SetTime1_mrpc(time_digi2);
          digi2->SetTime2_mrpc(time2_digi2);
          m_besTofDigitsCollection->insert(digi2);
 
          
 
        }//close i
 
 
 
       
    }//close if(m_THC)
 
} //Digitize()
 
 
 
 
 
//////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////
//////
//////  Function used within the Digitization process ////
//////
/////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////
 
 
 
G4double BesTofDigitizerEcV4_dbs::Smear_gaussian(G4double variable, G4double gauss_mean,G4double gauss_sigma)
{
  
  G4double gaussrand = G4RandGauss::shoot(gauss_mean,gauss_sigma);//returns x-value
  return (variable+gaussrand);
  
}
 
 
 
 
G4double BesTofDigitizerEcV4_dbs::Calculate_resulting_phi(G4double x_mm,G4double y_mm, G4int partId_f, G4int module_mrpc_f)
{
 
 
 
 
  G4double phi= atan2(y_mm,x_mm)*180./pi;
 
 
  if(-180<phi && phi<90){phi=phi+270;} // Enter a coordinate system, which has on its highest point the 0 degree, and increases to the left/right side! The increasing depends on the endcap
  else {phi=phi-90;}
 
 
  G4double result_phi=0;
  if(partId_f!=3 && partId_f!=6) result_phi=0;
  else result_phi=9.97312;  //This is the angle to which the first and second layer are shifted in degree!!!
 
  if(partId_f==3 || partId_f==4)  result_phi = -((19-module_mrpc_f-1)*20 - phi + result_phi);
  else result_phi = (module_mrpc_f-1)*20 - phi + result_phi;  // If > 0  Readoutstripnumber -> 2,4,6,8,10,12,14.16,18,20,22,24
 
  //The 18th and the 1st module are located between 350 and 10 degree. This has to be taken into account
  if(partId_f==4 && module_mrpc_f==18 && result_phi>350) result_phi=result_phi-360;
  if(partId_f==5 && module_mrpc_f==1 && result_phi<-350) result_phi=result_phi+360;
 
  return result_phi;
  // If > 0  Readoutstripnumber -> 2,4,6,8,10,12,14.16,18,20,22,24 
 
}
 
 
 
//This function caluclates depending on the energy deposition which strips cause a signal!!
G4int BesTofDigitizerEcV4_dbs::Calculate_Readoutstrip_number(G4double x_mm,G4double y_mm, G4int partId_f, G4int module_mrpc_f)
 
{
       
//Within this functions all numbers are in mm or in degree!!!!
     
  G4int readoutstripnumber=0;
       
  G4double result_radius= sqrt(x_mm*x_mm+y_mm*y_mm);
  G4double phi= atan2(y_mm,x_mm)*180./pi; 
       
       
  if(-180<phi && phi<90){phi=phi+270;} // Enter a coordinate system, which has on its highest point the 0 degree, and increases to the left/right side! The increasing depends on the endcap
  else {phi=phi-90;}
       
       
  G4double result_phi=0;
  if(partId_f!=3 && partId_f!=6) result_phi=0;
  else result_phi=9.97312;  //This is the angle to which the first and second layer are shifted in degree!!!
 
 
 
  if(partId_f==3 || partId_f==4)  result_phi = -((19-module_mrpc_f-1)*20 - phi + result_phi);
  else result_phi = (module_mrpc_f-1)*20 - phi + result_phi;  // If > 0  Readoutstripnumber -> 2,4,6,8,10,12,14.16,18,20,22,24
  
  //The 18th and the 1st module are located between 350 and 10 degree. This has to be taken into account
  if(partId_f==4 && module_mrpc_f==18 && result_phi>350) result_phi=result_phi-360;
  if(partId_f==5 && module_mrpc_f==1 && result_phi<-350) result_phi=result_phi+360;
  
 
 
  result_radius= result_radius*cos(result_phi*pi/180.0) - 634.56 + 159.5 - 3.82;    //634.56 is the origin of coordinates from the module, 159.5 is the distance from the transformation of the readout strips...
                                                           // 3.82 is the ks_trafo     
 
  double result_x=fabs((result_radius/cos(result_phi*pi/180.0)+634.56-159.5+3.82)*sin(result_phi*pi/180)); //Calculate the resulting position x!
 
              
 
  if( -12.5 <result_radius && result_radius < 12.5 && result_x <= 44 ) readoutstripnumber=1;
  if( (-12.5 +(25.0+4.0)*1 )  <result_radius && result_radius < (12.5+(25.0+4.0)*1) && result_x <=46) readoutstripnumber=2;
  if( (-12.5 +(25.0+4.0)*2 )  <result_radius && result_radius < (12.5+(25.0+4.0)*2) && result_x <=49) readoutstripnumber=3;
  if( (-12.5 +(25.0+4.0)*3 )  <result_radius && result_radius < (12.5+(25.0+4.0)*3) && result_x <=51) readoutstripnumber=4;
  if( (-12.5 +(25.0+4.0)*4 )  <result_radius && result_radius < (12.5+(25.0+4.0)*4) && result_x <=54) readoutstripnumber=5;
  if( (-12.5 +(25.0+4.0)*5 )  <result_radius && result_radius < (12.5+(25.0+4.0)*5) && result_x <=56) readoutstripnumber=6;
  if( (-12.5 +(25.0+4.0)*6 )  <result_radius && result_radius < (12.5+(25.0+4.0)*6) && result_x <=59) readoutstripnumber=7;
  if( (-12.5 +(25.0+4.0)*7 )  <result_radius && result_radius < (12.5+(25.0+4.0)*7) && result_x <=61) readoutstripnumber=8;
  if( (-12.5 +(25.0+4.0)*8 )  <result_radius && result_radius < (12.5+(25.0+4.0)*8) && result_x <=64) readoutstripnumber=9;
  if( (-12.5 +(25.0+4.0)*9 )  <result_radius && result_radius < (12.5+(25.0+4.0)*9) && result_x <=66) readoutstripnumber=10;
  if( (-12.5 +(25.0+4.0)*10 ) <result_radius && result_radius < (12.5+(25.0+4.0)*10) && result_x <=69) readoutstripnumber=11;
  if( (-12.5 +(25.0+4.0)*11 ) <result_radius && result_radius < (12.5+(25.0+4.0)*11) && result_x <=71) readoutstripnumber=12;
       
                 
  return readoutstripnumber;
                       
}
 
 
 
 
//This function is almost the same as above, but it does not consider the "dead" readout area! This is important for the reconstruction!
//So it is possible to determine the readoutstrip, even if no strip is hit and one can do the extrapolation for the MDC Tracks
G4int BesTofDigitizerEcV4_dbs::Calculate_Readoutstrip_number_continuum(G4double x_mm,G4double y_mm, G4int partId_f, G4int module_mrpc_f)
{
 
  //Within this functions all numbers are in mm or in degree!!!!                                                                                                                                                
 
  G4int readoutstripnumber=0;
 
  G4double result_radius= sqrt(x_mm*x_mm+y_mm*y_mm);
  G4double phi= atan2(y_mm,x_mm)*180./pi;
 
 
  if(-180<phi && phi<90){phi=phi+270;}//See funtion above!
  else {phi=phi-90;}
 
 
  G4double result_phi=0;
  if(partId_f!=3 && partId_f!=6) result_phi=0;
  else result_phi=9.97312;  //This is the angle to which the first and second layer are shifted in degree!!!                                                                                                    
  if(partId_f==3 || partId_f==4)  result_phi = -((19-module_mrpc_f-1)*20 - phi + result_phi);
  else result_phi = (module_mrpc_f-1)*20 - phi + result_phi;  // If > 0  Readoutstripnumber -> 2,4,6,8,10,12,14.16,18,20,22,24                                                                                       
 
  //The 18th and the 1st module are located between 350 and 10 degree. This has to be taken into account
  if(partId_f==4 && module_mrpc_f==18 && result_phi>350) result_phi=result_phi-360;
  if(partId_f==5 && module_mrpc_f==1 && result_phi<-350) result_phi=result_phi+360;
  
 
 
  result_radius= result_radius*cos(result_phi*pi/180.0) - 634.56 + 159.5-3.82;    //634.56 is the origin of coordinates from the module, 159.5 is the distance from the transformation of the readout strips...                               
 
  //The +-12.0 are chosen arbitrarily. In this region we assume reconstructed tracks as MRPC tracks!
  if( (-12.5 -12.0) <=result_radius && result_radius < 12.5+2.0 ) readoutstripnumber=1;
  if( (-12.5 +(25.0+4.0)*1-2.0 )  <=result_radius && result_radius < (12.5+(25.0+4.0)*1)+2.0 ) readoutstripnumber=2;
  if( (-12.5 +(25.0+4.0)*2-2.0 )  <=result_radius && result_radius < (12.5+(25.0+4.0)*2) +2.0) readoutstripnumber=3;
  if( (-12.5 +(25.0+4.0)*3-2.0 )  <=result_radius && result_radius < (12.5+(25.0+4.0)*3) +2.0) readoutstripnumber=4;
  if( (-12.5 +(25.0+4.0)*4-2.0 )  <=result_radius && result_radius < (12.5+(25.0+4.0)*4) +2.0) readoutstripnumber=5;
  if( (-12.5 +(25.0+4.0)*5-2.0 )  <=result_radius && result_radius < (12.5+(25.0+4.0)*5) +2.0) readoutstripnumber=6;
  if( (-12.5 +(25.0+4.0)*6-2.0 )  <=result_radius && result_radius < (12.5+(25.0+4.0)*6) +2.0) readoutstripnumber=7;
  if( (-12.5 +(25.0+4.0)*7-2.0 )  <=result_radius && result_radius < (12.5+(25.0+4.0)*7) +2.0) readoutstripnumber=8;
  if( (-12.5 +(25.0+4.0)*8-2.0 )  <=result_radius && result_radius < (12.5+(25.0+4.0)*8) +2.0) readoutstripnumber=9;
  if( (-12.5 +(25.0+4.0)*9-2.0 )  <=result_radius && result_radius < (12.5+(25.0+4.0)*9) +2.0) readoutstripnumber=10;
  if( (-12.5 +(25.0+4.0)*10-2.0 ) <=result_radius && result_radius < (12.5+(25.0+4.0)*10) +2.0) readoutstripnumber=11;
  if( (-12.5 +(25.0+4.0)*11-2.0 ) <=result_radius && result_radius < (12.5+(25.0+4.0)*11) +12.0) readoutstripnumber=12;
  
 
  return readoutstripnumber;
 
}
 
 
 
 
 
 
 
//Calculates the smaller one of the transition times
G4double BesTofDigitizerEcV4_dbs::Calculate_strip_transition_time_1(G4double x_mm,G4double y_mm, G4int partId_f, G4int module_mrpc_f)
{
  G4double signal_velocity_copper_mrpc_strip= 0.8*0.299792458;  // unit:        percentage * c_0 * mm/ps   
 
  //Within this functions all numbers are in mm or in degree!!!!
  G4double length_to_border_of_the_readoutstrip=0;
  G4double result_radius= sqrt(x_mm*x_mm+y_mm*y_mm); //Radius where the particle did hit the detector
  G4double phi= atan2(y_mm,x_mm)*180./pi;            //Phi angle of the particle hitting the detector
 
  if(-180<phi && phi<90){phi=phi+270;} // Transformation into a coordinate-system, whose coordinates do increase to the right/left side. Left or right does depend on the endcap!
  else {phi=phi-90;}
 
  G4double result_phi=0;
  if(partId_f!=3 && partId_f!=6) result_phi=0;
  else result_phi=9.97312;  //This is the angle to which the first and second layer are shifted in degree!!!
 
 
  if(partId_f==3 || partId_f==4)  result_phi = -((19-module_mrpc_f-1)*20 - phi + result_phi);
  else result_phi = (module_mrpc_f-1)*20 - phi + result_phi;  // If > 0  Readoutstripnumber -> 2,4,6,8,10,12,14.16,18,20,22,24
  //The 18th and the 1st module are located between 350 and 10 degree. This has to be taken into account
 
  if(partId_f==4 && module_mrpc_f==18 && result_phi>350) result_phi=result_phi-360;
  if(partId_f==5 && module_mrpc_f==1 && result_phi<-350) result_phi=result_phi+360;
 
  //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  //The Phi is now the is the phi angle for each module and not for the whole assembly 
  //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 
  result_radius= result_radius*cos(result_phi*pi/180.0) - 634.56 + 159.5 - 3.82;    //634.56 is the origin of coordinates from the module, 159.5 is the distance from the transformation of the readout strips... 
  // 3.82 is the ks_trafo
 
  //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  //The radius is now the the radius within one module. The point of origin for the
  //radius is the middle of the lowest readout cell
  //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 
  double result_x=fabs((result_radius/cos(result_phi*pi/180.0)+634.56-159.5+3.82)*sin(result_phi*pi/180)); //Calculate the resulting position x!
 
  
  if( -12.5 <result_radius && result_radius < 12.5 && result_x <= 44 )
    {
      length_to_border_of_the_readoutstrip=sqrt( (44-result_x)*(44-result_x) +  fabs((result_radius-0.)*(result_radius-0.)));
    }
  if( (-12.5 +(25.0+4.0)*1 )  <result_radius && result_radius < (12.5+(25.0+4.0)*1) && result_x <=46)
    {
      length_to_border_of_the_readoutstrip=sqrt( (46-result_x)*(46-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*1)));
    }
  if( (-12.5 +(25.0+4.0)*2 )  <result_radius && result_radius < (12.5+(25.0+4.0)*2) && result_x <=49)
    {
      length_to_border_of_the_readoutstrip=sqrt( (49-result_x)*(49-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*2)));
    }
  if( (-12.5 +(25.0+4.0)*3 )  <result_radius && result_radius < (12.5+(25.0+4.0)*3) && result_x <=51)
    {
      length_to_border_of_the_readoutstrip=sqrt( (51-result_x)*(51-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*3)));
    }
  if( (-12.5 +(25.0+4.0)*4 )  <result_radius && result_radius < (12.5+(25.0+4.0)*4) && result_x <=54)
    {
      length_to_border_of_the_readoutstrip=sqrt( (54-result_x)*(54-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*4)));
    }
  if( (-12.5 +(25.0+4.0)*5 )  <result_radius && result_radius < (12.5+(25.0+4.0)*5) && result_x <=56)
    {
      length_to_border_of_the_readoutstrip=sqrt( (56-result_x)*(56-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*5)));
    }
  if( (-12.5 +(25.0+4.0)*6 )  <result_radius && result_radius < (12.5+(25.0+4.0)*6) && result_x <=59)
    {
      length_to_border_of_the_readoutstrip=sqrt( (59-result_x)*(59-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*6)));
    }
  if( (-12.5 +(25.0+4.0)*7 )  <result_radius && result_radius < (12.5+(25.0+4.0)*7) && result_x <=61)
    {
      length_to_border_of_the_readoutstrip=sqrt( (61-result_x)*(61-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*7)));
    }
  if( (-12.5 +(25.0+4.0)*8 )  <result_radius && result_radius < (12.5+(25.0+4.0)*8) && result_x <=64)
    {
      length_to_border_of_the_readoutstrip=sqrt( (64-result_x)*(64-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*8)));
    }
  if( (-12.5 +(25.0+4.0)*9 )  <result_radius && result_radius < (12.5+(25.0+4.0)*9) && result_x <=66)
    {
      length_to_border_of_the_readoutstrip=sqrt( (66-result_x)*(66-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*9)));
    }
  if( (-12.5 +(25.0+4.0)*10 ) <result_radius && result_radius < (12.5+(25.0+4.0)*10) && result_x <=69)
    {
      length_to_border_of_the_readoutstrip=sqrt( (69-result_x)*(69-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*10)));
    }
  if( (-12.5 +(25.0+4.0)*11 ) <result_radius && result_radius < (12.5+(25.0+4.0)*11) && result_x <=71)
    {
      length_to_border_of_the_readoutstrip=sqrt( (71-result_x)*(71-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*11)));
    }
  double transition_time = length_to_border_of_the_readoutstrip/signal_velocity_copper_mrpc_strip/1000.;  //length is in mm , 3*10^8 m/s = 0.3 mm/ps  --> Unit of transitionstime is ns (due to division by 1000.)
  return transition_time;
 
}//close calculate_strip_transition_time_1
 
 
 
//Calculates the larger one of the transition times
G4double BesTofDigitizerEcV4_dbs::Calculate_strip_transition_time_2(G4double x_mm,G4double y_mm, G4int partId_f, G4int module_mrpc_f)
{
  G4double signal_velocity_copper_mrpc_strip= 0.8*0.299792458;  // unit:        percentage * c_0 * mm/ps   
 
  //Within this functions all numbers are in mm or in degree!!!!
  G4double length_to_border_of_the_readoutstrip=0;
  G4double result_radius= sqrt(x_mm*x_mm+y_mm*y_mm); //Radius where the particle did hit the detector
  G4double phi= atan2(y_mm,x_mm)*180./pi;            //Phi angle of the particle hitting the detector
 
  if(-180<phi && phi<90){phi=phi+270;} // Transformation into a coordinate-system, whose coordinates do increase to the right/left side. Left or right does depend on the endcap!
  else {phi=phi-90;}
 
  G4double result_phi=0;
  if(partId_f!=3 && partId_f!=6) result_phi=0;
  else result_phi=9.97312;  //This is the angle to which the first and second layer are shifted in degree!!!
 
 
  if(partId_f==3 || partId_f==4)  result_phi = -((19-module_mrpc_f-1)*20 - phi + result_phi);
  else result_phi = (module_mrpc_f-1)*20 - phi + result_phi;  // If > 0  Readoutstripnumber -> 2,4,6,8,10,12,14.16,18,20,22,24
  //The 18th and the 1st module are located between 350 and 10 degree. This has to be taken into account
 
  if(partId_f==4 && module_mrpc_f==18 && result_phi>350) result_phi=result_phi-360;
  if(partId_f==5 && module_mrpc_f==1 && result_phi<-350) result_phi=result_phi+360;
 
  //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  //The Phi is now the is the phi angle for each module and not for the whole assembly 
  //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 
  result_radius= result_radius*cos(result_phi*pi/180.0) - 634.56 + 159.5 - 3.82;    //634.56 is the origin of coordinates from the module, 159.5 is the distance from the transformation of the readout strips... 
  // 3.82 is the ks_trafo
 
  //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  //The radius is now the the radius within one module. The point of origin for the
  //radius is the middle of the lowest readout cell
  //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 
  double result_x=fabs((result_radius/cos(result_phi*pi/180.0)+634.56-159.5+3.82)*sin(result_phi*pi/180)); //Calculate the resulting position x!
 
  
  if( -12.5 <result_radius && result_radius < 12.5 && result_x <= 44 )
    {
      length_to_border_of_the_readoutstrip=sqrt( (48+result_x)*(48+result_x) +  fabs((result_radius-0.)*(result_radius-0.)));
    }
  if( (-12.5 +(25.0+4.0)*1 )  <result_radius && result_radius < (12.5+(25.0+4.0)*1) && result_x <=46)
    {
      length_to_border_of_the_readoutstrip=sqrt( (50+result_x)*(50+result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*1)));
    }
  if( (-12.5 +(25.0+4.0)*2 )  <result_radius && result_radius < (12.5+(25.0+4.0)*2) && result_x <=49)
    {
      length_to_border_of_the_readoutstrip=sqrt( (53+result_x)*(53+result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*2)));
    }
  if( (-12.5 +(25.0+4.0)*3 )  <result_radius && result_radius < (12.5+(25.0+4.0)*3) && result_x <=51)
    {
      length_to_border_of_the_readoutstrip=sqrt( (55+result_x)*(55+result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*3)));
    }
  if( (-12.5 +(25.0+4.0)*4 )  <result_radius && result_radius < (12.5+(25.0+4.0)*4) && result_x <=54)
    {
      length_to_border_of_the_readoutstrip=sqrt( (58+result_x)*(58+result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*4)));
    }
  if( (-12.5 +(25.0+4.0)*5 )  <result_radius && result_radius < (12.5+(25.0+4.0)*5) && result_x <=56)
    {
      length_to_border_of_the_readoutstrip=sqrt( (60+result_x)*(60+result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*5)));
    }
  if( (-12.5 +(25.0+4.0)*6 )  <result_radius && result_radius < (12.5+(25.0+4.0)*6) && result_x <=59)
    {
      length_to_border_of_the_readoutstrip=sqrt( (63+result_x)*(63+result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*6)));
    }
  if( (-12.5 +(25.0+4.0)*7 )  <result_radius && result_radius < (12.5+(25.0+4.0)*7) && result_x <=61)
    {
      length_to_border_of_the_readoutstrip=sqrt( (65+result_x)*(65+result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*7)));
    }
  if( (-12.5 +(25.0+4.0)*8 )  <result_radius && result_radius < (12.5+(25.0+4.0)*8) && result_x <=64)
    {
      length_to_border_of_the_readoutstrip=sqrt( (68+result_x)*(68+result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*8)));
    }
  if( (-12.5 +(25.0+4.0)*9 )  <result_radius && result_radius < (12.5+(25.0+4.0)*9) && result_x <=66)
    {
      length_to_border_of_the_readoutstrip=sqrt( (70+result_x)*(70+result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*9)));
    }
  if( (-12.5 +(25.0+4.0)*10 ) <result_radius && result_radius < (12.5+(25.0+4.0)*10) && result_x <=69)
    {
      length_to_border_of_the_readoutstrip=sqrt( (73+result_x)*(73+result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*10)));
    }
  if( (-12.5 +(25.0+4.0)*11 ) <result_radius && result_radius < (12.5+(25.0+4.0)*11) && result_x <=71)
    {
      length_to_border_of_the_readoutstrip=sqrt( (75+result_x)*(75+result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*11)));
    }
 
  double transition_time = length_to_border_of_the_readoutstrip/signal_velocity_copper_mrpc_strip/1000.;  //length is in mm , 3*10^8 m/s = 0.3 mm/ps  --> Unit of transitionstime is ns (due to division by 1000.)
  return transition_time;
 
}//close calculate_strip_transition_time_2
 
 
 
//////////////////////////////////////////////////////////////////////
// The following function for the transitiontime are used within    //
// the reconstruction! I did only put them here, that when          //
// the signal velocity is changed, it is not forgotten for those    //
// functions!                                                       //
////////////////////////////////////////////////////////////////////// 
 
 
 
double BesTofDigitizerEcV4_dbs::GetTransitionTime_dbsmatch(double delta_t,int my_strip)
{
  double signal_velocity_copper_mrpc_strip= 0.8*0.299792458;   // unit:        percentage * c_0 * mm/ps 
  double strip_length=0;
  
  if(my_strip==1 || my_strip==2 ) strip_length=88;
  else if(my_strip==3 || my_strip==4 ) strip_length=92;
  else if(my_strip==5 || my_strip==6 ) strip_length=98; 
  else if(my_strip==7 || my_strip==8 ) strip_length=102;
  else if(my_strip==9 || my_strip==10 ) strip_length=108;
  else if(my_strip==11 || my_strip==12 ) strip_length=112;
  else if(my_strip==13 || my_strip==14 ) strip_length=118;
  else if(my_strip==15 || my_strip==16 ) strip_length=122;
  else if(my_strip==17 || my_strip==18 ) strip_length=128;
  else if(my_strip==19|| my_strip==20 ) strip_length=132;
  else if(my_strip==21|| my_strip==22 ) strip_length=138;
  else if(my_strip==23 || my_strip==24 ) strip_length=142;
  
 
 
  double trans_time =(strip_length/signal_velocity_copper_mrpc_strip/1000+delta_t)/2.;
    
  
  return trans_time;
 
}
 
 
 
 
 
 
 
double BesTofDigitizerEcV4_dbs::GetTransitionTime_extrap_track(G4double x_mm, G4double y_mm,  int partId_f , int module_mrpc_f , int  my_strip)
{
 
  G4double signal_velocity_copper_mrpc_strip= 0.8*0.299792458;  // unit:        percentage * c_0 * mm/ps 
 
  //Within this functions all numbers are in mm or in degree!!!!
  G4double length_to_border_of_the_readoutstrip=0;
  G4double result_radius= sqrt(x_mm*x_mm+y_mm*y_mm); //Radius where the particle did hit the detector
  G4double phi= atan2(y_mm,x_mm)*180./pi;            //Phi angle of the particle hitting the detector
 
  if(-180<phi && phi<90){phi=phi+270;} // Transformation into a coordinate-system, whose coordinates do increase to the right/left side. Left or right does depend on the endcap!
  else {phi=phi-90;}
 
  G4double result_phi=0;
  if(partId_f!=3 && partId_f!=6) result_phi=0.;
  else result_phi=9.97312;  //This is the angle to which the first and second layer are shifted in degree!!!
 
  if(partId_f==3 || partId_f==4)  result_phi = -((19-module_mrpc_f-1)*20. - phi + result_phi);
  else result_phi = (module_mrpc_f-1)*20. - phi + result_phi;  // If > 0  Readoutstripnumber -> 2,4,6,8,10,12,14.16,18,20,22,24
  //The 18th and the 1st module are located between 350 and 10 degree. This has to be taken into account
 
  if(partId_f==4 && module_mrpc_f==18 && result_phi>350) result_phi=result_phi-360;
  if(partId_f==5 && module_mrpc_f==1 && result_phi<-350) result_phi=result_phi+360;
 
 
  result_radius= result_radius*cos(result_phi*pi/180.0) - 634.56 + 159.5 - 3.82;  
 
  double result_x=((result_radius/cos(result_phi*pi/180.0)+634.56-159.5+3.82)*sin(result_phi*pi/180)); 
 
 
 
  if( (my_strip%2)==0 )
   {
     if(my_strip==2)
       length_to_border_of_the_readoutstrip=sqrt( (44-result_x)*(44-result_x) +  fabs((result_radius-0.)*(result_radius-0.)));
     else if(my_strip==4)
       length_to_border_of_the_readoutstrip=sqrt( (46-result_x)*(46-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*1)));
     else if(my_strip==6)
       length_to_border_of_the_readoutstrip=sqrt( (49-result_x)*(49-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*2)));
     else  if(my_strip==8)
       length_to_border_of_the_readoutstrip=sqrt( (51-result_x)*(51-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*3)));
     else if(my_strip==10)
       length_to_border_of_the_readoutstrip=sqrt( (54-result_x)*(54-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*4)));
     else if(my_strip==12)
       length_to_border_of_the_readoutstrip=sqrt( (56-result_x)*(56-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*5)));
     else if(my_strip==14)
       length_to_border_of_the_readoutstrip=sqrt( (59-result_x)*(59-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*6)));
     else if(my_strip==16)
       length_to_border_of_the_readoutstrip=sqrt( (61-result_x)*(61-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*7)));
     else if(my_strip==18)
       length_to_border_of_the_readoutstrip=sqrt( (64-result_x)*(64-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*8)));
     else if(my_strip==20)
       length_to_border_of_the_readoutstrip=sqrt( (66-result_x)*(66-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*9)));
     else if(my_strip==22)
       length_to_border_of_the_readoutstrip=sqrt( (69-result_x)*(69-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*10)));
     else if(my_strip==24)
       length_to_border_of_the_readoutstrip=sqrt( (71-result_x)*(71-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*11)));
   }
  else if((my_strip%2)!=0)
   {
     if(my_strip==1)
       length_to_border_of_the_readoutstrip=sqrt( (-44-result_x)*(-44-result_x) +  fabs((result_radius-0.)*(result_radius-0.)));
     else if(my_strip==3)
       length_to_border_of_the_readoutstrip=sqrt( (-46-result_x)*(-46-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*1)));
     else if(my_strip==5)
       length_to_border_of_the_readoutstrip=sqrt( (-49-result_x)*(-49-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*2)));
     else if(my_strip==7)
       length_to_border_of_the_readoutstrip=sqrt( (-51-result_x)*(-51-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*3)));
     else if(my_strip==9)
       length_to_border_of_the_readoutstrip=sqrt( (-54-result_x)*(-54-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*4)));
     else if(my_strip==11)
       length_to_border_of_the_readoutstrip=sqrt( (-56-result_x)*(-56-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*5)));
     else if(my_strip==13)
       length_to_border_of_the_readoutstrip=sqrt( (-59-result_x)*(-59-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*6)));
     else if(my_strip==15)
       length_to_border_of_the_readoutstrip=sqrt( (-61-result_x)*(-61-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*7)));
     else if(my_strip==17)
       length_to_border_of_the_readoutstrip=sqrt( (-64-result_x)*(-64-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*8)));
     else if(my_strip==19)
       length_to_border_of_the_readoutstrip=sqrt( (-66-result_x)*(-66-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*9)));
     else if(my_strip==21)
       length_to_border_of_the_readoutstrip=sqrt( (-69-result_x)*(-69-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*10)));
     else if(my_strip==23)
       length_to_border_of_the_readoutstrip=sqrt( (-71-result_x)*(-71-result_x) +  fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*11)));
   }
 
 double transition_time = length_to_border_of_the_readoutstrip/signal_velocity_copper_mrpc_strip/1000.;  //length is in mm , 3*10^8 m/s = 0.3 mm/ps  --> Unit of transitionstime is ns (due to division by 1000.)
 return transition_time;
 
 
 
}
 
 
 
 
 
G4double BesTofDigitizerEcV4_dbs::Average_transition_time_dbs(G4int my_module)
{
  G4double signal_velocity_copper_mrpc_strip= 0.8*0.299792458;  // unit:        percentage * c_0 * mm/ps
                                                                                                                                                                                    
    G4double length_to_border_of_the_readoutstrip=0.;
 
    if(my_module==1 || my_module==2) length_to_border_of_the_readoutstrip=44;
    if(my_module==3 || my_module==4) length_to_border_of_the_readoutstrip=46;
    if(my_module==5 || my_module==6) length_to_border_of_the_readoutstrip=49;
    if(my_module==7 || my_module==8) length_to_border_of_the_readoutstrip=51;
    if(my_module==9 || my_module==10) length_to_border_of_the_readoutstrip=54;
    if(my_module==11 || my_module==12) length_to_border_of_the_readoutstrip=56;
    if(my_module==13 || my_module==14) length_to_border_of_the_readoutstrip=59;
    if(my_module==15 || my_module==16) length_to_border_of_the_readoutstrip=61;
    if(my_module==17 || my_module==18) length_to_border_of_the_readoutstrip=64;
    if(my_module==19 || my_module==20) length_to_border_of_the_readoutstrip=66;
    if(my_module==21 || my_module==22) length_to_border_of_the_readoutstrip=69;
    if(my_module==23 || my_module==24) length_to_border_of_the_readoutstrip=71;
 
    double transition_time = length_to_border_of_the_readoutstrip/signal_velocity_copper_mrpc_strip/1000.;
    return transition_time;
}
 
 
 
 
 
G4int BesTofDigitizerEcV4_dbs::Produce_unique_identifier(G4int module_mrpc_f, G4int readoutstripnumber_f)
{
  G4int returnvalue=0;
 
 
  returnvalue= module_mrpc_f*25 + readoutstripnumber_f;
  
  return returnvalue;  
}
 
 
 
 
G4int BesTofDigitizerEcV4_dbs::Get_module_mrpc_from_unique_identifier(G4int unique_identifier_f)
{
 
  G4int returnvalue=0;
 
  returnvalue= unique_identifier_f/25;
 
 
  return returnvalue;
}
 
 
G4int BesTofDigitizerEcV4_dbs::Get_stripnumber_from_unique_identifier(G4int unique_identifier_f)
{
 
  G4int returnvalue=0;
 
  returnvalue= unique_identifier_f%25;
 
  return returnvalue;
}
 
 
 
 
G4int BesTofDigitizerEcV4_dbs::Get_gapnumber(G4int partID, G4double z_position)
{
 
 
 
  G4double pos_west_1=-1343.5*mm;  // = partID 5                                                                                                              
  G4double pos_west_2=-1368.5*mm;  // = partID 6                                                                                                              
  G4double pos_east_1=+1343.5*mm; // = partID 4                                                                                                               
  G4double pos_east_2=+1368.5*mm; // = partID 3                                                                                                               
 
 
  //We transform our z-coordinate into the system of one MRPC module. PoV= Interactionpoint:                                                                  
  //The gaps 1-6 will have negative values, the gaps 7-12 will have positive values!                                                                          
  G4double z_zero=0;
  if(partID==3)
    {z_zero= z_position-pos_east_2;}
  else if(partID==4)
    {z_zero= z_position-pos_east_1;}
  else if(partID==5)
    {z_zero = pos_west_1 - z_position ;}//z_position is negative in this case!!!!                                                                             
  else if(partID==6)
    {z_zero = pos_west_2 - z_position ;}//z_position is negative in this case!!!!                                                                             
  else
    {G4cout << "TofDigitzer:    Error in part ID" <<G4endl; return 0;}
 
  //std::cout << "z_pos " << z_position << std::endl;                                                                                                         
  //std::cout << "z_zero " << std::setprecision(100) <<z_zero << std::endl;                                                                                   
 
 
  //There is a porblem with the numbers. The correct ones are the rounded ones. However in this case too often 0 is returned!                   
  if((-4.5999*mm) >= z_zero && z_zero>(-4.8201*mm) ) {return 1;}
  else if((-3.9799*mm) >= z_zero && z_zero>(-4.2001*mm) ){return 2;}
  else if((-3.3599*mm) >= z_zero && z_zero>(-3.5801*mm) ){return 3;}
  else if((-2.7399*mm) >= z_zero && z_zero>(-2.9601*mm) ){return 4;}
  else if((-2.1199*mm) >= z_zero && z_zero>(-2.3401*mm) ){return 5;}
  else if((-1.4999*mm) >= z_zero && z_zero>(-1.7201*mm)) {return 6;}
  else if((1.4999*mm) <= z_zero && z_zero<=(1.7201*mm) ) {return 7;}
  else if((2.1199*mm) <= z_zero && z_zero<=(2.3401*mm) ) {return 8;}
  else if((2.7399*mm) <= z_zero && z_zero<=(2.9601*mm) ) {return 9;}
  else if((3.3599*mm) <= z_zero && z_zero<=(3.5801*mm) ) {return 10;}
  else if((3.9799*mm) <= z_zero && z_zero<=(4.2001*mm) ) {return 11;}
  else if((4.5999*mm)  <= z_zero && z_zero<=(4.8201*mm) ){return 12;}
  else {G4cout << "TofDigitzer:  Error in Hit-Position!  " <<G4endl; return 0; }
 
}
 
 
 
G4double BesTofDigitizerEcV4_dbs::Distance_for_avalanche(G4int gap, G4double z_position, G4int partID)
{
 
  //All the values are in mm                                                                                                                                  
  G4double pos_west_1=-1343.5*mm;  // = partID 5                                                                                                              
  G4double pos_west_2=-1368.5*mm;  // = partID 6                                                                                                              
  G4double pos_east_1=+1343.5*mm; // = partID 4                                                                                                               
  G4double pos_east_2=+1368.5*mm; // = partID 3                                                                                                               
 
 
  //We transform our z-coordinate into the system of one MRPC module. PoV= Interactionpoint:                                                                  
  //The gaps 1-6 will have negative values, the gaps 7-12 will have positive values!                                                                          
  G4double z_zero=0;
  if(partID==3)
    {z_zero= z_position-pos_east_2;}
  else if(partID==4)
    {z_zero= z_position-pos_east_1;}
  else if(partID==5)
    {z_zero = pos_west_1 - z_position ;}//z_position is negative in this case!!!!                                                                             
  else if(partID==6)
    {z_zero = pos_west_2 - z_position ;}//z_position is negative in this case!!!!                                                                             
  else
    {G4cout << "TofDigitzer:    Error in part ID" <<G4endl; return 0;}
 
  //cout << "z_zero  "<< setprecision(100)<< z_zero << G4endl;                                                                                                
 
  /*
  if(gap==1) G4cout << (-z_zero-4.60*mm) << "    "<< gap << G4endl ;
  if(gap==2) G4cout  <<(-z_zero-3.98*mm)<< "    "<< gap << G4endl;
  if(gap==3) G4cout<< (-z_zero-3.36*mm)<< "    "<< gap << G4endl;
  if(gap==4) G4cout<< (-z_zero-2.74*mm)<< "    "<< gap << G4endl;
  if(gap==5) G4cout<< (-z_zero-2.12*mm)<< "    "<< gap << G4endl;
  if(gap==6) G4cout<< (-z_zero-1.5*mm)<< "    "<< gap << G4endl;
  if(gap==7) G4cout<< (1.72-z_zero*mm)<< "    "<< gap << G4endl;
  if(gap==8) G4cout<< (2.34-z_zero*mm)<< "    "<< gap << G4endl;
  if(gap==9) G4cout<< (2.96-z_zero*mm)<< "    "<< gap << G4endl;
  if(gap==10)G4cout<< (3.58-z_zero*mm)<< "    "<< gap << G4endl;
  if(gap==11)G4cout<< (4.20-z_zero*mm)<< "    "<< gap << G4endl;
  if(gap==12)G4cout<< (4.82-z_zero*mm)<< "    "<< gap << G4endl;
  */
 
 
 
 
 
  if(gap==1) return (-z_zero-4.60*mm);
  if(gap==2) return (-z_zero-3.98*mm);
  if(gap==3) return (-z_zero-3.36*mm);
  if(gap==4) return (-z_zero-2.74*mm);
  if(gap==5) return (-z_zero-2.12*mm);
  if(gap==6) return (-z_zero-1.5*mm);
  if(gap==7) return (1.72-z_zero*mm);
  if(gap==8) return (2.34-z_zero*mm);
  if(gap==9) return (2.96-z_zero*mm);
  if(gap==10)return (3.58-z_zero*mm);
  if(gap==11)return (4.20-z_zero*mm);
  if(gap==12)return (4.82-z_zero*mm);
 
 
}
 
//////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////
//////
//////  Function used for avalanche simulation       //// 
//////
/////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////
 
 
 
void BesTofDigitizerEcV4_dbs::simulate_avalanche(G4int actual_strip, vector<adc_info> * avalanche_info, double * avalanche_threshold_time, double *induced_charge)
{
 
  //Originally this function was a stand alone program. That's why the style of programming may be a little strange:
  //It depends on a method described in the Dissertation of Christian Lippmann and is known in there as 1-D Model.
  //Simulation with this model have been done for a chamber of the same type as implemented: Simulation study on the operation of a multi-gap resistive plate chamber - Meas. Sci. Technol. 17 (2006) 123-127
  //
  //  This function use SI values as input. However it returns charge in picoColomb. The time is returned in seconds!
  //
  //General sequence of the programm:
  //1.  Initialize all the values: Distribut the avalanche along the gap and set the number of primary electrons.
  //2.  Invoke stepping action to simulate the avalanche evolution through the gap.
  //3.  Construct the signal
 
 
 
  int number_of_gaps=12;
  double alpha=144800.; // Townsend coefficient: /m
  double eta=5013;     /// Attachment coefficient:   /m
 
  double thick_gap = 0.00022; // m  //gap thickness
  const int steps =200;//Into how many steps is the gap divided?
  double user_stepwidth=thick_gap/(double)steps;
  double v_drift  = 210.9e3; // m/s 
 
  double E_weight = 1./(6.*0.22e-3+(5.*0.4e-3+2*0.55e-3)/(8.)+2.*0.07e-3/(3.1)); //epsilon_glass= 8  / epsilon_mylar=3.1  /d_mylar=0.07e-3 , d_g= 0.22e-3 ,
  double physics_e= 1.602176462e-19;
  double threshold=93622.7;//Threshold in "number of electrons" when the signal will be detected:  value_in_fc / physics_e
 
  //double alpha=76750;  double eta=8152; double v_drift=163.1e3; //85 kV/cm
  //double alpha=82060; double eta=7728; double v_drift=167.3e3; //87.5 kV/cm
  //double alpha=87560;  double eta=7358; double v_drift=171.1e3;//90 kv/cm
  //double alpha=99020;  double eta=6725; double v_drift=179.7e3;//95 kv/cm
  //double alpha=115200; double eta=6033; double v_drift=191.2e3;//102 kv/cm
  //double alpha=124700; double eta=5683; double v_drift=197.8e3;//106 kv/cm
  //double alpha=134600; double eta=5374; double v_drift=204.3e3;//110 kv/cm
  //double alpha=144500; double eta=5028; double v_drift=210.9e3;//114kv/cm
  //double alpha=154700; double eta=4823; double v_drift=217.5e3;//118kv/cm
 
 
 
 
//Intialize
  vector<int> primary_clusters_in_gap(number_of_gaps);
  for(int i=0;i<number_of_gaps;i++)primary_clusters_in_gap[i]=0;
 
  vector< vector< vector<double> > > step_avalanche(number_of_gaps);
  vector< vector< vector<long int> > > elec_avalche(number_of_gaps);
 
  //How many cluster do we have for the different steps:
  for(int i=0;i<number_of_gaps;i++)
    {
      int counter_avalanches=0;
      for(int j=0;j<avalanche_info->size();j++)
    {
      if(actual_strip==(*avalanche_info)[j].readoutstrip &&  ((i+1)==(*avalanche_info)[j].gap) )
        {counter_avalanches++;}   
    }
 
      primary_clusters_in_gap[i]=counter_avalanches;
      step_avalanche[i].resize(primary_clusters_in_gap[i]);
      elec_avalche[i].resize(primary_clusters_in_gap[i]);
    }
 
 
  for(int j=0;j<number_of_gaps;j++)
    {
      for(int i=0;i<primary_clusters_in_gap[j];i++)
    {
  step_avalanche[j][i].resize(steps+2);
      elec_avalche[j][i].resize(steps+2);
    }
    }
 
  for(int h=0;h<number_of_gaps;h++) {for(int i=0;i<primary_clusters_in_gap[h];i++){for(int j=0;j<steps+2;j++){step_avalanche[h][i][j]=0.;elec_avalche[h][i][j]=0;}}} //Set all arrays equal to zero
 
  vector<int> avalanches_in_gap(number_of_gaps);
 
 
 
 
//Distribute the cluster along the gap and set the number of ions
  for(int h=0;h<number_of_gaps;h++)
    {
 
 
      int i=0;
      for(int j=0;j<avalanche_info->size();j++)
        {
          if(actual_strip==(*avalanche_info)[j].readoutstrip &&  ((h+1)==(*avalanche_info)[j].gap) )
        {
          step_avalanche[h][i][0]=thick_gap-((*avalanche_info)[j].zdistance/m);
          elec_avalche[h][i][0]=(*avalanche_info)[j].ions;
          i++;
        }
        }
 
      avalanches_in_gap[h]=primary_clusters_in_gap[h];
 
    }//close for h (number_of_gaps)
 
  /*
    for(int h=0;h<number_of_gaps;h++)
    {
      for(int j=0;j<primary_clusters_in_gap[h];j++)
    {
      G4cout << "Gap " << h+1 << "   Avalanche  " << j+1 << "   Pos in gap   " <<   step_avalanche[h][j][0] << G4endl;
    }
    }
  */
 
 
 for(int h=0;h<number_of_gaps;h++)
    {
      //      G4cout << "Avalacnhes in gap  " << h+1 << "  are  " << avalanches_in_gap[h] << G4endl;
      long int tot_elec_in_gap=0;
      if(avalanches_in_gap[h]!=0)
    {
 
      int i=0;
 
 
      vector<bool> is_saturated(avalanches_in_gap[h]);
          for(int e=0;e<avalanches_in_gap[h];e++) is_saturated[e]=false;
 
 
 
      do //Steppingaction: Moves the avalanche through the gap!
        {
 
          
          tot_elec_in_gap=0;
          for(int j=0;j<primary_clusters_in_gap[h];j++)
        {
 
          
 
          elec_avalche[h][j][i+1]=adcsignal_simulate_step(elec_avalche[h][j][i],user_stepwidth,alpha,eta,is_saturated[j]);
          step_avalanche[h][j][i+1]= step_avalanche[h][j][i] + user_stepwidth;
          tot_elec_in_gap=tot_elec_in_gap + elec_avalche[h][j][i+1];
          //if((h+1)==12) G4cout<< "Gap  "<<h+1<<" cluster   "<<j+1<<"  =  "<< elec_avalche[h][j][i+1]<<"    Pos   "<<step_avalanche[h][j][i+1]<< "   tot  " <<tot_elec_in_gap<<G4endl;   
 
          if(elec_avalche[h][j][i+1]>1.5e7){is_saturated[j]=true;}else{is_saturated[j]=false;}
        }
          
          
          i++;
 
          
          for(int j=0;j<primary_clusters_in_gap[h];j++)
        {
 
          if((step_avalanche[h][j][i])>=thick_gap) 
            {
              primary_clusters_in_gap[h]=(primary_clusters_in_gap[h]-1);
            }
        }
 
 
        }
      while(step_avalanche[h][0][i]<thick_gap); //If the first avalanche leaves the gap, quit the stepping action
 
 
    }//Close if(avalanches_in_gap!=0)
 }//close for over gaps
 
 
 
 
//Produce signal
 double signal=0;//stores the total sum of the induced signal!
 bool over_threshold = false;
 
 for(int i=0;i<steps;i++)
   {
     for(int h=0;h<number_of_gaps;h++)
       {
     for(int j=0;j<avalanches_in_gap[h];j++)
       {
         signal= signal + elec_avalche[h][j][i];
       }
     }
   
     if(over_threshold == false && signal>threshold) {over_threshold=true; *avalanche_threshold_time=(thick_gap/steps/v_drift*(i+1));}
 
   }//close the for loop over steps.
 
 
 
 
 //Return the value in pC:
 
 //G4cout <<" *avalanche_threshold_time " << *avalanche_threshold_time << G4endl;
 //G4cout << "Induced Charge "<< signal*E_weight*v_drift*physics_e*1e12 << G4endl;
 *induced_charge=signal*E_weight*v_drift*physics_e*1e12   *    thick_gap/steps/v_drift;
 
 
 
}//close simulate_avalanche
 
 
 
 
 
 
 
 
 
 
 
//This fuction return the number of electrons for the next step
long int BesTofDigitizerEcV4_dbs::adcsignal_simulate_step(int n_elec,double stepwidth,double alpha,double eta, bool saturated)
{
 
  if(n_elec<=150 && (saturated == false))
    {
 
      long int n_elec_step =0;
     
      for(int j=0;j<n_elec;j++)
    {
      double s=G4UniformRand();
      
      n_elec_step = n_elec_step + adcsignal_get_n_electron(s, stepwidth, alpha, eta);
      
    }
      return n_elec_step;
    }
 
  if(n_elec>150 && saturated == false)
    {
      double sigma=adcsignal_get_sigma(alpha,eta,stepwidth);
      double nbar = exp((alpha-eta)*stepwidth);
      double my_sigma = G4RandGauss::shoot(0,(sqrt(n_elec)*sigma));
    
      return (long int)(n_elec*nbar+ my_sigma);
    }
 
  if(saturated==true)
    {
 
      return n_elec;
    }//close if(is_saturated==true)
 
}//closesimulate step
 
 
 
 
 
 
 
 
 
//Return the number of electrons within the next step
long int BesTofDigitizerEcV4_dbs::adcsignal_get_n_electron(double s, double stepwidth, double alpha, double eta)
{
 
  double nbar = exp((alpha-eta)*stepwidth);
  double s_border=   eta/alpha*(nbar-1.)/(nbar-eta/alpha);
  double k = eta/alpha;
 
 
  if(s<s_border)
    {
      return 0;
 
    }
  else
    {
      return (long int) (1. + (1/log(((nbar-1)/(nbar-k))) * log(((nbar-k)*(s-1)/(k-1)/nbar))));
    }
 
 
}//close get_n_electron
 
 
double BesTofDigitizerEcV4_dbs::adcsignal_get_sigma(double alpha,double eta, double stepwidth)
{
 
  double k=eta/alpha;
  double nbar = exp((alpha-eta)*stepwidth);
 
  return (sqrt(   (1.+k)/(1.-k)  *   nbar *  (nbar-1.)    )) ;
 
}//close get_sigma
 
 
 
 
 
 
double BesTofDigitizerEcV4_dbs::chargetotime(double charge)
{
 
 
  double time=0.;
  if(charge>=13.) time=(-4.50565/log(charge*0.0812208)+16.6653)*ns;
 
  return time;
 
}//close chargetotime
 
 
 
 
 
 
 
double BesTofDigitizerEcV4_dbs::Get_NINO_pulselength_resolution(double induced_charge_fc)
{
 
  double result = (64.3326*exp(-induced_charge_fc/25.4638 + 0.944184)+19.4846 );
 
  return result;
 
}
 
 
 
double BesTofDigitizerEcV4_dbs::Get_NINO_leadingedge_resolution(double induced_charge_fc)
{
 
  double result =0;
 
  if(induced_charge_fc > 50.)
    {
      result =  67.6737*exp(-induced_charge_fc/50.9995-0.27755)+9.06223;
    }
  else
    {
      result = 176.322-2.98345*induced_charge_fc;
    }
 
  return result;
 
}