#include "KITutil_LDC00.h"
// Will be called by KIT.cc
// requires inout from myPGDB & MarlinMath
//Shits = superHits
void CreateAllShits2(LCCollection* colt,CellIDDecoder& id,vector* calo, int debug)
{
unsigned int n_elements = colt->getNumberOfElements();
int superHitCounterSampling1 = 0;
int superHitCounterSampling2 = 0;
for ( unsigned int elementCounter=0; elementCounter(colt->getElementAt(elementCounter));
unsigned int layer = (unsigned int)id(caloHit)["K-1"];
Superhit2 *superHit;
//requires information about detector geometry from Physics Geometry database myPGDB
//for 2 different ECAL samplings in barrel ECAL1 & ECAL2
if ( layer<=myPGDB[myPGdb::ECAL1_BAR].max_lay ) {//maxlayer from gear-file
superHit = new Superhit2(myPGDB[myPGdb::ECAL1_BAR].mip_whole,caloHit); //mip_whole is hardcoded via e_coeff in myPGDB == energy of a mip?
superHitCounterSampling1++;
} else {
superHit = new Superhit2(myPGDB[myPGdb::ECAL2_BAR].mip_whole,caloHit);
superHitCounterSampling2++;
}
calo[0].push_back(superHit);
}
if (debug > 3 ) {
std::cout<<"DEBUG KITutil: * from "<< n_elements <<" calorimeter hits created "
<< superHitCounterSampling1 <<" superhits in sampling 1 and "
<< superHitCounterSampling2 <<" superhits in sampling 2"<* calo,CellIDDecoder& id,unsigned int nelem,int Ncut,int debug)
{
if ( debug > 3 )
std::cout<<"DEBUG KITutil: - precalculation for "< tmpsh[9]; //fields 0 to 7 are superhits in stave 0 to 7, field 8 holds the endcaps
int barrelCounter = 0;
int endcapCounter = 0;
for ( unsigned int hitCounter=0; hitCounterchit;
int stave = id(chh)["S-1"]; //which octant
int module = id(chh)["M"]; //which module
int layer = id(chh)["K-1"]; //which layer
// some remaining superhit members initialised
(calo[0])[hitCounter]->S=stave;
(calo[0])[hitCounter]->M=module;
(calo[0])[hitCounter]->K=layer;
(calo[0])[hitCounter]->isCalorimeter=1; //is ecal
if ( module !=0 && module!=6 ) { //if hit not in endcaps
if ( layer!=0 ) { //and not in first layer
tmpsh[stave].push_back((calo[0])[hitCounter]);
} else {
//for hits in first layer, phi = angle between hit and y axis
//the following section determines if the first layer hit is a neighbour of the next stave.
//If so it is assigned to both stave vectors
double phi = atan2(-(calo[0])[hitCounter]->chit->getPosition()[0],(calo[0])[hitCounter]->chit->getPosition()[1]);
if ( phi< 0.0 ) {
if ( stave!=0 ) phi = 2.0*M_PI + phi; //for staves 1 - 7; M_PI = pi = 3.14
else phi = fabs(phi);
}
double diff = phi - M_PI_4*stave; //M_PI_4 = pi/4 = 0.79
// 45 degrees time is the angle between to corners of the octagon
// now the hit lies in stave 0 and if fabs(phi)<45 DEG, it is not overlapping, see below
//The hit is written to its original stave anyway
tmpsh[stave].push_back((calo[0])[hitCounter]);
if ( fabs(diff) stored a 2nd time ?
tmpsh[stave].push_back((calo[0])[hitCounter]);
} else { // If its in the overlap area, put it in both stove collections
(calo[0])[hitCounter]->connect = true;
tmpsh[stave].push_back((calo[0])[hitCounter]);
int stavem1 = stave - 1;
// store hit in additional stave
if ( stavem1>=0 )
tmpsh[stavem1].push_back((calo[0])[hitCounter]);
else
tmpsh[7].push_back((calo[0])[hitCounter]);
}
} // special treatment for first layer (= layer 0)
barrelCounter++;
} // hit in barrel
else { //hits in endcap
tmpsh[8].push_back((calo[0])[hitCounter]);
endcapCounter++;
}
}//all elements
if ( debug > 3 )
std::cout<<"DEBUG KITutil: * found "<top = calo[0][ii]->neighbours.size();
for ( unsigned int ii=0; ii<(calo[0]).size(); ii++) {
if(debug > 3)
std::cout<<"DEBUG KITutil: - found hit with "<<(calo[0])[ii]->top<<" neighbours"<top !=0 ) {//if hit has neighbours
if ( (calo[0])[ii]->top >= 1 && (calo[0])[ii]->top < Ncut ) {
(calo[2]).push_back((calo[0])[ii]);// if hit with less than required neighbours -> no cluster
} else { //for topological cleaning
(calo[4]).push_back((calo[0])[ii]); // hit with required neighbours -> topological cleaning possible
}
} else { // if calo[0] hit has no neighbours at all
(calo[6]).push_back((calo[0])[ii]); //-> no cluster
}
}
}
void Precalc2 ( vector< Superhit2* >& shvec,double r, double z, double cell, double dist,bool isCalorimeter,int stave,int module,CellIDDecoder& id,int debug ) {
unsigned int nPoints = shvec.size();
ANNpointArray dataPoints;
ANNpoint queryPoint; //Annpoint: double[3]
ANNidxArray nnIndex;
ANNdistArray distance;
ANNkd_tree* kdTree;
queryPoint = annAllocPt(3); // 3d point
dataPoints = annAllocPts(nPoints,3);
nnIndex = new ANNidx[36];
distance = new ANNdist[36];
float xxx[3]; // transformed position in GridTransform
float Ecalradius = r;
float Ecalhalfz = z;
float Ecalcellsize = cell;
if ( debug>3 )
std::cout<<"DEBUG KITutil: shifting hits to new positions"<chit,Ecalradius,Ecalhalfz,Ecalcellsize,xxx,isCalorimeter,stave,module,id);
//dataPoints are filled with transformed coordinates for search alg.
dataPoints[hitCounter][0] = xxx[0];
dataPoints[hitCounter][1] = xxx[1];
dataPoints[hitCounter][2] = xxx[2];
//point[] is the transformed position in class SuperHit
shvec[hitCounter]->point[0] = xxx[0];
shvec[hitCounter]->point[1] = xxx[1];
shvec[hitCounter]->point[2] = xxx[2];
}
//search for nearest neighbours with ANN code, assumed to work correctly see ANN.cpp
kdTree = new ANNkd_tree(dataPoints,nPoints,3);
ANNdist trd = dist;
for ( unsigned int hitCounter=0; hitCounterannkFRSearch(dataPoints[hitCounter],trd,26,nnIndex,distance,0.0);
if ( N_NearestNeighbours>26 )
N_NearestNeighbours = 26;
if( isCalorimeter ) { // inside calorimeter: isCalorimeter(ecal)=1, isCalorimeter(hcal)=2
if( shvec[hitCounter]->connect ) { // hit in barrel and connected to two staves
for ( unsigned int neighbourCounter=0; neighbourCounterneighbours.end()
== find(shvec[hitCounter]->neighbours.begin(),shvec[hitCounter]->neighbours.end(), shvec[nnIndex[neighbourCounter]]) ) {
//add neighbour to SuperHit
shvec[hitCounter]->neighbours.push_back(shvec[nnIndex[neighbourCounter]]);
shvec[hitCounter]->top++;
}
if ( shvec[nnIndex[neighbourCounter]]->neighbours.end()
== find(shvec[nnIndex[neighbourCounter]]->neighbours.begin(),shvec[nnIndex[neighbourCounter]]->neighbours.end(), shvec[hitCounter]) ) {
shvec[nnIndex[neighbourCounter]]->neighbours.push_back(shvec[hitCounter]);
shvec[nnIndex[neighbourCounter]]->top++;
}
}//all neighbours
} else if ( N_NearestNeighbours !=0 ) {
shvec[hitCounter]->top = N_NearestNeighbours;
for ( unsigned int neighbourCounter=0; neighbourCounterneighbours.push_back( shvec[nnIndex[neighbourCounter]]);
}
} else if ( N_NearestNeighbours !=0 ) {
shvec[hitCounter]->top=N_NearestNeighbours;
for ( unsigned int neighbourCounter=0; neighbourCounterneighbours.push_back( shvec[nnIndex[neighbourCounter]]);
}
}//all hits
delete [] nnIndex;
delete [] distance;
annDeallocPts(dataPoints);
delete kdTree;
annClose();
}
//called by precalc2
void GridTransform2 ( CalorimeterHit* caloHit,float& radius, float& halfz, float& cellsize,float*X,bool isCalorimeter,int stave,int module,CellIDDecoder& id ) {
unsigned int layer;
double coss;
double sinn;
if ( !isCalorimeter ) { // !calorimeter -> tracker??
stave = id(caloHit)["S-1"];
module = id(caloHit)["M"]; //module is back to 6 or 1 for endcap?
layer = id(caloHit)["K-1"];
} else {
int stave2 = id(caloHit)["S-1"];
// The following calcs are done for the secondary stave the hit is assigned to, where it floats "freely"in space until now.
if ( stave2!=stave ) {
//define rotation matrix, hits are rotated to stave 0
coss = cos(M_PI_4*stave);
sinn = sin(M_PI_4*stave);
//rotate stave to position 0 (only y-axis)
float tmp_y = coss*caloHit->getPosition()[1] - sinn*caloHit->getPosition()[0];
//divide distance of hit from ecal face by sampling layer thickness->layer
layer = (unsigned int) ((tmp_y-radius)/myPGDB[myPGdb::ECAL1_BAR].sampling); // largest int number: 2.3->2, therefore n_sampl-1 in next line
if ( layer> myPGDB[myPGdb::ECAL1_BAR].n_sampl - 1 ) { // last layer of first ECAl structure
layer = myPGDB[myPGdb::ECAL1_BAR].n_sampl - 1
+ (unsigned int)(( tmp_y-radius- myPGDB[myPGdb::ECAL1_BAR].n_sampl*
myPGDB[myPGdb::ECAL1_BAR].sampling)/myPGDB[myPGdb::ECAL2_BAR].sampling);
}
} else { // (stave2 == stave) original layer for hits in their original stave
layer = id(caloHit)["K-1"];
}
//Now we know of every hit's (new) layer
}
// out of if
coss = cos(M_PI_4*stave);
sinn = sin(M_PI_4*stave);
if ( module !=0 && module!=6 ) {
float tmp_x = coss*caloHit->getPosition()[0] + sinn*caloHit->getPosition()[1]; // rotate, find new y position
float tmp_y = radius + (layer+1)*cellsize; // center of layer position.
// y using the variable cellsize, the hits are equally spaced in 3d
//rotate back
X[0] = coss*tmp_x - sinn*tmp_y;
X[1] = coss*tmp_y + sinn*tmp_x;
X[2] = caloHit->getPosition()[2];
} else { // endcap
X[0] = caloHit->getPosition()[0];
X[1] = caloHit->getPosition()[1];
X[2] = (-1+module/3)*(halfz+layer*cellsize);
}
}
// searching for cluster cores
// wil be called by KIT
// secal1 = superhits in calo[4] with topological cleaning
// secal1 = superhits in calo[0] without topological cleaning
void FindCores2 (Shitvec2* secal1, Tmpclvec2* clusterVector , vector * photonSeed2, unsigned int Nlevels, vector mipstep, CoreCut Ccut,int debug )
{
if (debug >2 )
std::cout<<" ***** started main loop: searching for cluster cores in "<size()<<" hits ***"< test;
int psl_plun_global=0;
double Diagby2 = myPGDB[myPGdb::ECAL1_BAR].cell_size*sqrt(2.0)/2.0;
double barrel_rInner = myPGDB[myPGdb::ECAL1_BAR].r_inner; //inner ECAL radius
double barrel_zMax = myPGDB[myPGdb::ECAL1_CAP].z_inner; //maximum z barrel
double Xatf[3];
vector mySuperhits[Nlevels]; //lcio does not allow for more than 16 levels
for ( unsigned int hitCounter = 0; hitCountersize(); hitCounter++ ) {
int myLevel = 0;
//search for cluster level thresholds: energy per hit > mipstep[level]
for ( unsigned int levelCounter=0; levelCountermip > mipstep[levelCounter] )
myLevel = levelCounter;
else
break;
}
//write back hits into correct levels
for ( int levelCounter=myLevel; levelCounter>=0; levelCounter-- ) {
mySuperhits[levelCounter].push_back((*secal1)[hitCounter]);
if ( debug >3 )
std::cout<<"DEBUG KITutil: * writing hits of level "< testvector;
//assigning neighbouring hits from superHitVector to temporary cluster for each level
for ( unsigned int levelCounter=0; levelCountercalcCenter(); //energy weighted cluster center
if (debug >4 )
std::cout<<"calculated cluster center for level "<findInertia();//cluster start
if (debug >4 )
std::cout<<"calculated cluster start for level "< photonSeed;
for (unsigned int clusterCounter=0; clusterCounterhits.size()>Ccut.MinHit0 ) {
//yet another protseed
PROTSEED2 km;
km.cl = clusterVector[0][clusterCounter];
km.level = 0;
km.X[0] = 0.0;
km.X[1] = 0.0;
km.X[2] = 0.0;
km.active = true;
photonSeed.push_back(km);
}
}
for ( unsigned int seedCounter=0; seedCountergetCenter(), photonSeed[seedCounter].cl->direction,barrel_rInner,barrel_zMax,Xatf,debug);//cluster center, cluster direction, inner barrel radius, inner z barrel,
if (debug >4 ) {
std::cout<<"calculated cluster - IP intersect for seed "<getCenter()<direction<push_back(photonSeed[seedCounter]);
if (debug >5 )
std::cout<<"assign shifted hits to new photonSeed of size"<size()<size(); seed2Counter++) {
if ( (*photonSeed2)[seed2Counter].active ) {
Tmpclvec2 tempCluster; //clusters inside clusters
ClusterInCluster2((*photonSeed2)[seed2Counter].cl,clusterVector[levelCounter],tempCluster, debug);
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
//FIXME: CAREFULL photonSeed.Size increases 1 for each call of ClusterInCluster2 ! -> Infinite loop !
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
if (debug >4 )
std::cout<<"after ClusterInCluster2: photonSeed size = "<size()<4 )
std::cout<<"did not find cluster in cluster "<4 )
std::cout<<"level "< largerDistanceCluster;
vector < int > smallerDistanceCluster;
//check included clusters
for ( unsigned int clusterCounter1=0; clusterCounter1hits.size()>=Ccut.MinHitSplit //more than minimum numgber of hits in included cluster
&& tempCluster[clusterCounter2]->hits.size()>=Ccut.MinHitSplit ) {
double XIP[3];
XIP[0] = 0.0;
XIP[1] = 0.0;
XIP[2] = 0.0;
if ( (LinePointDistance2(XIP,(*photonSeed2)[seed2Counter].cl->center,tempCluster[clusterCounter2]->center) > Diagby2
|| LinePointDistance2(XIP,(*photonSeed2)[seed2Counter].cl->center,tempCluster[clusterCounter1]->center) > Diagby2 ) ) {//distance of tempCluster to photonSeed2 > cellSize *sqrt(2) / 2 ->clear cluster split
if ( find(largerDistanceCluster.begin(),largerDistanceCluster.end(),(int)clusterCounter2)==largerDistanceCluster.end() )
largerDistanceCluster.push_back(clusterCounter2);
if ( find(largerDistanceCluster.begin(),largerDistanceCluster.end(),(int)clusterCounter1)==largerDistanceCluster.end() )
largerDistanceCluster.push_back(clusterCounter1);
} else { //small distance between clusters -> split unclear
if (debug >5 )
std::cout<<"not so clearly seperated small distant cluster -> Split unclear"<5 )
std::cout<<"clearly seperated large distant cluster "<6 )
std::cout<<"for clearly seperated large distant clusters: photonSeed size = "<size()<hits.size();
for ( unsigned int vectorEntry=0; vectorEntryhits.size();
double ratio = ((double)clusterSize1) / ((double)clusterSize3);
(*photonSeed2)[seed2Counter].active = false;
if ( ratio > Ccut.rCut ) { //require minimum ratio of cluster entries to supress fluctuations
for ( unsigned int clusterEntry=0; clusterEntry5 )
std::cout<<"cluster ratio larger than"<getCenter(),tempCluster[smallerDistanceCluster[clusterEntry]]->direction,barrel_rInner,barrel_zMax,Xatf,debug);
clusterSeed.X[0] = Xatf[0];
clusterSeed.X[1] = Xatf[1];
clusterSeed.X[2] = Xatf[2];
clusterSeed.active = true;
bool split = true;
for ( unsigned int seed2Counter=0; seed2Countersize(); seed2Counter++ ) {
if ( (*photonSeed2)[seed2Counter].active ) {
if ( (*photonSeed2)[seed2Counter].cl==clusterSeed.cl )
split=false;
if (debug >5 )
std::cout<<"do not split cluster"<push_back(clusterSeed);
if (debug >5 )
std::cout<<"splitting cluster"< ratio cut
for ( unsigned int clusterEntry=0; clusterEntry4 )
std::cout<<"cluster larger than "<getCenter(),tempCluster[largerDistanceCluster[clusterEntry]]->direction,barrel_rInner,barrel_zMax,Xatf,debug);
clusterSeed.X[0] = Xatf[0];
clusterSeed.X[1] = Xatf[1];
clusterSeed.X[2] = Xatf[2];
clusterSeed.active = true;
bool split = true;
for ( unsigned int seed2Counter=0; seed2Countersize(); seed2Counter++ ) {
if ( (*photonSeed2)[seed2Counter].active ) {
if ( (*photonSeed2)[seed2Counter].cl==clusterSeed.cl )
split = false;
}
}
if ( split )
photonSeed2->push_back(clusterSeed);
}
}//found split clusterSeed at large distance
if ( debug >6 ) std::cout<<"exit large distance cluster"<6 ) std::cout<<"exit cluster in cluster"<6 ) std::cout<<"exit photonseed.active"<6 ) std::cout<<"for seed "<size()<6 ) std::cout<<"exit loop over photonseed"<6 ) std::cout<<"exit loop over threhold limits"< mergeVector;
for ( unsigned int hitCounter=0; hitCountersize(); hitCounter++ ) {
if( (*photonSeed2)[hitCounter].active ) {
double Xa[3];
ModuleNormal2((*photonSeed2)[hitCounter].cl->center,barrel_zMax,Xa); //shift hits in barrel
if (debug >3 )
std::cout<<" start merging: shifting barrel hits"<size(); hitCounter2++) {
//if hitcounters do not agree, get difference between these hits
if ( hitCounter2!=hitCounter && (*photonSeed2)[hitCounter2].active ) {
double dir_diff[3];
dir_diff[0] = (*photonSeed2)[hitCounter].cl->center[0]
- (*photonSeed2)[hitCounter2].cl->center[0];
dir_diff[1] = (*photonSeed2)[hitCounter].cl->center[1]
- (*photonSeed2)[hitCounter2].cl->center[1];
dir_diff[2] = (*photonSeed2)[hitCounter].cl->center[2]
- (*photonSeed2)[hitCounter2].cl->center[2];
//MERGE CONDITION for clusters:
//D_cl_cl2 returns smallest distance between cluster1 & cluster2
//Dot2 returns angle between clusters center & split candidate ??
if ( D_cl_cl2((*photonSeed2)[hitCounter].cl,(*photonSeed2)[hitCounter2].cl) < Ccut.distCut //minimum distance between clusters < distance cut
&& fabs( Dot2((*photonSeed2)[hitCounter2].cl->center,dir_diff)) > Ccut.cosCut) { //minimum angle between clusters
test cluster1,cluster2; //pair of cluster1 & cluster2
cluster1.first = hitCounter;
cluster1.second = hitCounter2;
cluster2.first = hitCounter2;
cluster2.second = hitCounter;
//write result in test vector
if ( find(mergeVector.begin(),mergeVector.end(),cluster1)==mergeVector.end()
&& find(mergeVector.begin(),mergeVector.end(),cluster2)==mergeVector.end() ) {
mergeVector.push_back(cluster1);
if (debug >3 )
std::cout<<" merging clusters"< differences in hits
if ( mergeVector.size()!=0) {
int benefit[mergeVector.size()];
for ( unsigned int vectorCounter=0; vectorCounter trinityVector; //holds entries from merged clusters
trinityVector.push_back(mergeVector[vectorCounter].first);
trinityVector.push_back(mergeVector[vectorCounter].second);
benefit[vectorCounter] = 1;
for(unsigned int vectorCounter2=0; vectorCounter2vectorCounter1 != vectorCounter
benefit[vectorCounter2] = 1;
if ( (*photonSeed2)[mergeVector[vectorCounter2].first].active )
trinityVector.push_back(mergeVector[vectorCounter2].first);
}
} else if ( benefit[vectorCounter2] != 1 ) {
if( (*photonSeed2)[mergeVector[vectorCounter2].second].active)
trinityVector.push_back(mergeVector[vectorCounter2].second);
benefit[vectorCounter2] = 1;
}
}
}
int levelCheck = 1;
for ( unsigned int vectorEntry=1; vectorEntryhits.size(); hitCounter++ ) {
cl->hits.push_back((*photonSeed2)[trinityVector[vectorEntry]].cl->hits[hitCounter]);
}
(*photonSeed2)[trinityVector[vectorEntry]].active = false;
}
}
if ( cl->hits.size()!=0 ) {//found merged cluster
cl->calcCenter();
if (debug >4 )
std::cout<<" calculate center of merged cluster"<findInertia();
if (debug >4 )
std::cout<<" calculate start of merged cluster"<getCenter(),cl->direction,barrel_rInner,barrel_zMax,Xatf,debug);
clusterSeed.X[0] = Xatf[0];
clusterSeed.X[1] = Xatf[1];
clusterSeed.X[2] = Xatf[2];
clusterSeed.active = true;
bool split = true;
for ( unsigned int seedEntry=0; seedEntrysize(); seedEntry++ ) {
if ( (*photonSeed2)[seedEntry].active ) {
if( (*photonSeed2)[seedEntry].cl==clusterSeed.cl)
split = false;
}
}
if(split)
photonSeed2->push_back(clusterSeed);
} else {//no merged clusters
delete cl;
}
} else {//levels do not agree
if (debug >5 )
std::cout<<" clusters to merge have different levels -> create new vector"< newVector;
for ( unsigned int j=0; jgetEnergy() >=(*photonSeed2)[trinityVector[1]].cl->getEnergy() ) {
(*photonSeed2)[trinityVector[1]].active = false;
} else {
(*photonSeed2)[trinityVector[0]].active = false;
}
} else {//trinityVector.size !=2
vector selectedHits;
int sumedUp[trinityVector.size()];
for ( unsigned int trinityEntry=0; trinityEntry sumUpVector;
for ( unsigned int trinityEntry=0; trinityEntry1 ) {
Tmpcl2* cl = new Tmpcl2();
for ( unsigned int sumUpVectorEntry=0; sumUpVectorEntryhits.size(); hitCounter++ )
cl->hits.push_back((*photonSeed2)[sumUpVector[sumUpVectorEntry]].cl->hits[hitCounter]);
if (debug >5 )
std::cout<<" writing new, merged cluster"<calcCenter();
cl->findInertia();
clusterVector[newVector[newVectorEntry]].push_back(cl);
Xatf[0] = 0.0;
Xatf[1] = 0.0;
Xatf[2] = 0.0;
PROTSEED2 clusterSeed;
clusterSeed.cl = cl;
clusterSeed.level = (*photonSeed2)[newVector[newVectorEntry]].level;
LineCaloIntersectD2(cl->getCenter(),cl->direction,barrel_rInner,barrel_zMax,Xatf,debug);
clusterSeed.X[0] = Xatf[0];
clusterSeed.X[1] = Xatf[1];
clusterSeed.X[2] = Xatf[2];
clusterSeed.active = true;
bool split = true;
for ( unsigned int im=0; imsize(); im++ ) {
if ( (*photonSeed2)[im].active ) {
if ( (*photonSeed2)[im].cl==clusterSeed.cl )
split=false;
}
}
if ( split ) {
photonSeed2->push_back(clusterSeed);
selectedHits.push_back(photonSeed2->size()-1);
} else {
delete cl;
}
} else {
selectedHits.push_back(sumUpVector[0]);
}
}//all newVector entries
if ( selectedHits.size()>1 ) {
if (debug >5 )
std::cout<<" pairing merged cluster"< DpairI;
vector selectedPairsVector;
for ( unsigned int j=0; jgetEnergy();
tmp.second = selectedHits[j];
selectedPairsVector.push_back(tmp);
}
sort(selectedPairsVector.begin(),selectedPairsVector.end());
for ( unsigned int j=0; j2 )
std::cout<<" ***** FindCores2 is done ***"<* superHitVector, vector* clusterVector, int debug )
{
if ( (*superHitVector).size()!=0 ) {
if ( debug > 5 )
std::cout<<"DEBUG KITutil: * clustering neighbouring hits of same level "<is_assigned = false;
//fill superHitVector into temporary clusterVector
if ( (*superHitVector)[hitCounter]->top==0 ) {
Tmpcl2* tempCluster = new Tmpcl2();
tempCluster->hits.push_back((*superHitVector)[hitCounter]);
(*superHitVector)[hitCounter]->is_assigned = true;
clusterVector->push_back(tempCluster);
continue;
} else if ( (*superHitVector)[hitCounter]->is_assigned==false ) {
vector sshv;
stack superhitStack;
sshv.push_back((*superHitVector)[hitCounter]);
//assigning neighbouring hits
for ( unsigned int neighbourCounter=0; neighbourCounter<(*superHitVector)[hitCounter]->neighbours.size(); neighbourCounter++ ) {
if ( !((*superHitVector)[hitCounter]->neighbours[neighbourCounter]->is_assigned)
&& superHitVector->end() != find(superHitVector->begin(),superHitVector->end(),(*superHitVector)[hitCounter]->neighbours[neighbourCounter]) ) {
sshv.push_back((*superHitVector)[hitCounter]->neighbours[neighbourCounter]);
(*superHitVector)[hitCounter]->neighbours[neighbourCounter]->is_assigned = true;
superhitStack.push((*superHitVector)[hitCounter]->neighbours[neighbourCounter]);
}
}//all neighbours
while ( !superhitStack.empty() ) {
Superhit2* superHit = superhitStack.top();
superhitStack.pop();
for ( unsigned int neighbourCounter=0; neighbourCounterneighbours.size(); neighbourCounter++ ) {
if ( superHit->neighbours[neighbourCounter]!=0 ) {
if ( ! superHit->neighbours[neighbourCounter]->is_assigned &&
superHitVector->end() != find(superHitVector->begin(),superHitVector->end(),superHit->neighbours[neighbourCounter]) )
if ( sshv.end() == find(sshv.begin(),sshv.end(),superHit->neighbours[neighbourCounter]) ) {
sshv.push_back(superHit->neighbours[neighbourCounter]);
superHit->neighbours[neighbourCounter]->is_assigned=true;
superhitStack.push(superHit->neighbours[neighbourCounter]);
}
} // neighbours not 0
} // all neighbours
if ( debug >5 )
std::cout<<"DEBUG KITutil: * found neighbours "<neighbours.size()<hits.push_back(sshv[sshvCounter]);
}
clusterVector->push_back(tempCluster);
}
} // no superhit assigned
} // over all hits
for ( unsigned int superHitVectorCounter=0; superHitVectorCountersize(); superHitVectorCounter++ )
if ( (*superHitVector)[superHitVectorCounter]->is_assigned == false ) {
Tmpcl2* tempCluster = new Tmpcl2();
tempCluster->hits.push_back((*superHitVector)[superHitVectorCounter]);
clusterVector->push_back(tempCluster);
if ( debug >5 )
std::cout<<"DEBUG KITutil: * found cluster "<getPosition()[positionCounter];
}
mip = chit->getEnergy() / E; //mip calibration
connect = false;
is_assigned = false;
mipE = E;
top = 0;
cl = 0;
S = 0; //stave
M = 0; //module
K = 0; //layer
isCalorimeter = 0; // ecal=1, hcal=2, else=0
}
mySuperhit::~mySuperhit () {}
mySuperhit::mySuperhit ( float E, CalorimeterHit* chit_in ) {
chit = chit_in;
for ( unsigned int positionCounter=0; positionCounter<3; positionCounter++ ) {
point[positionCounter] = 0.0;
pointt[positionCounter] = chit_in->getPosition()[positionCounter];
}
mip = chit->getEnergy() / E; //mip calibration
connect = false;
is_assigned = false;
mipE = E;
top = 0;
cl = 0;
S = 0; //stave
M = 0; //module
K = 0; //layer
}
double* myTmpcl:: getCenter()
{
return center;
}
myTmpcl::myTmpcl(){
energy=0.0;
}
myTmpcl::~myTmpcl(){}
double myTmpcl::getEnergy(){
energy =0.0;
for ( unsigned int i=0; ichit->getEnergy();
}
return energy;
}
//calculating energy weighted center of cluster
void myTmpcl::calcCenter()
{
center[0] = 0.0;
center[1] = 0.0;
center[2] = 0.0;
double energySum = 0.0;
double energy = 0.0;
for ( unsigned int hitCounter=0; hitCounterchit->getEnergy();
center[0] += hits[hitCounter]->pointt[0] * energy;
center[1] += hits[hitCounter]->pointt[1] * energy;
center[2] += hits[hitCounter]->pointt[2] * energy;
energySum += energy;
}
center[0] = center[0] / energySum;
center[1] = center[1] / energySum;
center[2] = center[2] / energySum;
}
//finding cluster start
void myTmpcl::findInertia()
{
double aIne[3][3];
for ( int xPosCounter(0); xPosCounter<3; ++xPosCounter )
for ( int yPosCounter(0); yPosCounter < 3; ++yPosCounter ) {
aIne[xPosCounter][yPosCounter] = 0.0;
}
for ( unsigned int hitCounter=0; hitCounterpointt[0] - center[0];
double distY = hits[hitCounter]->pointt[1] - center[1];
double distZ = hits[hitCounter]->pointt[2] - center[2];
double eHit = hits[hitCounter]->chit->getEnergy();
eHit = 1.0;
aIne[0][0] += eHit * (distY*distY + distZ*distZ);
aIne[1][1] += eHit * (distX*distX + distZ*distZ);
aIne[2][2] += eHit * (distX*distX + distY*distY);
aIne[0][1] -= eHit * distX*distY;
aIne[0][2] -= eHit * distX*distZ;
aIne[1][2] -= eHit * distY*distZ;
}
for (int i(0); i<2; ++i ) {
for (int j=i+1; j<3; ++j ) {
aIne[j][i] = aIne[i][j];
}
}
gsl_matrix_view aMatrix = gsl_matrix_view_array((double*)aIne,3,3);
gsl_vector* aVector = gsl_vector_alloc(3);
gsl_matrix* aEigenVec = gsl_matrix_alloc(3,3);
gsl_eigen_symmv_workspace* wa = gsl_eigen_symmv_alloc(3);
gsl_eigen_symmv(&aMatrix.matrix,aVector,aEigenVec,wa);
gsl_eigen_symmv_free(wa);
gsl_eigen_symmv_sort(aVector,aEigenVec,GSL_EIGEN_SORT_ABS_ASC);
double norm = 0.0;
for ( int i(0); i < 3; i++ ) {
inteigen[i] = gsl_vector_get(aVector,i);
norm += inteigen[i] * inteigen[i];
for (int j(0); j<3; j++ ) {
inteigenvec[i+3*j] = gsl_matrix_get(aEigenVec,i,j);
}
}
norm = sqrt(norm);
for ( int i(0); i<3; i++ )
inteigen[i] = inteigen[i]/norm;
direction[0] = inteigenvec[0];
direction[1] = inteigenvec[1];
direction[2] = inteigenvec[2];
gsl_vector_free(aVector);
gsl_matrix_free(aEigenVec);
}//find inertia
double* Tmpcl2:: getCenter()
{
return center;
}
Tmpcl2::Tmpcl2(){
energy=0.0;
}
Tmpcl2::~Tmpcl2(){}
double Tmpcl2::getEnergy(){
energy =0.0;
for ( unsigned int i=0; ichit->getEnergy();
}
return energy;
}
//calculating energy weighted center of cluster
void Tmpcl2::calcCenter()
{
center[0] = 0.0;
center[1] = 0.0;
center[2] = 0.0;
double energySum = 0.0;
double energy = 0.0;
for ( unsigned int hitCounter=0; hitCounterchit->getEnergy();
center[0] += hits[hitCounter]->pointt[0] * energy;
center[1] += hits[hitCounter]->pointt[1] * energy;
center[2] += hits[hitCounter]->pointt[2] * energy;
energySum += energy;
}
center[0] = center[0] / energySum;
center[1] = center[1] / energySum;
center[2] = center[2] / energySum;
}
//finding cluster start
void Tmpcl2::findInertia()
{
double aIne[3][3];
for ( int xPosCounter(0); xPosCounter<3; ++xPosCounter )
for ( int yPosCounter(0); yPosCounter < 3; ++yPosCounter ) {
aIne[xPosCounter][yPosCounter] = 0.0;
}
for ( unsigned int hitCounter=0; hitCounterpointt[0] - center[0];
double distY = hits[hitCounter]->pointt[1] - center[1];
double distZ = hits[hitCounter]->pointt[2] - center[2];
double eHit = hits[hitCounter]->chit->getEnergy();
eHit = 1.0;
aIne[0][0] += eHit * (distY*distY + distZ*distZ);
aIne[1][1] += eHit * (distX*distX + distZ*distZ);
aIne[2][2] += eHit * (distX*distX + distY*distY);
aIne[0][1] -= eHit * distX*distY;
aIne[0][2] -= eHit * distX*distZ;
aIne[1][2] -= eHit * distY*distZ;
}
for (int i(0); i<2; ++i ) {
for (int j=i+1; j<3; ++j ) {
aIne[j][i] = aIne[i][j];
}
}
gsl_matrix_view aMatrix = gsl_matrix_view_array((double*)aIne,3,3);
gsl_vector* aVector = gsl_vector_alloc(3);
gsl_matrix* aEigenVec = gsl_matrix_alloc(3,3);
gsl_eigen_symmv_workspace* wa = gsl_eigen_symmv_alloc(3);
gsl_eigen_symmv(&aMatrix.matrix,aVector,aEigenVec,wa);
gsl_eigen_symmv_free(wa);
gsl_eigen_symmv_sort(aVector,aEigenVec,GSL_EIGEN_SORT_ABS_ASC);
double norm = 0.0;
for ( int i(0); i < 3; i++ ) {
inteigen[i] = gsl_vector_get(aVector,i);
norm += inteigen[i] * inteigen[i];
for (int j(0); j<3; j++ ) {
inteigenvec[i+3*j] = gsl_matrix_get(aEigenVec,i,j);
}
}
norm = sqrt(norm);
for ( int i(0); i<3; i++ )
inteigen[i] = inteigen[i]/norm;
direction[0] = inteigenvec[0];
direction[1] = inteigenvec[1];
direction[2] = inteigenvec[2];
gsl_vector_free(aVector);
gsl_matrix_free(aEigenVec);
}//find inertia
//find shower start in first ECAL layer?
//where does 400 come from?
void LineCaloIntersectD2( double* X1,double* dir,double& r_inner,double& zmax,double* X,int debug)
{ //X1: cluster center; dir: cluster direction; r_inner: inner barrel radius, zmax: maximum barrel z-value; X: 0
if(debug >4)
std::cout<<"---- FindCores2 called LineCaloIntersect -----"<4)
std::cout<<"---- Finished LineCaloIntersect ----"< assumes 8fold detector geometry + 2 endcaps!
//searches for cluster closest to IP
void LineCaloIntersect2(double* X1, double* X2,double& r_inner,double& zmax,double* X, int debug)
{
double n[10][3];
double sqrt2_half = M_SQRT2/2.0; //= sqrt(2)/2 = 0.71
//n[stave][xyz] ->stave geometry
n[0][0] = 0.0;
n[0][1] = 1.0;
n[0][2] = 0.0;
n[1][0] = sqrt2_half;
n[1][1] = sqrt2_half;
n[1][2] = 0.0;
n[2][0] = 1.0;
n[2][1] = 0.0;
n[2][2] = 0.0;
n[3][0] = -1.0;
n[3][1] = 0.0;
n[3][2] = 0.0;
n[4][0] = 0.0;
n[4][1] = -1.0;
n[4][2] = 0.0;
n[5][0] = -sqrt2_half;
n[5][1] = -sqrt2_half;
n[5][2] = 0.0;
n[6][0] = sqrt2_half;
n[6][1] = -sqrt2_half;
n[6][2] = 0.0;
n[7][0] = -sqrt2_half;
n[7][1] = sqrt2_half;
n[7][2] = 0.0;
//endcaps
n[8][0] = 0.0;
n[8][1] = 0.0;
n[8][2] = 1.0;
n[9][0] = 0.0;
n[9][1] = 0.0;
n[9][2] = -1.0;
if ( abs(X2[2])5)
std::cout<<"---- LineCaloIntersect2 found hit in barrel ----"< > test;
for ( unsigned int staveCounter=0; staveCounter<8; staveCounter++ ) {
double tmp = - (n[staveCounter][0]*n[staveCounter][0] + n[staveCounter][1]*n[staveCounter][1] + n[staveCounter][2]*n[staveCounter][2]) * r_inner
+ n[staveCounter][0]*X1[0] + n[staveCounter][1]*X1[1] + n[staveCounter][2]*X1[2];
double tmp2 = n[staveCounter][0]*(X1[0]-X2[0])
+ n[staveCounter][1]*(X1[1]-X2[1])
+ n[staveCounter][2]*(X1[2]-X2[2]);
double t = -2.0;
if ( tmp2!=0.0 ) {
t = tmp/tmp2;
//for stave[0]: t = (-r_inner+clusterCenterY)/(400*clusterDirectionY)
}
if (t>0.0 && t<=1.0 ) {
//for stave[0]: if (clusterCenterY>r_inner && -r_inner+clusterCenterY<400*clusterDirectionY) ->hit inside ECAL
vector tmpi;
tmpi.push_back(X1[0]+(X2[0]-X1[0])*t); //2*clusterCenterX-r_inner ->shifts cluster center 2*depper into the ECAL stave
tmpi.push_back(X1[1]+(X2[1]-X1[1])*t);
tmpi.push_back(X1[2]+(X2[2]-X1[2])*t);
test.push_back(tmpi);
}
}//all staves
if (test.size()==1 ) {
X[0] = test[0][0];
X[1] = test[0][1];
X[2] = test[0][2];
double rts = sqrt(X[0]*X[0] + X[1]*X[1]);
if ( rts< ((r_inner+2.0)*1.0823922002924) ) { //why +2 ?? why *1.08...??
return; //check wether hit still in ECAL??
} else {
X[0] = 0.0;
X[1] = 0.0;
X[2] = 0.0;
return;//if so, seteverything back
}
} //if test positive for exactly one stave
else {
double rmin = 1e+22;
unsigned int imin = 999;
//loop over all positive staves
for ( unsigned int i=0; i5)
std::cout<<"---- LineCaloIntersect2 found hit in endcap ----"<0.0 ) {
n[0][0] = 0.0;
n[0][1] = 0.0;
n[0][2] = 1.0;
} else {
n[0][0] = 0.0;
n[0][1] = 0.0;
n[0][2] = -1.0;
}
double tmp = -(n[0][0]*n[0][0] + n[0][1]*n[0][1] + n[0][2]*n[0][2]) * zmax
+ n[0][0]*X1[0] + n[0][1]*X1[1] + n[0][2]*X1[2];
double tmp2 = n[0][0]*(X1[0]-X2[0])
+ n[0][1]*(X1[1]-X2[1])
+ n[0][2]*(X1[2]-X2[2]);
double t = 0.0;
if ( tmp2!=0.0 ) t=tmp/tmp2;
if ( t>0.0 && t<=1.0 ) {
X[0] = X1[0]+(X2[0]-X1[0])*t ;
X[1] = X1[1]+(X2[1]-X1[1])*t ;
X[2] = X1[2]+(X2[2]-X1[2])*t ;
return;
} else {
X[0] = 0.0;
X[1] = 0.0;
X[2] = 0.0;
return;
}
}//endcap end
}
void ClusterInCluster(myTmpcl* cl, vector& clv,vector& clout)
{
sort(cl->hits.begin(),cl->hits.end());//sort hits in first cluster
for (unsigned int i=0; ihits.begin(),clv[i]->hits.end());//sort hits in 2nd cluster
if ( includes(cl->hits.begin(),cl->hits.end(),clv[i]->hits.begin(),clv[i]->hits.end()) )
//every element of clusterVector[i] included in cl
clout.push_back(clv[i]);//return included cluster
}
// std::cout<<"ClusterInCluster loop over "< filled tempCluster"<& clusterVector, int debug)
{
sort(cl->hits.begin(),cl->hits.end());
for ( unsigned int vectorCounter=0; vectorCounterhits.begin(),clusterVector[vectorCounter]->hits.end());
if ( includes(cl->hits.begin(),cl->hits.end(),clusterVector[vectorCounter]->hits.begin(),clusterVector[vectorCounter]->hits.end()) ) {
cl->daughters.push_back(clusterVector[vectorCounter]);
clusterVector[vectorCounter]->parents.push_back(cl);
}
}
}
//called by FindCores2
//returns distance of included cluster hit to surrounding cluster hit ??
double LinePointDistance2( double* X1, double* X2, double* X0){
double tmp1 = X1[1]*X2[0] - X1[0]*X2[1] - X1[1]*X0[0] + X2[1]*X0[0] + X1[0]*X0[1] - X2[0]*X0[1];
double tmp2 = X1[2]*X2[0] - X1[0]*X2[2] - X1[2]*X0[0] + X2[2]*X0[0] + X1[0]*X0[2] - X2[0]*X0[2];
double tmp3 = X1[2]*X2[1] - X1[1]*X2[2] - X1[2]*X0[1] + X2[2]*X0[1] + X1[1]*X0[2] - X2[1]*X0[2];
double tmp4 = (X1[0]-X2[0])*(X1[0]-X2[0]) + (X1[1]-X2[1])*(X1[1]-X2[1]) + (X1[2]-X2[2])*(X1[2]-X2[2]);
double tmp5 = tmp1*tmp1 + tmp2*tmp2 + tmp3*tmp3;
tmp5 = tmp5/tmp4;
return sqrt(tmp5);
}
//called by findCores2
void ModuleNormal2(double* X1,double& zmax, double* X0)
{
if ( abs(X1[2]) very similar to precalculation
double phi = atan2(X1[0],X1[1]);
double phi0 = M_PI_4/2.0; //pi/4/2 = 0.39
int i;
if ( phi>=0. )
i = (int)((phi+phi0) / M_PI_4);
else
i = (int)((phi-phi0) / M_PI_4);
double phi_p = i * M_PI_4;
if ( i!=0 ) {
X0[0] = sin(phi_p);
X0[1] = cos(phi_p);
X0[2] = 0.0;
} else {
X0[0] = 0.0;
X0[1] = 1.0;
X0[2] = 0.0;
}
} else {
if ( X1[2]>0. ) {
X0[0] = 0.0;
X0[1] = 0.0;
X0[2] = 1.0;
} else {
X0[0] = 0.0;
X0[1] = 0.0;
X0[2] = -1.0;
}
}
}
double D_cl_cl(myTmpcl* cl1,myTmpcl* cl2)
{
vector dist;
for ( unsigned int cluster1entry=0; cluster1entryhits.size(); cluster1entry++ ) {
double x1 = cl1->hits[cluster1entry]->pointt[0];
double y1 = cl1->hits[cluster1entry]->pointt[1];
double z1 = cl1->hits[cluster1entry]->pointt[2];
for ( unsigned int cluster2entry=0; cluster2entryhits.size(); cluster2entry++ ) {
double tmp = ( (x1 - cl2->hits[cluster2entry]->pointt[0]) * (x1 - cl2->hits[cluster2entry]->pointt[0])
+ (y1 - cl2->hits[cluster2entry]->pointt[1]) * (y1 - cl2->hits[cluster2entry]->pointt[1])
+ (z1 - cl2->hits[cluster2entry]->pointt[2]) * (z1 - cl2->hits[cluster2entry]->pointt[2]) ) ;
dist.push_back(tmp);
}
}
sort(dist.begin(),dist.end());
return sqrt(dist[0]); //return smallest distance between cluster1 & cluster2
}
//calculate distance between 2 clusters called by FindCores
double D_cl_cl2(Tmpcl2* cl1,Tmpcl2* cl2)
{
vector dist;
for ( unsigned int cluster1entry=0; cluster1entryhits.size(); cluster1entry++ ) {
double x1 = cl1->hits[cluster1entry]->pointt[0];
double y1 = cl1->hits[cluster1entry]->pointt[1];
double z1 = cl1->hits[cluster1entry]->pointt[2];
for ( unsigned int cluster2entry=0; cluster2entryhits.size(); cluster2entry++ ) {
double tmp = ( (x1 - cl2->hits[cluster2entry]->pointt[0]) * (x1 - cl2->hits[cluster2entry]->pointt[0])
+ (y1 - cl2->hits[cluster2entry]->pointt[1]) * (y1 - cl2->hits[cluster2entry]->pointt[1])
+ (z1 - cl2->hits[cluster2entry]->pointt[2]) * (z1 - cl2->hits[cluster2entry]->pointt[2]) ) ;
dist.push_back(tmp);
}
}
sort(dist.begin(),dist.end());
return sqrt(dist[0]); //return smallest distance between cluster1 & cluster2
}
// called by FindCores
inline double Dot2(double* X1,double* X2)
{
double n1 = sqrt(X1[0]*X1[0] + X1[1]*X1[1] + X1[2]*X1[2]);
double n2 = sqrt(X2[0]*X2[0] + X2[1]*X2[1] + X2[2]*X2[2]);
return (X1[0]*X2[0] + X1[1]*X2[1] + X1[2]*X2[2]) / (n1*n2);
}
void ClusterInCluster2(Tmpcl2* cl, vector& clusterVector, vector& clout, int debug)
{
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
//FIXME: The main problem seems to be this sort routine, but why?
//This should not touch the content, only the order of the cluster!
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
//sort(cl->hits.begin(),cl->hits.end());//sort hits in first cluster
for (unsigned int i=0; ihits.begin(),clusterVector[i]->hits.end());//sort hits in 2nd cluster
if ( includes(cl->hits.begin(),cl->hits.end(),clusterVector[i]->hits.begin(),clusterVector[i]->hits.end()) )
//every element of clusterVector[i] included in cl
clout.push_back(clusterVector[i]);//return included cluster
}
if ( debug >2)
std::cout<<"ClusterInCluster loop over "< filled tempCluster"<getPosition(),X);
double TP[3];
TP[0] = ch->getPosition()[0] - start[0];
TP[1] = ch->getPosition()[1] - start[1];
TP[2] = ch->getPosition()[2] - start[2];
// distance must have direction
if ( Dot2(TP,dir)<0.0 ) {
out[0] = 0.0;
out[1] = 0.0;
return ;
}
double pos_X0 = sqrt( (start[0]-X[0])*(start[0]-X[0])+
(start[1]-X[1])*(start[1]-X[1])+
(start[2]-X[2])*(start[2]-X[2]))/x0eff;
double tau = pos_X0/Tsam;
double pos[3];
pos[0] = ch->getPosition()[0];
pos[1] = ch->getPosition()[1];
pos[2] = ch->getPosition()[2];
double radius = LinePointDistance2(start,dir,pos);
// average radial profiles
double RChom = z1 + z2*tau;
double RThom = k1 * (exp(k3*(tau-k2)) + exp(k4*(tau-k2)));
double phom = p1 * exp((p2-tau)/p3 - exp((p2-tau)/p3));
// average radial profiles sampling
double RCsam = RChom - 0.0203*(1.0-eprime) + 0.0397*exp(-tau)/Fs;
double RTsam = RThom - 0.14*(1.0-eprime) - 0.495*exp(-tau)/Fs;
double psam = phom + (1.0-eprime)*(0.348-0.642*exp(-pow(tau-1.0,2.0))/Fs);
double rb = radius/Rm;
if ( rb<20.0 && pos_X0<35.0 ) {
double tmp = pow(betasam*pos_X0,alfasam-1.0)*betasam*exp(-betasam*pos_X0)/gsl_sf_gamma(alfasam);
tmp = Ee*tmp*(psam*2.0*rb*RCsam*RCsam/pow(rb*rb+RCsam*RCsam,2.0) +
(1-psam)*(2.0*rb*RTsam*RTsam/pow(rb*rb+RTsam*RTsam,2.0)))/rb;
if ( tmp>cut ) {
out[0] = tmp;
out[1] = radius;
// std::cout<<"set probability to "< set to 0."< set to 0."< cc)
{
for ( unsigned int i=0; i "<Ecore ) { //FIXME: does find more than one appropriate calibration!
/*
std::cout << "For level " << level
<< " core with energy: " << Ecore
<< " : estimated photon energy = cc[i].b: " << cc[i].b
<< " + cc[i].a: " << cc[i].a
<< " * Ecore = " << cc[i].b + cc[i].a*Ecore << std::endl;
*/
if ( cc[i].b==0.0 || cc[i].a==0.0 )
return Ecore; //if no valid calibration, return uncalibrated value
else
return cc[i].b + cc[i].a*Ecore;
}
else
return Ecore; //if no valid calibration, return uncalibrated value
}
return 0.0;
}
// !!! has to be adopted for different detector models !!!
// relates measured energy in cluster to estimated true energy of the cluster,
// depending on the cluster level
// numbers = energy cut in the enery per hit spectrum per level, for various energies
// E_true = aa*E_cluster + b; E_min < E_cluster < E_max
// can be obtained with KALIBRATOR.cc
void CreateCalibration_LDC00(vector* cc)
{
double aEnom[9] = {0.5,1.0,2.0,3.0,5.0,8.0,12.0,20.0,40.0};
double aa[9][10] = {{0.283,0.357,0.337,0.398,0.000,0.000,0.000,0.000,0.000,0.000},
{0.703,0.346,0.814,0.785,0.939,0.000,0.000,0.000,0.000,0.000},
{0.505,0.821,1.000,1.112,1.527,1.860,0.000,0.000,0.000,0.000},
{0.809,0.892,1.999,1.554,2.170,2.870,0.000,0.000,0.000,0.000},
{1.075,0.827,1.290,2.350,2.123,3.958,4.814,0.000,0.000,0.000},
{1.140,1.310,1.765,3.652,4.030,5.600,6.670,0.000,0.000,0.000},
{0.760,1.845,3.444,2.850,5.150,6.760,10.72,0.000,0.000,0.000},
{3.770,4.460,2.270,5.930,8.730,11.00,17.40,17.77,0.000,0.000},
{1.890,3.560,5.120,7.320,9.650,11.71,21.74,0.000,0.000,0.000}};
double bb[9][10] = {{0.730,0.595,0.772,0.683,0.000,0.000,0.000,0.000,0.000,0.000},
{0.493,1.079,0.458,0.592,0.447,0.000,0.000,0.000,0.000,0.000},
{0.941,0.813,0.766,0.782,0.559,0.405,0.000,0.000,0.000,0.000},
{0.892,0.910,0.550,0.809,0.593,0.347,0.000,0.000,0.000,0.000},
{0.918,1.015,0.967,0.800,0.954,0.510,0.356,0.000,0.000,0.000},
{0.974,0.986,0.970,0.778,0.784,0.593,0.509,0.000,0.000,0.000},
{1.032,0.971,0.869,0.986,0.835,0.742,0.339,0.000,0.000,0.000},
{0.891,0.888,1.025,0.884,0.782,0.620,0.314,0.361,0.000,0.000},
{1.020,0.992,0.977,0.962,0.945,0.943,0.724,0.000,0.000,0.000}};
double aEmin[9][10] = {{0.250,0.200,0.150,0.060,0.000,0.000,0.000,0.000,0.000,0.000},
{0.600,0.500,0.350,0.200,0.100,0.000,0.000,0.000,0.000,0.000},
{1.500,1.200,1.100,0.900,0.650,0.300,0.000,0.000,0.000,0.000},
{2.400,2.100,1.900,1.600,1.300,0.600,0.000,0.000,0.000,0.000},
{4.000,4.100,3.900,3.200,2.900,2.400,1.100,0.000,0.000,0.000},
{7.000,6.800,6.200,5.500,5.000,4.200,2.600,0.000,0.000,0.000},
{10.50,10.00,9.500,8.600,8.400,8.300,5.200,0.000,0.000,0.000},
{18.50,18.20,17.40,16.00,15.50,13.00,11.50,9.000,0.000,0.000},
{36.00,35.00,34.00,33.50,32.50,31.00,28.00,0.000,0.000,0.000}};
double aEmax[9][10]={{0.600,0.500,0.350,0.300,0.000,0.000,0.000,0.000,0.000,0.000},
{1.200,1.100,1.000,0.900,0.700,0.000,0.000,0.000,0.000,0.000},
{2.500,2.200,1.900,1.800,1.400,1.200,0.000,0.000,0.000,0.000},
{3.300,3.100,3.000,2.700,2.400,2.000,0.000,0.000,0.000,0.000},
{5.600,5.500,5.200,4.800,4.300,3.900,2.900,0.000,0.000,0.000},
{9.000,8.800,8.300,7.600,7.200,6.400,5.200,0.000,0.000,0.000},
{13.50,13.00,12.50,11.70,11.10,11.00,8.300,0.000,0.000,0.000},
{22.20,22.80,21.50,20.50,19.00,16.30,15.50,13.50,0.000,0.000},
{44.00,43.50,42.50,41.50,39.20,37.50,33.50,0.000,0.000,0.000}};
for ( unsigned int levelCounter=0; levelCounter<10; levelCounter++ ) {
for ( unsigned int energyCounter=0; energyCounter<12; energyCounter++ ) {
if ( aEmin[energyCounter][levelCounter]!=0.0 && aEmax[energyCounter][levelCounter]!=0.0 ) {
CoreCalib2 cck;
cck.level = levelCounter;
cck.Enom = aEnom[energyCounter];
cck.a = aa[energyCounter][levelCounter];
cck.b = bb[energyCounter][levelCounter];
cck.Emin = aEmin[energyCounter][levelCounter];
cck.Emax = aEmax[energyCounter][levelCounter];
cc->push_back(cck);
}
}
}
}
void CreateCalibration_LDCPrime02Sc(vector* cc)
{
// fit straight line: E_true = aa*E_reco + bb for each energy and level
double aEnom[13] = {0.5,1.0,2.0,3.0,5.0,8.0,12.0,20.0,40.0,70.0,100.0,150.0,200.0}; //nominal photon energy
double aa[13][10] = {{0.390,0.386,0.875,0.327,0.000,0.000,0.000,0.000,0.000,0.000},
{0.574,0.437,0.365,0.380,0.378,0.000,0.000,0.000,0.000,0.000},
{0.693,0.644,0.735,0.701,0.663,0.386,0.000,0.000,0.000,0.000},
{0.992,1.003,0.703,0.920,0.705,0.666,0.000,0.000,0.000,0.000},
{0.942,0.894,0.857,0.685,0.851,0.719,0.000,0.000,0.000,0.000},
{0.850,0.876,1.036,0.806,0.746,0.532,0.000,0.000,0.000,0.000},
{0.911,0.901,0.851,0.847,0.796,0.693,0.468,0.464,0.000,0.000},
{0.832,1.005,1.006,0.920,0.975,0.865,0.655,0.441,0.472,0.000},
{0.969,0.899,0.854,0.939,0.858,0.762,0.716,0.378,0.417,0.000},
{1.056,1.073,1.046,0.894,0.879,1.127,1.140,1.199,0.308,0.650},
{0.738,1.080,2.044,1.089,1.032,0.862,0.885,0.758,1.294,0.422},
{0.879,0.879,0.879,0.879,0.879,0.879,0.879,0.879,0.000,0.000},
{0.836,0.836,0.836,0.836,0.836,0.836,0.386,0.836,0.000,0.000}};
double bb[13][10] = {{0.412,0.436,0.383,0.471,0.000,0.000,0.000,0.000,0.000,0.000},
{0.690,0.784,0.870,0.886,0.926,0.000,0.000,0.000,0.000,0.000},
{1.095,1.209,1.155,1.281,1.422,1.748,0.000,0.000,0.000,0.000},
{0.894,0.965,1.675,1.454,1.986,2.251,0.000,0.000,0.000,0.000},
{1.466,1.804,2.109,2.900,2.760,3.442,0.000,0.000,0.000,0.000},
{2.658,2.722,2.135,3.841,4.589,5.984,0.000,0.000,0.000,0.000},
{3.158,3.588,4.480,5.106,6.226,7.719,9.918,10.47,0.000,0.000},
{6.274,3.616,4.306,6.645,7.095,9.971,14.18,17.14,17.78,0.000},
{5.938,9.247,11.98,10.65,14.84,20.11,25.15,33.79,34.32,0.000},
{3.345,3.711,7.038,18.97,22.79,13.86,21.53,-6.09,62.73,54.05},
{33.77,3.409,-81.3,8.968,17.70,36.57,43.36,60.02,35.92,85.20},
{29.30,29.30,29.30,29.30,29.30,29.30,29.30,29.30,0.000,0.000},
{46.13,46.13,46.13,46.13,46.13,46.13,46.13,46.13,0.000,0.000}};
double aEmin[13][10] = {{0.250,0.200,0.150,0.060,0.000,0.000,0.000,0.000,0.000,0.000},
{0.600,0.500,0.350,0.200,0.100,0.000,0.000,0.000,0.000,0.000},
{1.500,1.200,1.100,0.900,0.650,0.300,0.000,0.000,0.000,0.000},
{2.400,2.100,1.900,1.600,1.300,0.600,0.000,0.000,0.000,0.000},
{4.000,4.100,3.900,3.200,2.900,2.400,0.000,0.000,0.000,0.000},
{4.279,3.848,3.501,4.618,2.439,1.695,0.000,0.000,0.000,0.000},
{9.189,6.657,6.144,5.337,4.811,3.780,2.749,0.954,0.000,0.000},
{14.29,13.64,12.69,11.56,10.76,8.683,6.213,3.714,3.040,0.000},
{31.50,30.82,28.61,28.23,25.75,22.29,18.14,13.61,9.857,0.000},
{57.10,56.41,51.97,51.54,48.47,32.96,39.48,28.08,25.26,23.19},
{62.30,84.58,71.61,68.98,59.43,62.63,53.12,33.60,35.83,34.79},
{134.9,134.9,134.9,194.9,134.9,134.9,134.9,134.9,0.000,0.000},
{182.7,182.7,182.7,182.7,182.7,165.0,135.0,112.0,0.000,0.000}};
double aEmax[13][10]={{0.600,0.500,0.350,0.300,0.000,0.000,0.000,0.000,0.000,0.000},
{1.200,1.100,1.000,0.900,0.700,0.000,0.000,0.000,0.000,0.000},
{2.500,2.200,1.900,1.800,1.400,1.200,0.000,0.000,0.000,0.000},
{3.300,3.100,3.000,2.700,2.400,2.000,0.000,0.000,0.000,0.000},
{5.600,5.500,5.200,4.800,4.300,3.900,0.000,0.000,0.000,0.000},
{8.464,8.442,8.051,6.774,6.945,6.120,0.000,0.000,0.000,0.000},
{12.79,12.45,12.01,11.35,10.21,9.137,7.234,6.474,0.000,0.000},
{20.40,19.61,19.38,18.38,16.55,15.58,12.83,11.91,9.463,0.000},
{40.94,39.63,39.58,36.69,34.23,32.34,27.23,24.58,22.26,0.000},
{71.53,70.61,73.92,66.78,63.47,71.53,49.70,50.04,42.17,34.28},
{136.1,100.0,106.0,139.5,140.0,92.21,85.62,91.04,72.17,132.5},
{152.6,152.6,152.6,152.6,152.6,152.6,152.6,152.6,0.000,0.000},
{204.1,204.1,204.1,204.1,204.1,204.1,204.1,204.1,0.000,0.000}};
for ( unsigned int levelCounter=0; levelCounter<10; levelCounter++ ) {
for ( unsigned int energyCounter=0; energyCounter<12; energyCounter++ ) {
if ( aEmin[energyCounter][levelCounter]!=0.0 && aEmax[energyCounter][levelCounter]!=0.0 ) {
CoreCalib2 cck;
cck.level = levelCounter;
cck.Enom = aEnom[energyCounter];
cck.a = aa[energyCounter][levelCounter];
cck.b = bb[energyCounter][levelCounter];
cck.Emin = aEmin[energyCounter][levelCounter];
cck.Emax = aEmax[energyCounter][levelCounter];
cc->push_back(cck);
}
}
}
}//creatCalibration_LDCPrime02Sc