| 1 |
pam-fi |
1.1 |
#include <CaloHough.h> |
| 2 |
|
|
//////////////////////////////////////////////////////////// |
| 3 |
|
|
/// INITIALIZATIONS & GLOBAL PARAMETERS |
| 4 |
|
|
//////////////////////////////////////////////////////////// |
| 5 |
|
|
|
| 6 |
|
|
CaloHough_parameters * CaloHough_parameters::_parameters = 0; |
| 7 |
|
|
|
| 8 |
|
|
/** |
| 9 |
|
|
* Set external parameters to default values |
| 10 |
|
|
*/ |
| 11 |
|
|
void CaloHough_parameters::SetDefault(){ |
| 12 |
|
|
|
| 13 |
|
|
cout << "CaloHough --> SET DEFAULT PARAMETERS" << endl; |
| 14 |
|
|
|
| 15 |
|
|
isimu = 0; |
| 16 |
|
|
|
| 17 |
|
|
// Q_UP = 4.;// mip |
| 18 |
|
|
// Q_DN = 0.; |
| 19 |
|
|
// M_UP = 3; |
| 20 |
|
|
// M_DN = 0; |
| 21 |
|
|
Q_UP = 10000.; |
| 22 |
|
|
Q_DN = 100.; |
| 23 |
|
|
M_UP = 3000; |
| 24 |
|
|
M_DN = 0; |
| 25 |
|
|
|
| 26 |
|
|
h_qtot[0] = new TH2F("houghqtot_0","Hit pattern",96,-0.5,95.5,22,-0.5,21.5); |
| 27 |
|
|
h_qtot[1] = new TH2F("houghqtot_1","Hit pattern",96,-0.5,95.5,22,-0.5,21.5); |
| 28 |
|
|
|
| 29 |
|
|
int n1 = 150; |
| 30 |
|
|
int n2 = 150; |
| 31 |
|
|
|
| 32 |
|
|
h_par[0] = new TH2F("houghpar_0","par1 vs par2",n1,-40.,40.,n2,-10.,10.); |
| 33 |
|
|
h_par[1] = new TH2F("houghpar_1","par1 vs par2",n1,-40.,40.,n2,-10.,10.); |
| 34 |
|
|
|
| 35 |
|
|
h_par_fine[0] = new TH2F("houghparfine_0","par1 vs par2",10*n1,-40.,40.,10*n2,-10.,10.); |
| 36 |
|
|
h_par_fine[1] = new TH2F("houghparfine_1","par1 vs par2",10*n1,-40.,40.,10*n2,-10.,10.); |
| 37 |
|
|
|
| 38 |
|
|
h_par1[0] = new TH1F("houghpar1_0","par1",n1,-10.,10.); |
| 39 |
|
|
h_par1[1] = new TH1F("houghpar1_1","par1",n1,-10.,10.); |
| 40 |
|
|
|
| 41 |
|
|
h_par2[0] = new TH1F("houghpar2_0","par2",n2,-40.,40.); |
| 42 |
|
|
h_par2[1] = new TH1F("houghpar2_1","par2",n2,-40.,40.); |
| 43 |
|
|
}; |
| 44 |
|
|
//////////////////////////////////////////////////////////// |
| 45 |
|
|
/// CLASS IMPLEMENTATION |
| 46 |
|
|
//////////////////////////////////////////////////////////// |
| 47 |
|
|
/** |
| 48 |
|
|
* Reset event |
| 49 |
|
|
* it has to be called for each event |
| 50 |
|
|
*/ |
| 51 |
|
|
void CaloHough::Reset(){ |
| 52 |
|
|
return; |
| 53 |
|
|
} |
| 54 |
|
|
|
| 55 |
|
|
/** |
| 56 |
|
|
* \brief Set event |
| 57 |
|
|
*/ |
| 58 |
|
|
bool CaloHough::Set(CaloLevel1 *cl1){ |
| 59 |
|
|
|
| 60 |
|
|
// bool debug = CaloHough_parameters::Get()->debug; |
| 61 |
|
|
float Q_UP = CaloHough_parameters::Get()->Q_UP; |
| 62 |
|
|
int M_UP = CaloHough_parameters::Get()->M_UP; |
| 63 |
|
|
float Q_DN = CaloHough_parameters::Get()->Q_DN; |
| 64 |
|
|
int M_DN = CaloHough_parameters::Get()->M_DN; |
| 65 |
|
|
|
| 66 |
|
|
cout << "bool CaloHough::Set(CaloLevel1*)"<<endl; |
| 67 |
|
|
|
| 68 |
|
|
if(!cl1)return false; |
| 69 |
|
|
|
| 70 |
|
|
TH2F* h_qtot[2]; |
| 71 |
|
|
TH2F* h_par[2]; |
| 72 |
|
|
TH2F* h_par_fine[2]; |
| 73 |
|
|
TH1F* h_par1[2]; |
| 74 |
|
|
TH1F* h_par2[2]; |
| 75 |
|
|
for(int iv=0; iv<2; iv++){ |
| 76 |
|
|
h_qtot[iv] = CaloHough_parameters::Get()->h_qtot[iv]; |
| 77 |
|
|
h_par[iv] = CaloHough_parameters::Get()->h_par[iv]; |
| 78 |
|
|
h_par_fine[iv] = CaloHough_parameters::Get()->h_par_fine[iv]; |
| 79 |
|
|
h_par1[iv] = CaloHough_parameters::Get()->h_par1[iv]; |
| 80 |
|
|
h_par2[iv] = CaloHough_parameters::Get()->h_par2[iv]; |
| 81 |
|
|
if(!h_qtot[iv] || !h_par[iv] || !h_par1[iv] || !h_par2[iv] || !h_par_fine[iv])return false; |
| 82 |
|
|
h_qtot[iv]->Reset(); |
| 83 |
|
|
h_par[iv]->Reset(); |
| 84 |
|
|
h_par_fine[iv]->Reset(); |
| 85 |
|
|
h_par1[iv]->Reset(); |
| 86 |
|
|
h_par2[iv]->Reset(); |
| 87 |
|
|
} |
| 88 |
|
|
|
| 89 |
|
|
// ----------------------------------- |
| 90 |
|
|
// set charge collected in each plane |
| 91 |
|
|
// ----------------------------------- |
| 92 |
|
|
// |
| 93 |
|
|
CaloStrip cstrip; |
| 94 |
|
|
int isimu = CaloHough_parameters::Get()->isimu; |
| 95 |
|
|
if(isimu==0)cstrip = CaloStrip(cl1,false); |
| 96 |
|
|
if(isimu==1)cstrip = CaloStrip(cl1,true); |
| 97 |
|
|
// |
| 98 |
|
|
vector<int> ids[22][2]; //id strip 0-95 |
| 99 |
|
|
vector<float> xy[22][2]; //cluster coordinate |
| 100 |
|
|
int islast = -1; |
| 101 |
|
|
int iplast = -1; |
| 102 |
|
|
int mult = 0; |
| 103 |
|
|
float qsum = 0.; |
| 104 |
|
|
float xyqsum = 0.; |
| 105 |
|
|
float issum = 0.; |
| 106 |
|
|
int ncluster[] = {0,0}; |
| 107 |
|
|
for(int ih=0;ih<cl1->istrip;ih++){ |
| 108 |
|
|
// |
| 109 |
|
|
int iv=-1; |
| 110 |
|
|
int ip=-1; |
| 111 |
|
|
int is=-1; |
| 112 |
|
|
cl1->DecodeEstrip(ih,iv,ip,is); |
| 113 |
|
|
float xyz[3]; |
| 114 |
|
|
cstrip.Set(iv,ip,is); |
| 115 |
|
|
// cout << " * "<<setw(3)<<iv<<setw(3)<<ip<<setw(3)<<is<<endl; |
| 116 |
|
|
xyz[0] = cstrip.GetX(); |
| 117 |
|
|
xyz[1] = cstrip.GetY(); |
| 118 |
|
|
xyz[2] = cstrip.GetZ(); |
| 119 |
|
|
// |
| 120 |
|
|
// ...................................... metodo 2 |
| 121 |
|
|
if( is != islast+1 || ip != iplast ){ // new cluster |
| 122 |
|
|
if( |
| 123 |
|
|
qsum < Q_UP && |
| 124 |
|
|
qsum > Q_DN && |
| 125 |
|
|
mult < M_UP && |
| 126 |
|
|
mult > M_DN && |
| 127 |
|
|
true){ // store the previous cluster |
| 128 |
|
|
xy[ip][iv].push_back( xyqsum/qsum ); |
| 129 |
|
|
ids[ip][iv].push_back( (int)(0.5+issum/qsum) ); |
| 130 |
|
|
h_qtot[iv]->Fill( issum/qsum, (float)(21-ip), qsum); |
| 131 |
|
|
ncluster[iv]++; |
| 132 |
|
|
cout << " >> "<<qsum<<endl; |
| 133 |
|
|
}// ... and reset counters |
| 134 |
|
|
mult = 0; |
| 135 |
|
|
qsum = 0.; |
| 136 |
|
|
xyqsum = 0.; |
| 137 |
|
|
issum = 0.; |
| 138 |
|
|
} |
| 139 |
|
|
|
| 140 |
|
|
mult++; |
| 141 |
|
|
qsum += cstrip.GetE(); |
| 142 |
|
|
xyqsum += xyz[iv] * cstrip.GetE(); |
| 143 |
|
|
issum += (float)is * cstrip.GetE(); |
| 144 |
|
|
islast = is; |
| 145 |
|
|
iplast = ip; |
| 146 |
|
|
|
| 147 |
|
|
// ...................................... metodo 1 |
| 148 |
|
|
// if( cstrip.GetE() < Q_UP ){ |
| 149 |
|
|
// h_qtot[iv]->Fill( (float)is, (float)(21-ip), cstrip.GetE()); |
| 150 |
|
|
// ids[ip][iv].push_back(is); |
| 151 |
|
|
// } |
| 152 |
|
|
|
| 153 |
|
|
} |
| 154 |
|
|
cout << "n.clusters X "<<ncluster[0]<<endl; |
| 155 |
|
|
cout << "n.clusters Y "<<ncluster[0]<<endl; |
| 156 |
|
|
// |
| 157 |
|
|
// cstrip.Set(0,0,0); |
| 158 |
|
|
// cstrip.Set(0,21,0); |
| 159 |
|
|
cstrip.Set(0,11,0); |
| 160 |
|
|
float zref = cstrip.GetZ();//reference plane |
| 161 |
|
|
// |
| 162 |
|
|
for(int iv=0; iv<2; iv++){ |
| 163 |
|
|
for(int ip1=0; ip1<22-1; ip1++){ |
| 164 |
|
|
// for(int ip2=ip1+1; ip2<22; ip2++){ |
| 165 |
|
|
for(int ip2=ip1+1; ip2<min(22,ip1+3); ip2++){ |
| 166 |
|
|
// cout << "ip1 "<<ip1<<" "<<ids[ip1][iv].size()<<endl; |
| 167 |
|
|
// cout << "ip2 "<<ip2<<" "<<ids[ip2][iv].size()<<endl; |
| 168 |
|
|
for(int is1=0; is1<ids[ip1][iv].size(); is1++){ |
| 169 |
|
|
for(int is2=0; is2<ids[ip2][iv].size(); is2++){ |
| 170 |
|
|
|
| 171 |
|
|
// ...................................... metodo 1 |
| 172 |
|
|
// int ids1 = ids[ip1][iv][is1]; |
| 173 |
|
|
// int ids2 = ids[ip2][iv][is2]; |
| 174 |
|
|
// float xyz1[3]; |
| 175 |
|
|
// float xyz2[3]; |
| 176 |
|
|
// cstrip.Set(iv,ip1,ids1); |
| 177 |
|
|
// xyz1[0] = cstrip.GetX(); |
| 178 |
|
|
// xyz1[1] = cstrip.GetY(); |
| 179 |
|
|
// xyz1[2] = cstrip.GetZ(); |
| 180 |
|
|
// cstrip.Set(iv,ip2,ids2); |
| 181 |
|
|
// xyz2[0] = cstrip.GetX(); |
| 182 |
|
|
// xyz2[1] = cstrip.GetY(); |
| 183 |
|
|
// xyz2[2] = cstrip.GetZ(); |
| 184 |
|
|
|
| 185 |
|
|
// ...................................... metodo 2 |
| 186 |
|
|
int ids1 = ids[ip1][iv][is1]; |
| 187 |
|
|
int ids2 = ids[ip2][iv][is2]; |
| 188 |
|
|
float xyz1[3]; |
| 189 |
|
|
float xyz2[3]; |
| 190 |
|
|
cstrip.Set(iv,ip1,is1); |
| 191 |
|
|
xyz1[iv] = xy[ip1][iv][is1]; |
| 192 |
|
|
xyz1[2] = cstrip.GetZ(); |
| 193 |
|
|
cstrip.Set(iv,ip2,ids2); |
| 194 |
|
|
xyz2[iv] = xy[ip2][iv][is2]; |
| 195 |
|
|
xyz2[2] = cstrip.GetZ(); |
| 196 |
|
|
// |
| 197 |
|
|
float par1 = 0.; |
| 198 |
|
|
if(fabs(xyz1[2]-xyz2[2])>0.)par1 = (xyz1[iv]-xyz2[iv])/(xyz1[2]-xyz2[2]); |
| 199 |
|
|
float par2 = xyz1[iv] - par1*xyz1[2]; |
| 200 |
|
|
par2 = par1*zref + par2; |
| 201 |
|
|
h_par[iv]->Fill(par2,par1); // parameter space |
| 202 |
|
|
h_par_fine[iv]->Fill(par2,par1); // parameter space |
| 203 |
|
|
h_par1[iv]->Fill(par1); |
| 204 |
|
|
h_par2[iv]->Fill(par2); |
| 205 |
|
|
// cout << setw(15)<<par1<< setw(15)<<par2<<endl; |
| 206 |
|
|
} |
| 207 |
|
|
} |
| 208 |
|
|
} |
| 209 |
|
|
} |
| 210 |
|
|
// |
| 211 |
|
|
for( int ibx=1; ibx <= h_par[iv]->GetXaxis()->GetNbins(); ibx++){ |
| 212 |
|
|
for( int iby=1; iby <= h_par[iv]->GetYaxis()->GetNbins(); iby++){ |
| 213 |
|
|
// get bin content |
| 214 |
|
|
float content = h_par[iv]->GetBinContent(ibx,iby); |
| 215 |
|
|
if( content > 6){ |
| 216 |
|
|
double xmi = h_par[iv]->GetXaxis()->GetBinLowEdge(ibx); |
| 217 |
|
|
double xma = h_par[iv]->GetXaxis()->GetBinUpEdge(ibx); |
| 218 |
|
|
double ymi = h_par[iv]->GetYaxis()->GetBinLowEdge(iby); |
| 219 |
|
|
double yma = h_par[iv]->GetYaxis()->GetBinUpEdge(iby); |
| 220 |
|
|
h_par_fine[iv]->GetXaxis()->SetRangeUser(xmi,xma); |
| 221 |
|
|
h_par_fine[iv]->GetYaxis()->SetRangeUser(ymi,yma); |
| 222 |
|
|
double par1 = h_par_fine[iv]->GetMean(1); |
| 223 |
|
|
double par2 = h_par_fine[iv]->GetMean(2); |
| 224 |
|
|
cout << iv << " <><><><><><><>< "<<setw(15)<<par1<<setw(15)<<par2<<endl; |
| 225 |
|
|
} |
| 226 |
|
|
} |
| 227 |
|
|
} |
| 228 |
|
|
|
| 229 |
|
|
// |
| 230 |
|
|
for(int ip=0; ip<22; ip++){ |
| 231 |
|
|
ids[ip][iv].clear(); |
| 232 |
|
|
xy[ip][iv].clear(); |
| 233 |
|
|
} |
| 234 |
|
|
}//end loop on views |
| 235 |
|
|
|
| 236 |
|
|
return true; |
| 237 |
|
|
} |
| 238 |
|
|
ClassImp(CaloHough); |
| 239 |
|
|
ClassImp(CaloHough_parameters); |
| 240 |
|
|
|
Parent Directory
| 
Revision Log