diff --git a/src/DDCaloDigi_BIB.cc b/src/DDCaloDigi_BIB.cc index ea0c7bc..b385ec1 100644 --- a/src/DDCaloDigi_BIB.cc +++ b/src/DDCaloDigi_BIB.cc @@ -51,6 +51,47 @@ dd4hep::rec::LayeredCalorimeterData * getExtension(unsigned int includeFlag, uns DDCaloDigi_BIB aDDCaloDigi_BIB ; +double GetThresholdCrilin(double x,double y,double z,int nsigma){ + + float r = sqrt(x*x+y*y); + + double phi_rel= acos(x/r); + if (y<0) phi_rel = -phi_rel+2*3.14159; + int n = phi_rel/(2*3.14159/24); + phi_rel = phi_rel - n*2*3.14159/24; + double x_rel = r*cos(phi_rel); + + double emean[5][3]={161.1, 125.4, 68.55, + 63.75, 50.38, 27.59, + 26.02, 20.85, 11.65, + 14.29, 12.01, 7.568, + 10.34, 9.356, 6.108}; + + double estd[5][3]={62.45, 53.89, 34.94, + 51.8, 43.57, 28, + 21.77, 18.92, 13.58, + 15.27, 13.78, 10.93, + 13.05, 12.22, 9.632}; + + int nx; + + if (x_rel<1550.) nx=0; + if (x_rel>1550. && x_rel< 1580.) nx=1; + if (x_rel>1580. && x_rel< 1630.) nx=2; + if (x_rel>1630. && x_rel< 1670.) nx=3; + if (x_rel>1670.) nx=4; + + int nz = abs(z/(2210./3.)); + + if (nz>2) nz=2; + + double eth=(emean[nx][nz]+nsigma*estd[nx][nz])*0.001; + + return eth; + + } + + // helper struct for string comparision struct XToLower{ int operator() ( int ch ) { @@ -59,100 +100,72 @@ struct XToLower{ }; //Thresholds optimized for BIB at 1.5 TeV -bool applyDifferentialThresholdEcalBIB(CalorimeterHitImpl * calhit, bool useCLIC, bool useCrilin) { +bool applyDifferentialThresholdEcalBIB(CalorimeterHitImpl * calhit, bool useCrilinBarrel) { bool pass = false; - if (useCLIC==false && useCrilin==false) pass = true; + if (useCrilinBarrel==true){ - if (useCLIC==true){ - float x = calhit->getPosition()[0]; float y = calhit->getPosition()[1]; float z = calhit->getPosition()[2]; - float d = sqrt(x*x+y*y+z*z); - float R = sqrt(x*x+y*y); - float theta = acos(z/d); - - float theta_bins[3+1] = {1.57-1.02,1.57-0.27,1.57+0.27,1.57+1.02}; - float R_bins[5+1] = {1500,1570,1600,1650,1710,1770}; - - float th_barrel[5][3] = {18.2731, 17.1456, 18.4469, - 24.7357, 21.9555, 24.6477, - 33.8206, 30.3771, 33.8971, - 43.646, 42.9506, 43.3244, - 45.6661, 47.0891, 46.0237}; - - float std_barrel[5][3] = {32.9455, 31.5414, 34.2989, - 55.6836, 48.5732, 55.6375, - 77.1442, 71.8559, 76.0541, - 93.3592, 94.8317, 91.6914, - 96.7721, 96.9367, 94.6892}; - - for (int i_theta=0; i_theta<3; i_theta++){ - for (int i_R=0; i_R<5; i_R++){ - if (theta>theta_bins[i_theta] && thetaR_bins[i_R] && RgetEnergy()>(th_barrel[i_R][i_theta]+2*std_barrel[i_R][i_theta])*0.001) {pass = true; - calhit->setEnergy(calhit->getEnergy()-th_barrel[i_theta][i_R]*0.001); - } - } - } - } + float t = calhit->getTime(); - } + //cout << calhit->getTime() << endl; - if (useCrilin==true){ + float r = sqrt(x*x+y*y); - float x = calhit->getPosition()[0]; - float y = calhit->getPosition()[1]; - float z = calhit->getPosition()[2]; + double phi_rel= acos(x/r); + if (y<0) phi_rel = -phi_rel+2*3.14159; + int n = phi_rel/(2*3.14159/24); + phi_rel = phi_rel - n*2*3.14159/24; + double x_rel = r*cos(phi_rel); - //Rotation of the dodecaedra + int nx; - float phi = acos(x/sqrt(x*x+y*y)); + if (x_rel<1550.) nx=0; + if (x_rel>1550. && x_rel< 1580.) nx=1; + if (x_rel>1580. && x_rel< 1630.) nx=2; + if (x_rel>1630. && x_rel< 1670.) nx=3; + if (x_rel>1670.) nx=4; - if (y<0) phi = -phi; + int nz = abs(z/(2210./3.)); - float xphi = phi + 3.14159/12.; + if (nz>2) nz=2; - if (xphi<0) xphi = 2*3.14159+xphi; + //Efficiency for charge collection depending from integration time + + double eff = 1.0; - int nphi = xphi*6./3.14159; + double t_max = 25; //integration time + double t_int = t_max-t; - float delta_phi = 3.14159/2.-2*3.14159/12.*nphi; + if (t_int<15) eff=0.67; + if (t_int>=15 && t_int<25) eff=0.73; + if (t_int>=25 && t_int<35) eff=0.82; + if (t_int>=35 && t_int<45) eff=0.88; + if (t_int>=45 && t_int<55) eff=0.91; + if (t_int>=55 && t_int<70) eff=0.93; + if (t_int>=70 && t_int<90) eff=0.96; + if (t_int>=90 && t_int<110) eff=0.99; + if (t_int>=110) eff=1.0; - //float xprime = x*cos(delta_phi)-y*sin(delta_phi); - float yprime = x*sin(delta_phi)+y*cos(delta_phi); + //double eff_BIB = 0.99; //t=100 ns + //double eff_BIB = 0.91; //t=50 ns + double eff_BIB = 0.78; //t=25 ns - float m_energyc[5] = {14.1607,17.6719,23.0927,24.9074,28.0001}; - float m_energyf[5] = {20.3787,18.6048,19.045,20.4087,23.8197}; - float s_energyc[5] = {18.0425,24.4413,36.4643,33.8902,36.6783}; - float s_energyf[5] = {20.9336,20.1998,22.8853,26.9295,35.3377}; - float tenergyc[5]; - float tenergyf[5]; + double eth=GetThresholdCrilin(x,y,z,3); - for (int k=0; k<5; k++){ - tenergyc[k]=m_energyc[k]+2*s_energyc[k]; - tenergyf[k]=m_energyf[k]+2*s_energyf[k]; - } - if (z>-500 && z<500){ //Central region - if (yprime<1550) if (calhit->getEnergy()>tenergyc[0]*0.001) pass = true; - if (yprime>1550 && yprime<1600 ) if (calhit->getEnergy()>tenergyc[1]*0.001) pass = true; - if (yprime>1600 && yprime<1630) if (calhit->getEnergy()>tenergyc[2]*0.001) pass = true; - if (yprime>1630 && yprime<1660) if (calhit->getEnergy()>tenergyc[3]*0.001) pass = true; - if (yprime>1660) if (calhit->getEnergy()>tenergyc[4]*0.001) pass = true; - } + //if (eff*calhit->getEnergy()>eff_BIB*eth && abs(t)<0.250) cout << "energy = " << eff_BIB*calhit->getEnergy()*1000 << " threshold = " << eth*1000 << " nx = " << nx << " nz = " << nz << " z = " << abs(z) << endl; + + if (eff*calhit->getEnergy()>eff_BIB*eth && abs(t)<0.250) pass = true; - if (z<-500 || z>500){ //Forward region - if (yprime<1550) if (calhit->getEnergy()>tenergyf[0]*0.001) pass = true; - if (yprime>1550 && yprime<1600 ) if (calhit->getEnergy()>tenergyf[1]*0.001) pass = true; - if (yprime>1600 && yprime<1630) if (calhit->getEnergy()>tenergyf[2]*0.001) pass = true; - if (yprime>1630 && yprime<1660) if (calhit->getEnergy()>tenergyf[3]*0.001) pass = true; - if (yprime>1660) if (calhit->getEnergy()>tenergyf[4]*0.001) pass = true; - } + cout << nx << " " << nz << " " << calhit->getEnergy()*1000 << " " << x << " " << y << " " << z << " " << t << endl; + double energy = eff*calhit->getEnergy(); + if (energy>0) calhit->setEnergy(energy); } @@ -496,14 +509,9 @@ DDCaloDigi_BIB::DDCaloDigi_BIB() : Processor("DDCaloDigi_BIB") { _deadCellEcal_keep, (bool)false); - registerProcessorParameter("ECAL_use_CLIC" , - "Use CLIC BIB subtraction " , - _useCLIC, - (bool)false); - registerProcessorParameter("ECAL_use_Crilin" , - "Use Crilin BIB subtraction " , - _useCrilin, + "Use Crilin Barrel BIB subtraction " , + _useCrilinBarrel, (bool)false); registerProcessorParameter("ECAL_strip_absorbtionLength", @@ -793,7 +801,7 @@ void DDCaloDigi_BIB::processRunHeader( LCRunHeader* /*run*/) { } void DDCaloDigi_BIB::processEvent( LCEvent * evt ) { - + // create the output collections LCCollectionVec *relcol = new LCCollectionVec(LCIO::LCRELATION); @@ -854,12 +862,25 @@ void DDCaloDigi_BIB::processEvent( LCEvent * evt ) { int numElements = col->getNumberOfElements(); streamlog_out ( DEBUG ) << colName << " number of elements = " << numElements << endl; + cout << "//////////////////Number of elements = " << numElements << endl; //Lorenzo + for (int j(0); j < numElements; ++j) { SimCalorimeterHit * hit = dynamic_cast( col->getElementAt( j ) ) ; float energy = hit->getEnergy(); + float hitx = hit->getPosition()[0]; + float hity = hit->getPosition()[1]; + float hitz = hit->getPosition()[2]; + + double eth=GetThresholdCrilin(hitx,hity,hitz,3); + + //Correction factors due to integration time effects + double eff_min=0.67; + double eff_BIB=0.78; + // apply threshold cut - if (energy > _thresholdEcal) { + //if (eff_min*energy>eff_BIB*eth) { + if (1) { int cellid = hit->getCellID0(); int cellid1 = hit->getCellID1(); int layer = idDecoder(hit)["layer"]; @@ -1011,7 +1032,7 @@ void DDCaloDigi_BIB::processEvent( LCEvent * evt ) { calhit->setPosition(hit->getPosition()); calhit->setType( CHT( CHT::em, CHT::ecal , caloLayout , layer ) ); calhit->setRawHit(hit); - if (applyDifferentialThresholdEcalBIB(calhit,_useCLIC,_useCrilin)){ + if (applyDifferentialThresholdEcalBIB(calhit, _useCrilinBarrel) || (caloLayout!=CHT::barrel)){ ecalcol->addElement(calhit); LCRelationImpl *rel = new LCRelationImpl(calhit,hit,1.0); relcol->addElement( rel ); @@ -1044,7 +1065,7 @@ void DDCaloDigi_BIB::processEvent( LCEvent * evt ) { calhit->setPosition(hit->getPosition()); calhit->setType( CHT( CHT::em, CHT::ecal , caloLayout , layer ) ); calhit->setRawHit(hit); - if (applyDifferentialThresholdEcalBIB(calhit,_useCLIC,_useCrilin)){ + if (applyDifferentialThresholdEcalBIB(calhit, _useCrilinBarrel) || (caloLayout!=CHT::barrel)){ ecalcol->addElement(calhit); LCRelationImpl *rel = new LCRelationImpl(calhit,hit,1.0); relcol->addElement( rel );