TrackerMap/bin/tkScanBuild.C

0001 // Root macro to create the 41 images of a Tracker Scan

0002 // Code developed with the help of Raffaella Radogna physics student of

0003 //  Bari University

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){ // separated modules

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.;}//inner

0051     }
0052   }else{ //overlayed modules

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){//endcap

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){//pix endcap

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 { //barrel

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   //for(int j=0; j<5; j++){std::cout<<xpol[j]<<" "<<ypol[j]<<std::endl;}

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   //Store data about tracker modules

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)))){//use this to select stereo modules

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