Skip to content

Commit

Permalink
Updated DDCaloDigi_BIB for Crilin
Browse files Browse the repository at this point in the history
  • Loading branch information
@sestini88
sestini88 committed May 8, 2024
1 parent 2f5f3ba commit 0cb4695
Showing 1 changed file with 102 additions and 81 deletions.
183 changes: 102 additions & 81 deletions src/DDCaloDigi_BIB.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 ) {
Expand All @@ -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] && theta<theta_bins[i_theta+1] && R>R_bins[i_R] && R<R_bins[i_R+1]) {
if (calhit->getEnergy()>(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);

}

Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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);

Expand Down Expand Up @@ -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<SimCalorimeterHit*>( 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"];
Expand Down Expand Up @@ -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 );
Expand Down Expand Up @@ -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 );
Expand Down

0 comments on commit 0cb4695

Please sign in to comment.