ChargeDividerFP420.cc

Go to the documentation of this file.

// File: ChargeDividerFP420
// Date: 12.2006
// Description: ChargeDividerFP420  for FP420
// Modifications: std::
#include "SimRomanPot/SimFP420/interface/ChargeDividerFP420.h"
#include "DataFormats/GeometryVector/interface/LocalPoint.h"
#include "DataFormats/GeometryVector/interface/LocalVector.h"

//#include "SimDataFormats/TrackingHit/interface/PSimHit.h"
//#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"

using namespace std;
#include<vector>

// unpacking variable - zside - Left or Right planes

//DigitizerFP420::DigitizerFP420(const edm::ParameterSet& conf):conf_(conf),stripDigitizer_(new FP420DigiMain(conf)) 
ChargeDividerFP420::ChargeDividerFP420(double pit, double az420, double azD2, double azD3,int ver){

  verbosity=ver;
  //                           pit - is really moduleThickness here !!!
  if(verbosity>0) {
    std::cout << "ChargeDividerFP420.h: constructor" << std::endl;
    std::cout << "peakMode = " << peakMode << "fluctuateCharge=   "<< fluctuateCharge <<  "chargedivisionsPerHit = "  << chargedivisionsPerHit << "deltaCut=   "<< deltaCut << std::endl;
  }
  
  // Run APV in peak instead of deconvolution mode, which degrades the time resolution
  //  peakMode=true ; //     APVpeakmode
  peakMode=false; //  peakMode=true -->  APVconvolutionmode
  decoMode=false;//  decoMode=true -->  deconvolution mode
  // Enable interstrip Landau fluctuations within a cluster.
  fluctuateCharge=true;   
  
  // Number of segments per strip into which charge is divided during simulation.
  // If large the precision of simulation improves.
  chargedivisionsPerHit=10; // = or =20
  
  // delta cutoff in MeV, has to be same as in OSCAR (0.120425 MeV corresponding // to 100um range for electrons)
  //SimpleConfigurable<double>  ChargeDividerFP420::deltaCut(0.120425,
  deltaCut=0.120425;  //  DeltaProductionCut
  
  pitchcur= pit;// pitchcur - is really moduleThickness here !!!
  
  // but position before Stations:
  z420 = az420;  // dist between centers of 1st and 2nd stations
  zD2 = azD2;  // dist between centers of 1st and 2nd stations
  zD3 = azD3;  // dist between centers of 1st and 3rd stations
  
  
  
  //                                                                                                                           .
  //                                                                                                                           .
  //  -300     -209.2             -150              -90.8                        0                                           +300
  //                                                                                                                           .
  //            X  | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | X                        station                                          .
  //                   8*13.3+ 2*6 = 118.4                                    center                                           .
  //                                                                                                                           .
  //zStationBegPos[0] = -150. - (118.4+10.)/2 + z420; // 10. -arbitrary
  zStationBegPos[0] = -40. + z420; // 5 superplanes per station 79.7mm: -40.- left edge of Station
  zStationBegPos[1] = zStationBegPos[0]+zD2;
  zStationBegPos[2] = zStationBegPos[0]+zD3;
  zStationBegPos[3] = zStationBegPos[0]+2*zD3;
  
}

// Virtual destructor needed.
ChargeDividerFP420::~ChargeDividerFP420() { 
  //  if(verbosity>0) {
  //    std::cout << "Destroying a ChargeDividerFP420" << std::endl;
  //  }
  
}  
CDividerFP420::ionization_type 
ChargeDividerFP420::divide(
               const PSimHit& hit, const double& pitchcur) {
  // !!!
  //                                       pitchcur - is really moduleThickness here !!!
  // !!!
  
  
  
  // sign "-" mean not the same as "+" for middle point !!!
  //  G4ThreeVector direction = hit.getExitLocalP() - hit.getEntryLocalP();
  LocalVector direction = hit.exitPoint() - hit.entryPoint();
  //   G4ThreeVector direction = hit.exitPoint() - hit.entryPoint();
  
  //  LocalVector direction = hit.exitPoint() - hit.entryPoint();
  // direction.mag() - length or (size of path) of the hit; direction/direction.mag() - cosines of direction
  
  if(verbosity>0) {
    std::cout << " CDividerFP420::ChargeDividerFP420:divide: direction= " << direction << std::endl;
    std::cout << " CDividerFP420::ChargeDividerFP420:divide: direction.mag = " << direction.mag() << std::endl;
    std::cout << " obtained as ExitLocalP = " << hit.exitPoint() << "  -  "<<  "  EntryLocalP = "  << hit.entryPoint() << std::endl;
    std::cout << "  pitchcur= " << pitchcur << std::endl;
    std::cout << "  peakMode = " << peakMode << "  decoMode = " << decoMode << "  fluctuateCharge=   "<< fluctuateCharge <<  "  chargedivisionsPerHit = "  << chargedivisionsPerHit << "  deltaCut=   "<< deltaCut << std::endl;
  }
  
  int NumberOfSegmentation =  
    
    //     (int)(1+chargedivisionsPerHit*fabs(direction.x())/pitchcur); // equidistant in X
    //    (int)(1+chargedivisionsPerHit*fabs(direction.z())/pitchcur); // equidistant in Z, but why?
    
    (int)(1+chargedivisionsPerHit*direction.mag()/pitchcur); // equidistant over hit path
  
  
  if(verbosity>0) {
    std::cout << "NumberOfSegmentation= " << NumberOfSegmentation << std::endl;
  }
  
  float eLoss = hit.energyLoss();  // Eloss in GeV
  //  float eLoss = hit.getEnergyLoss();  // Eloss in GeV
  
  if(verbosity>0) {
    std::cout << "CDividerFP420::ChargeDividerFP420:divide: eLoss=  " << eLoss << std::endl;
  }
  
  //
  //     return the energyLoss weighted CR-RC shape peaked at t0.(PeakShape)
  //     return the energyLoss weighted with a gaussian centered at t0 (DeconvolutionShape)
  float decSignal = TimeResponse(hit);
  if(verbosity>0) {
    std::cout << "CDividerFP420::ChargeDividerFP420:divide: decSignal=  " << decSignal << std::endl;
  }
  
  ionization_type _ionization_points;
  
  _ionization_points.resize(NumberOfSegmentation);
  
  float energy;
  
  // Fluctuate charge in track subsegments
  float* eLossVector = new float[NumberOfSegmentation];
  
  
  if(verbosity>0) {
    std::cout << "CDividerFP420::ChargeDividerFP420:divide:  resize done; then, fluctuateCharge ? = " << fluctuateCharge << std::endl;
  }
  if( fluctuateCharge ) {
    //    int pid = hit.getParticleType();
    //    float momentum = hit.getPabs();
    int pid = hit.particleType();
    float momentum = hit.pabs();
    float length = direction.mag();  // length or (size of path) of the hit;
    
    if(verbosity>0) {
      std::cout << "pid= " << pid << "momentum= " << momentum << "eLoss= " << eLoss << "length= " << length << std::endl;
    }
    fluctuateEloss(pid, momentum, eLoss, length, NumberOfSegmentation, eLossVector);   
  }
  
  for ( int i = 0; i != NumberOfSegmentation; ++i) {
    if( fluctuateCharge ) {
      energy=eLossVector[i]*decSignal/eLoss;
      EnergySegmentFP420 edu(energy,hit.entryPoint()+float((i+0.5)/NumberOfSegmentation)*direction);//take energy value from vector eLossVector  
      //   EnergySegmentFP420 edu(energy,hit.getEntryLocalP()+float((i+0.5)/NumberOfSegmentation)*direction);//take energy value from vector eLossVector  
      _ionization_points[i] = edu; //save
    }else{
      energy=decSignal/float(NumberOfSegmentation);
      EnergySegmentFP420 edu(energy,hit.entryPoint()+float((i+0.5)/NumberOfSegmentation)*direction);//take energy value from eLoss average over n.segments 
      // EnergySegmentFP420 edu(energy,hit.getEntryLocalP()+float((i+0.5)/NumberOfSegmentation)*direction);//take energy value from eLoss average over n.segments 
      _ionization_points[i] = edu; //save
    }
  }
  
  if(verbosity>0) {
    std::cout << "CDividerFP420::ChargeDividerFP420:divide:  !!!  RESULT !!!" << std::endl;
    std::cout <<  " _ionization_points size = " << _ionization_points.size() << std::endl;
    for(unsigned int i = 0; i < _ionization_points.size(); ++i ) {
      std::cout <<  " eLossVector[i] i = " << i << eLossVector[i] << std::endl;
    }
  }
  
  delete[] eLossVector;
  return _ionization_points;
}

void ChargeDividerFP420::fluctuateEloss(int pid, float particleMomentum, 
                    float eloss, float length, 
                    int NumberOfSegs,float elossVector[]) {
  
  if(verbosity>0) {
    std::cout << "fluctuateEloss: eloss=  " << eloss << "length=  " << length << "NumberOfSegs=  " << NumberOfSegs << std::endl;
  }
  
  //  double particleMass = 139.57; // Mass in MeV, Assume pion
  double particleMass = 938.271; // Mass in MeV, Assume proton   ----  AZ
  //  if( particleTable->getParticleData(pid) ) {  // Get mass from the PDTable
  //    particleMass = 1000. * particleTable->getParticleData(pid)->mass(); //Conv. GeV to MeV
  //  }
  pid = abs(pid);
  if(pid==11) particleMass = 0.511;         // Mass in MeV
  else if(pid==13) particleMass = 105.658;
  else if(pid==211) particleMass = 139.570;
  //  else if(pid==2212) particleMass = 938.271;
  
  float segmentLength = length/NumberOfSegs;
  
  // Generate charge fluctuations.
  float de=0.;
  float sum=0.;
  double segmentEloss = (1000.*eloss)/NumberOfSegs; //eloss in MeV
  if(verbosity>0) {
    std::cout << "segmentLength=  " << segmentLength << "segmentEloss=  " << segmentEloss << std::endl;
  }
  
  for (int i=0;i<NumberOfSegs;++i) {
    // The G4 routine needs momentum in MeV, mass in Mev, delta-cut in MeV,
    // track segment length in mm(!!!), segment eloss in MeV 
    // Returns fluctuated eloss in MeV
    //    double deltaCutoff = deltaCut.value(); // the cutoff is sometimes redefined inside, so fix it.
    double deltaCutoff = deltaCut;
    de = fluctuate.SampleFluctuations(double(particleMomentum*1000.),
                      particleMass, deltaCutoff, 
                      double(segmentLength),
                      segmentEloss )/1000.; //convert to GeV
    elossVector[i]=de;
    sum +=de;
  }
  
  if(verbosity>0) {
    std::cout << "sum=  " << sum << std::endl;
  }
  if(sum>0.) {  // If fluctuations give eloss>0.
    // Rescale to the same total eloss
    float ratio = eloss/sum;
    for (int ii=0;ii<NumberOfSegs;++ii) elossVector[ii]= ratio*elossVector[ii];
  } else {  // If fluctuations gives 0 eloss
    float averageEloss = eloss/NumberOfSegs;
    for (int ii=0;ii<NumberOfSegs;++ii) elossVector[ii]= averageEloss; 
  }
  return;
}

float ChargeDividerFP420::TimeResponse( const PSimHit& hit ) {
  if (peakMode) {
    
    if(verbosity>0) {
      std::cout << "ChargeDividerFP420:TimeResponse: call of PeakShape" << std::endl;
    }
    return this->PeakShape( hit );
  } else if (decoMode) {
    
    if(verbosity>0) {
      std::cout << "ChargeDividerFP420:TimeResponse: call of DeconvolutionShape" << std::endl;
    }
    return this->DeconvolutionShape( hit );
  } else {
    
    if(verbosity>0) {
      std::cout << "ChargeDividerFP420:TimeResponse: no any  Shape" << std::endl;
    }
    //    return hit.getEnergyLoss();
    return hit.energyLoss();
    
  }
}
float ChargeDividerFP420::PeakShape(const PSimHit& hit){
  // 
  // Aim:     return the energyLoss weighted CR-RC shape peaked at t0.
  // 
  
  //   float xEntry = hit.getX() - hit.getVx();
  //  float yEntry = hit.getY() - hit.getVy();
  //  float zEntry = hit.getZ() - hit.getVz();
  float xEntry = 0.5;
  float yEntry = 0.5;
  float zEntry = 1000.;
  
  // unsigned int unitID = hit.getUnitID();
  unsigned int unitID = hit.detUnitId();
  //      int sScale = 20;
  int det, zside, sector, zmodule;
  FP420NumberingScheme::unpackFP420Index(unitID, det, zside, sector, zmodule);
  // intindex is a continues numbering of FP420
  //      int zScale=2;  unsigned int intindex = sScale*(sector - 1)+zScale*(zmodule - 1)+zside;
  //     int zScale=10; unsigned int intindex = sScale*(sector - 1)+zScale*(zside - 1)+zmodule;
  
  float RRR      = sqrt(xEntry*xEntry + yEntry*yEntry + zEntry*zEntry);
  float costheta = zEntry / RRR ;
  //  float theta    = acos(min(max(costheta,float(-1.)),float(1.)));
  //  float dist = hit.det().position().mag();
  //  float dist = hit.localPosition().mag();//AZ
  //  float dist = hit.getEntry().mag();
  //  float dist = hit.getEntryLocalP().mag();
  float dist     = (zStationBegPos[sector-1] - 420000.) / costheta;
  //  float dist     = (zStationBegPos[sector-1] - hit.getVz()) / costheta;
  dist     = dist/10.;// mm  --> cm as light velocity = 30 cm/ns
  
  if(verbosity>0) {
    std::cout << "sector=" << sector << std::endl;
    std::cout << "zmodule=" << zmodule << std::endl;
    std::cout << "zStationBegPos[sector-1]=" << zStationBegPos[sector-1] << std::endl;
    std::cout << "RRR=" << RRR << std::endl;
    std::cout << "costheta=" << costheta << std::endl;
    std::cout << "unitID=" << unitID << std::endl;
    //std::cout << "thetaEntry=" << thetaEntry << std::endl;
    //std::cout << "my theta=" << theta*180./3.1415927 << std::endl;
    std::cout << "dist found =" << dist << std::endl;
  }
  
  // Time when read out relative to time hit produced.
  float t0 = dist/30.;  // light velocity = 30 cm/ns
  float SigmaShape = 52.17; 
  //  float tofNorm = (hit.getTof() - t0)/SigmaShape;
  float tofNorm = (hit.tof() - t0)/SigmaShape;
  
  float readTimeNorm = -tofNorm;
  // return the energyLoss weighted with CR-RC shape peaked at t0.
  
  if(verbosity>0) {
    std::cout << "ChargeDividerFP420:PeakShape::dist=" << dist << std::endl;
    std::cout << "t0=" <<t0  << std::endl;
    std::cout << "hit.getTof()=" << hit.tof() << std::endl;
    std::cout << "tofNorm=" << tofNorm << std::endl;
    std::cout << "1 + readTimeNorm=" << 1 + readTimeNorm << std::endl;
    std::cout << "hit.getEnergyLoss()=" << hit.energyLoss()  << std::endl;
    std::cout << "(1 + readTimeNorm)*exp(-readTimeNorm)=" << (1 + readTimeNorm)*exp(-readTimeNorm)  << std::endl;
    std::cout << "return=" << hit.energyLoss()*(1 + readTimeNorm)*exp(-readTimeNorm)  << std::endl;
  }
  if (1 + readTimeNorm > 0) {
    //    return hit.energyLoss()*(1 + readTimeNorm)*exp(-readTimeNorm);
    return hit.energyLoss()*(1 + readTimeNorm)*exp(-readTimeNorm);
    //  return hit.getEnergyLoss()*(1 + readTimeNorm)*exp(-readTimeNorm);
  } else {
    return 0.;
  }
}

float ChargeDividerFP420::DeconvolutionShape(const PSimHit& hit){
  // 
  // Aim:     return the energyLoss weighted with a gaussian centered at t0 
  // 
  //   float xEntry = hit.getX() - hit.getVx();
  //  float yEntry = hit.getY() - hit.getVy();
  //  float zEntry = hit.getZ() - hit.getVz();
  float xEntry = 0.5;
  float yEntry = 0.5;
  float zEntry = 1000.;
  
  // unsigned int unitID = hit.getUnitID();
  unsigned int unitID = hit.detUnitId();
  //      int sScale = 20;
  int det, zside, sector, zmodule;
  FP420NumberingScheme::unpackFP420Index(unitID, det, zside, sector, zmodule);
  // intindex is a continues numbering of FP420
  //      int zScale=2;  unsigned int intindex = sScale*(sector - 1)+zScale*(zmodule - 1)+zside;
  //     int zScale=10; unsigned int intindex = sScale*(sector - 1)+zScale*(zside - 1)+zmodule;
  
  float RRR      = sqrt(xEntry*xEntry + yEntry*yEntry + zEntry*zEntry);
  float costheta = zEntry / RRR ;
  //  float theta    = acos(min(max(costheta,float(-1.)),float(1.)));
  //  float dist = hit.det().position().mag();
  //  float dist = hit.localPosition().mag();//AZ
  //  float dist = hit.getEntry().mag();
  //  float dist = hit.getEntryLocalP().mag();
  float dist     = (zStationBegPos[sector-1] - 420000.) / costheta;
  //  float dist     = (zStationBegPos[sector-1] - hit.getVz()) / costheta;
  dist     = dist/10.;// mm  --> cm as light velocity = 30 cm/ns
  
  if(verbosity>0) {
    std::cout << "sector=" << sector << std::endl;
    std::cout << "zmodule=" << zmodule << std::endl;
    std::cout << "zStationBegPos[sector-1]=" << zStationBegPos[sector-1] << std::endl;
    std::cout << "RRR=" << RRR << std::endl;
    std::cout << "costheta=" << costheta << std::endl;
    std::cout << "unitID=" << unitID << std::endl;
    //std::cout << "thetaEntry=" << thetaEntry << std::endl;
    //std::cout << "my theta=" << theta*180./3.1415927 << std::endl;
    std::cout << "dist found =" << dist << std::endl;
  }
  
  float t0 = dist/30.;  // light velocity = 30 cm/ns
  float SigmaShape = 12.; 
  //fun/pl 1*exp(-0.5*((0.1/30-x)/0.1)**2) 0. 0.08
  //  float SigmaShape = 22.; 
  //  float tofNorm = (hit.tof() - t0)/SigmaShape;
  float tofNorm = (hit.tof() - t0)/SigmaShape;
  // Time when read out relative to time hit produced.
  float readTimeNorm = -tofNorm;
  // return the energyLoss weighted with a gaussian centered at t0 
  //  return hit.energyLoss()*exp(-0.5*readTimeNorm*readTimeNorm);
  
  if(verbosity>0) {
    std::cout << "ChargeDividerFP420:DeconvolutionShape::dist=" << dist << std::endl;
    std::cout << "t0=" <<t0  << std::endl;
    std::cout << "hit.getTof()=" << hit.tof() << std::endl;
    std::cout << "tofNorm=" << tofNorm << std::endl;
    std::cout << "hit.getEnergyLoss()=" << hit.energyLoss()  << std::endl;
    std::cout << "exp(-0.5*readTimeNorm*readTimeNorm)=" << exp(-0.5*readTimeNorm*readTimeNorm)  << std::endl;
    std::cout << "return=" << hit.energyLoss()*exp(-0.5*readTimeNorm*readTimeNorm)  << std::endl;
  }
  return hit.energyLoss()*exp(-0.5*readTimeNorm*readTimeNorm);
  //  return hit.getEnergyLoss();
}