File indexing completed on 2024-04-06 12:01:12
0001
0002
0003
0004 #include <map>
0005 #include <math.h>
0006 #include <iostream>
0007 #include <fstream>
0008 #include <sstream>
0009 #include <string>
0010 #include <vector>
0011 #include "TCanvas.h"
0012 #include "TPad.h"
0013 #include "TPolyLine.h"
0014 #include "TPolyMarker.h"
0015 #include <map>
0016 #include "TROOT.h"
0017
0018
0019 bool posrel;
0020 bool saveAsSingleLayer;
0021 int nlay,key;
0022 double xmin,xmax,ymin,ymax;
0023 int xsize,ysize,ix,iy;
0024 int cwidth,cheight;
0025 TCanvas * MyGC;
0026
0027 struct coordinateModulo{
0028 int part;
0029 int subpart;
0030 int layer;
0031 int ring;
0032 int nmod;
0033 double posx, posy, posz;
0034 double length, width, thickness, widthAtHalfLength;
0035 int detid;
0036 std::string nome;
0037 } mod[16588];
0038
0039 void defwindow(int num_lay){
0040 nlay = num_lay;
0041 if(posrel){
0042 xmin=-2.;ymin=-2.;xmax=2.;ymax=2.;
0043 if(nlay >12 && nlay < 19){
0044 xmin=-.40;xmax=.40;ymin=-.40;ymax=.40;
0045 }
0046 if(nlay>30){
0047 xmin=-0.1;xmax=3.;ymin=-0.1;ymax=8.5;
0048 if(nlay<34){xmin=-0.3;xmax=1.0;}
0049 if(nlay>33&&nlay<38){xmax=2.0;}
0050 if(nlay>37){ymax=8.;}
0051 }
0052 }else{
0053 xmin=-1.3;ymin=-1.3;xmax=1.3;ymax=1.3;
0054 if(nlay >12 && nlay < 19){
0055 xmin=-.20;xmax=.20;ymin=-.20;ymax=.20;
0056 }
0057 if(nlay>30){
0058 xmin=-1.5;xmax=1.5;ymin=-1.;ymax=28.;
0059 if(nlay<34){xmin=-0.5;xmax=0.5;}
0060 if(nlay>33&&nlay<38){xmin=-1.;xmax=1.;}
0061 }
0062
0063 }
0064 if(nlay<16){
0065 ix=0;iy=(15-nlay)*ysize;}
0066 if(nlay>15&&nlay<31){
0067 ix=3*xsize;iy=(nlay-16)*ysize;}
0068 if(nlay>30){
0069 if(nlay==31){ix=(int)(1.5*xsize);iy=0;}
0070 if(nlay>31 && nlay%2==0){int il=(nlay-30)/2;ix=xsize;iy=il*2*ysize;}
0071 if(nlay>31 && nlay%2!=0){int il=(nlay-30)/2;ix=2*xsize;iy=il*2*ysize;}
0072 }
0073 }
0074
0075 double xdpixel(double x){
0076 double res;
0077 if(saveAsSingleLayer)res= ((x-xmin)/(xmax-xmin)*xsize);
0078 else res= ((x-xmin)/(xmax-xmin)*xsize)+ix;
0079 return res;
0080 }
0081
0082 double ydpixel(double y){
0083 double res;
0084 double y1;
0085 y1 = (y-ymin)/(ymax-ymin);
0086 if(nlay>30) res= 2*ysize - (y1*2*ysize);
0087 else res= ysize - (y1*ysize);
0088 if(!saveAsSingleLayer) res=res+iy;
0089 return res;
0090 }
0091
0092 double phival(double x, double y){
0093 double phi;
0094 double phi1=atan(y/x);
0095 phi = phi1;
0096 if(y<0. && x>0) phi = phi1+2.*M_PI;
0097 if(x<0.)phi=phi1+M_PI;
0098 if(fabs(y)<0.000001 && x>0)phi=0;
0099 if(fabs(y)<0.000001&&x<0)phi=M_PI;
0100 if(fabs(x)<0.000001&&y>0)phi=M_PI/2.;
0101 if(fabs(x)<0.000001&&y<0)phi=3.*M_PI/2.;
0102 return phi;
0103 }
0104 bool isRingStereo(int key){
0105 int layer=key/100000;
0106 int ring = key - layer*100000;
0107 ring = ring/1000;
0108 if(layer==34 || layer==35 || layer==38 || layer==39) return true;
0109 if(layer<13 || (layer>18&&layer<31))
0110 if(ring==1 || ring==2 || ring==5)return true;
0111 return false;
0112 }
0113
0114
0115 void computemodule(int modulo,int nlay,int *npoints,double *xpol, double *ypol ){
0116 double phi,r;
0117 double xp[4],yp[4],xp1,yp1;
0118 double vhbot,vhtop,vhapo;
0119 double xt1,yt1,xs1=0.,ys1=0.,xt2,yt2,xs2,ys2,pv1,pv2;
0120 *npoints=5;
0121 phi = phival(mod[modulo].posx,mod[modulo].posy);
0122 r = sqrt(mod[modulo].posx*mod[modulo].posx+mod[modulo].posy*mod[modulo].posy);
0123 vhbot = mod[modulo].width;
0124 vhtop=mod[modulo].width;
0125 vhapo=mod[modulo].length;
0126 if(nlay < 31){
0127 vhbot = mod[modulo].widthAtHalfLength/2.-(mod[modulo].width/2.-mod[modulo].widthAtHalfLength/2.);
0128 vhtop=mod[modulo].width/2.;
0129 vhapo=mod[modulo].length/2.;
0130 if(nlay >12 && nlay <19){
0131 xp[0]=r-vhtop;yp[0]=-vhapo;
0132 xp[1]=r+vhtop;yp[1]=-vhapo;
0133 xp[2]=r+vhtop;yp[2]=vhapo;
0134 xp[3]=r-vhtop;yp[3]=vhapo;
0135 }else{
0136 xp[0]=r-vhapo;yp[0]=-vhbot;
0137 xp[1]=r+vhapo;yp[1]=-vhtop;
0138 xp[2]=r+vhapo;yp[2]=vhtop;
0139 xp[3]=r-vhapo;yp[3]=vhbot;
0140 }
0141 for(int j=0;j<4;j++){
0142 xp1 = xp[j]*cos(phi)-yp[j]*sin(phi);
0143 yp1 = xp[j]*sin(phi)+yp[j]*cos(phi);
0144 xp[j] = xp1;yp[j]=yp1;
0145 }
0146 } else {
0147 int numod;
0148 numod=mod[modulo].nmod;if(numod>100)numod=numod-100;
0149 xt1=r; yt1=-vhtop/2.;
0150 xs1 = xt1*cos(phi)-yt1*sin(phi);
0151 ys1 = xt1*sin(phi)+yt1*cos(phi);
0152 xt2=r; yt2=vhtop/2.;
0153 xs2 = xt2*cos(phi)-yt2*sin(phi);
0154 ys2 = xt2*sin(phi)+yt2*cos(phi);
0155 pv1=phival(xs1,ys1);
0156 pv2=phival(xs2,ys2);
0157 if(fabs(pv1-pv2)>M_PI && numod==1)pv1=pv1-2.*M_PI;
0158 if(fabs(pv1-pv2)>M_PI && numod!=1)pv2=pv2+2.*M_PI;
0159 xp[0]=mod[modulo].posz-vhapo/2.;yp[0]=4.2*pv1;
0160 xp[1]=mod[modulo].posz+vhapo/2.;yp[1]=4.2*pv1;
0161 xp[2]=mod[modulo].posz+vhapo/2. ;yp[2]=4.2*pv2;
0162 xp[3]=mod[modulo].posz-vhapo/2.;yp[3]=4.2*pv2;
0163 }
0164 for(int j=0;j<4;j++){
0165 ypol[j]=xdpixel(xp[j]);xpol[j]=ydpixel(yp[j]);
0166 }
0167 ypol[4]=ypol[0];xpol[4]=xpol[0];
0168
0169 }
0170
0171
0172 int main(int argc, char *argv[]){
0173
0174 char *inputFile, *outputFile;
0175 if(argc>1)
0176 inputFile=argv[1];
0177 if(argc>2)
0178 outputFile=argv[2];
0179 xsize=340; ysize=200;
0180 cwidth=2000; cheight=1000;
0181 TPolyMarker *PM[43];
0182 std::string layerName[43]={"TECM9.png","TECM8.png","TECM7.png","TECM6.png","TECM5.png","TECM4.png","TECM3.png","TECM2.png","TECM1.png","TIDM3.png","TIDM2.png","TIDM1.png"," ","PIXEM2.png","PIXEM1.png","PIXEB1.png","PIXEB2.png"," ","TIDP1.png","TIDP2.png","TIDP3.png","TECP1.png","TECP2.png","TECP3.png","TECP4.png","TECP5.png","TECP6.png","TECP7.png","TECP8.png","TECP9.png","PIXB1.png","PIXB2.png","PIXB3.png","TIB1.png","TIB2.png","TIB3.png","TIB4.png","TOB1.png","TOB2.png","TOB3.png","TOB4.png","TOB5.png","TOB6.png"};
0183
0184
0185 std::ifstream *cfile;
0186 int cont;
0187
0188 std::string ciccio;
0189
0190 std::map<int,int>indice_moduli;
0191
0192 cfile = new std::ifstream("tracker.dat",std::ios::in);
0193
0194 Int_t np = 0;
0195
0196 while(!cfile->eof()) {
0197 *cfile >> cont >> mod[np].part >> mod[np].subpart
0198 >> mod[np].layer >> mod[np].ring >> mod[np].nmod
0199 >> mod[np].posx >> mod[np].posy >> mod[np].posz
0200 >> mod[np].length >> mod[np].width >> mod[np].thickness
0201 >> mod[np].widthAtHalfLength >> mod[np].detid;
0202 getline(*cfile,ciccio);
0203 getline(*cfile,mod[np].nome);
0204 indice_moduli[mod[np].detid]=np;;
0205
0206 np++;
0207
0208 }
0209 std::cout << "tracker.dat processed " << std::endl;
0210 MyGC = new TCanvas("MyGC","Tracker Layer",cwidth,cheight);
0211 for(int i=0;i<43;i++){
0212 PM[i] = new TPolyMarker; PM[i]->SetMarkerStyle(1); PM[i]->SetMarkerColor(1);
0213 }
0214 int modulo;
0215 double xp,yp,zp;
0216 int id;
0217 double * x = new double[5];
0218 double * y = new double[5];
0219 int npoints;
0220 std::ifstream *bfile;
0221 gPad->SetFillColor(10);
0222 gPad->Range(-10,0,400,340);
0223
0224
0225 bfile = new std::ifstream(inputFile,std::ios::in);
0226
0227 np=0;
0228 saveAsSingleLayer=true;
0229 posrel=false;
0230
0231 while(!bfile->eof()) {
0232 *bfile >> id >> xp >> yp >> zp ;
0233 np++;
0234
0235 modulo= indice_moduli[id];
0236 if(indice_moduli.find(id)!=indice_moduli.end() && mod[modulo].layer>0 && mod[modulo].layer<44){
0237 key=mod[modulo].layer*100000+mod[modulo].ring*1000+mod[modulo].nmod;
0238 if( (mod[modulo].nmod>100||
0239 (mod[modulo].nmod<100&&!isRingStereo(key))||
0240 (mod[modulo].nmod<100&&
0241 (mod[modulo].layer==36||mod[modulo].layer==37||mod[modulo].layer==40||mod[modulo].layer==41||mod[modulo].layer==42||mod[modulo].layer==43)))){
0242
0243 defwindow(mod[modulo].layer);
0244 if(mod[modulo].layer<31)PM[mod[modulo].layer -1]->SetNextPoint(ydpixel(yp/100.),xdpixel(xp/100.));
0245 else {PM[mod[modulo].layer -1]->SetNextPoint(ydpixel(4.2*phival(xp/100.,yp/100.)),xdpixel(zp/100.));
0246
0247 }
0248 }
0249 if (np% 1000000==0)
0250 std::cout<<np<<" points processed" << std::endl;
0251 }
0252 }
0253 bfile->close();
0254
0255 int sel_layer;
0256
0257 for(int i=0;i<43;i++){
0258 sel_layer=i+1;
0259 defwindow(sel_layer);
0260 for (int k=0; k<16588; k++){
0261
0262 if(mod[k].layer==sel_layer) {
0263 computemodule(k,sel_layer,&npoints,x,y);
0264 TPolyLine* pline = new TPolyLine(npoints,x,y);
0265 pline->Draw();
0266 }
0267
0268 }
0269 PM[i]->Draw();
0270 MyGC->Update();
0271 if(layerName[i]!=" ")MyGC->Print(layerName[i].c_str());
0272 MyGC->Clear();
0273 delete PM[i];
0274 }
0275 return 0;
0276
0277 }
0278
0279