CMSSW Reference Manual

00001 #include "FWCore/Framework/interface/EventSetup.h"
00002 #include "FWCore/Framework/interface/ESHandle.h"
00003 #include "VisReco/VisCustomTracker/interface/VisCuTkBuilder.h"
00004 #include "VisReco/VisCustomTracker/interface/VisCuTracker.h"
00005 #include "VisReco/VisCustomTracker/interface/VisCuTkModuleMap.h"
00006 #include "VisReco/VisCustomTracker/interface/VisCuCmsTracker.h"
00007 #include "VisReco/VisCustomTracker/interface/VisCuTkModule.h"
00008 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
00009 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00010 #include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
00011 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00012 
00013 #include <fstream>
00014 #include <iostream>
00015 
00016 VisCuCmsTracker* VisCuTkBuilder::tracker()
00017 {
00018         return theTracker;
00019 }
00020 
00021 void  VisCuTkBuilder::create(){
00022     ntotmod = 0;
00023     theTracker = new VisCuCmsTracker();
00024 }
00025 
00026 void VisCuTkBuilder::getPos(int pix_sil,int fow_bar,int layer, int ring,int nmod,int &idetector,int &ipart,int &ilayer, int &iring,int &module)
00027 {
00028   int spicchif[] ={24,24,40,56,40,56,80};
00029   int spicchib[] ={20,32,44,30,38,46,56,42,48,54,60,66,74};
00030   
00031   idetector = pix_sil;  
00032   if(idetector==2){
00033     if (((layer<=12 && layer>=10) || (layer<=21 && layer>=19) || (layer<=37 && layer>=34))){
00034       idetector=2;
00035     }
00036     else {idetector=3;}
00037   }
00038   
00039   ipart = fow_bar;
00040   ilayer = layer;
00041   //inner silicon encap
00042   if(idetector==2 && ipart==1){
00043     ilayer=ilayer-9;
00044   } 
00045   if(idetector==2 && ipart==3){
00046     ilayer=ilayer-18;
00047   }      
00048   //pixel encap
00049   if(idetector==1 && ipart==1){
00050     ilayer=ilayer-12;
00051   }
00052   if(idetector==1 && ipart==3){
00053     ilayer=ilayer-15; 
00054   }
00055   //outer silicon encap
00056   if(idetector==3 && ipart==3){
00057     ilayer=ilayer-21;
00058   }
00059   //pixel barrel
00060   if(idetector==1 && ipart==2){
00061     ilayer=ilayer-30;
00062   }
00063   //inner silicon barrel
00064   if(idetector==2 && ipart==2){
00065     ilayer=ilayer-33;
00066   }
00067   //outer silicon barrel
00068   if(idetector==3 && ipart==2){
00069     ilayer=ilayer-37;
00070   }
00071   iring = ring;
00072   //outer silicon endcap
00073   if(idetector==3 && (ipart==3 || ipart==1)){
00074     if(layer==1 || layer==30) iring =ring-3;
00075     if(layer==2 || layer==29) iring =ring-2;
00076     if(layer==3 || layer==28) iring =ring-2;
00077     if(layer==4 || layer==27) iring =ring-1;
00078     if(layer==5 || layer==26) iring =ring-1;
00079     if(layer==6 || layer==25) iring =ring-1;
00080     if(layer==7 || layer==24) iring =ring;
00081     if(layer==8 || layer==23) iring =ring;
00082     if(layer==9 || layer==22) iring =ring;
00083   }
00084   if(idetector>1 && ipart==2 && ring < 7)iring=7-ring; //TIB,TOB ring numbering in barrel -z should be reversed
00085 if(idetector==2 && ipart==2 && ring < 7){ if((iring%2)==0)iring=iring-1; else iring=iring+1;}//more corrections for the TIB numbering
00086  if(theTracker->type==2)//correction only for MTCC  
00087   {
00088 if(idetector==2 && ipart==2 && ring > 6){ if((iring%2)==0)iring=iring-1; else iring=iring+1;}
00089   }
00090   module = nmod;
00091   if(module >100){
00092     if(ipart==2){ //barrel
00093       module=(module-100)+spicchib[layer-31];
00094       if(theTracker->type==2&&layer==35)module=(module-100)+spicchib[layer-30];
00095     }
00096     if(ipart!=2){//endcap
00097       module=(module-100)+spicchif[ring-1];
00098     }
00099   }
00100 }
00101 
00102 
00103 void VisCuTkBuilder::fill(const edm::ESHandle<TrackerGeometry> pDD){
00104   std::cout << "domenico ESHANDle " << pDD.product() << std::endl;
00105 
00106   fill( pDD.product() );
00107 }
00108 
00109 void VisCuTkBuilder::fill(const TrackerGeometry* pDD){
00110   int spicchif[] ={24,24,40,56,40,56,80};
00111   int spicchib[] ={20,32,44,30,38,46,56,42,48,54,60,66,74};
00112   pset=VisConfigurationService::pSet();
00113   if(!pset)
00114     {
00115       edm::ParameterSet p;
00116       VisConfigurationService *visService = new VisConfigurationService(p);
00117       pset = visService->pSet();
00118       delete visService;
00119     } 
00120   std::string saveTrackerdat = pset->getUntrackedParameter<std::string>("saveTrackerdat", "false");
00121    std::ofstream *  output = 0;
00122    if(saveTrackerdat=="true")output = new std::ofstream("tracker.dat",std::ios::out);
00123   // edm::ESHandle<TrackingGeometry> pDD;
00124   // iSetup.get<TrackerDigiGeometryRecord>().get( pDD );
00125   std::cout << " Geometry node for CmsDigiTracker is  "<<pDD<<std::endl;
00126   std::cout << " theTracker "<<theTracker->components()<<std::endl;
00127   
00128   std::vector<GeomDetUnit*>::const_iterator begin = pDD->detUnits().begin();
00129   std::vector<GeomDetUnit*>::const_iterator end = pDD->detUnits().end();
00130   cout <<"ntotmod="<<  pDD->detUnits().size() << endl;
00131 
00132   // Find where each module should be loaded
00133   buildMap(pDD);
00134 
00135   int iLayer=0;
00136   int iModule;
00137   int old_layer=0;
00138   int ring = 0;
00139   int nmod = 0;
00140   int ntotmod = 0;
00141   float r;
00142   int bar_fow = 1;
00143   float phi,phi1;
00144   float rmedioS[]={0.27665, 0.3671, 0.4474, 0.5617, 0.6768, 0.8189, 0.9907};
00145   float rmedioP[]={0.0623081, 0.074111,  0.0870344, 0.103416, 0.115766, 0.132728, 0.140506}; 
00146  
00147   int det_type;
00148   int nlay=0;
00149   int layer,subdet,leftright=0,ringno,petalno,moduleno,isStereo;
00150   subdet = (*begin)->geographicalId().subdetId();
00151   for ( ; begin != end; ++begin) {
00152     ntotmod++;
00153     subdet = (*begin)->geographicalId().subdetId();
00154     if(subdet==1){//pixel barrel
00155       layer = ((*begin)->geographicalId().rawId()>>16)&0xF;
00156       leftright=0;
00157       }
00158     if(subdet==3){//TIB
00159       layer = ((*begin)->geographicalId().rawId()>>14)&0x7;
00160       leftright=0;
00161       }
00162     if(subdet==5){//TOB
00163       layer = ((*begin)->geographicalId().rawId()>>14)&0x7;
00164       leftright=0;
00165       }
00166     if(subdet==2){//Pixel endcap
00167       leftright = ((*begin)->geographicalId().rawId()>>23)&0x3;
00168       layer = ((*begin)->geographicalId().rawId()>>16)&0xF;
00169     }
00170     if(subdet==4){//TID
00171       leftright = ((*begin)->geographicalId().rawId()>>13)&0x3;
00172       layer = ((*begin)->geographicalId().rawId()>>11)&0x3;
00173     }
00174     if(subdet==6){//TEC
00175       leftright = ((*begin)->geographicalId().rawId()>>18)&0x3;
00176       layer = ((*begin)->geographicalId().rawId()>>14)&0xF;
00177     }
00178     isStereo=(*begin)->geographicalId().rawId()&0x3;
00179     nlay = layerno(subdet,leftright,layer);
00180     ringno=0;
00181     petalno=0;
00182     moduleno=0;
00183     if(subdet==1){ringno=((*begin)->geographicalId().rawId()>>8)&0xFF;}
00184     if(subdet==2){ringno=((*begin)->geographicalId().rawId()>>8)&0xFF;}
00185     if(subdet==3){ringno=((*begin)->geographicalId().rawId()>>4)&0x3F;}
00186     if(subdet==5){ringno=((*begin)->geographicalId().rawId()>>5)&0x7F;}
00187     if(subdet==4){ringno=((*begin)->geographicalId().rawId()>>9)&0x3;}
00188     if(subdet==6){ringno=((*begin)->geographicalId().rawId()>>5)&0x7;}
00189     if(subdet==6){ petalno=((*begin)->geographicalId().rawId()>>8)&0xF;}
00190     if(subdet==2){ petalno=((*begin)->geographicalId().rawId()>>10)&0x3F;}
00191     if(subdet==1){moduleno=((*begin)->geographicalId().rawId() >>2)&0x3F; }
00192     if(subdet==3){moduleno=((*begin)->geographicalId().rawId() >>2)&0x3; }
00193     if(subdet==5){moduleno=((*begin)->geographicalId().rawId() >>2)&0x7; }
00194     if(subdet==4){moduleno=((*begin)->geographicalId().rawId() >>2)&0x1F; }
00195     if(subdet==6){moduleno=((*begin)->geographicalId().rawId() >>2)&0x7; }
00196     if(subdet==2){moduleno=((*begin)->geographicalId().rawId() >>2)&0x3F; }
00197     iLayer = nlay;
00198     if(iLayer !=old_layer){old_layer=iLayer;iModule=0;}
00199     float posx =  (*begin)->surface().position().x()/ 100.0;  // cm -> m
00200     float posy =  (*begin)->surface().position().y()/ 100.0;  // cm -> m
00201     float posz =  (*begin)->surface().position().z()/ 100.0;  // cm -> m
00202     bar_fow=2;
00203     if(nlay<16)bar_fow=1;
00204     if(nlay>15&&nlay<31)bar_fow=3;
00205     r=pow(((double)(posx*posx) + posy*posy),0.5);
00206     phi1=atan(posy/posx);
00207     phi=phi1;
00208     if(posy < 0.&&posx>0)phi=phi1+2.*M_PI;
00209     if(posx < 0.)phi=phi1+M_PI;
00210     if(fabs(posy)<0.000001&&posx>0)phi=0;
00211     if(fabs(posy)<0.000001&&posx<0)phi=M_PI;
00212     if(fabs(posx)<0.000001&&posy>0)phi=M_PI/2.;
00213     if(fabs(posx)<0.000001&&posy<0)phi=M_PI + M_PI/2.;
00214    
00215     if(bar_fow==2){ //barrel
00216       if(subdet==1)ring=moduleno;
00217       if(subdet==5){
00218         ring=moduleno;
00219         if((((*begin)->geographicalId().rawId()>>12)&0x3)==2)ring=ring+6;
00220       }
00221       if(subdet==3){
00222         ring=moduleno;
00223         if(layer==2||layer==4){
00224           if((((*begin)->geographicalId().rawId()>>10)&0x3)==2){ring=ring*2;} else {ring=ring*2-1;}
00225         }else{
00226           if((((*begin)->geographicalId().rawId()>>10)&0x3)==2){ring=ring*2-1;} else {ring=ring*2;}
00227         }
00228         if((((*begin)->geographicalId().rawId()>>12)&0x3)==2)ring=ring+6;
00229       }
00230 
00231       nmod=(int)((phi/(2.*M_PI))*spicchib[nlay-31]+.1)+1;
00232       if(theTracker->type==2&&nlay==35)nmod=(int)((phi/(2.*M_PI))*spicchib[nlay-30]+.1)+1;
00233       if(nlay==40)nmod = nmod-1;
00234       if(subdet==1)nmod=ringno;
00235     } else { // endcap
00236       if (subdet==4||subdet==6) { //endcap silicon strip
00237         for (int i=0;i<7; i++){
00238           if (fabs(r-rmedioS[i])<0.015){
00239             ring =i+1;
00240             break;
00241           }
00242         }
00243         nmod=(int)((phi/(2.*M_PI))*spicchif[ring-1]+.1)+1;
00244       } else{ // endcap pixel
00245         for (int i=0; i<7; i++){
00246           if(fabs(r -rmedioP[i])<0.0001){
00247             ring=i+1;
00248             break;
00249           }
00250         }
00251         nmod=(int)((phi/(2.*M_PI))*24+.1)+1;
00252       }
00253     } //end of endcap part
00254     if (isStereo==1) nmod=nmod+100;
00255 
00256     float  length = (*begin)->surface ().bounds ().length () / 100.0; // cm -> m
00257     float  width = (*begin)->surface ().bounds ().width () / 100.0;   // cm -> mo
00258     float  thickness = (*begin)->surface ().bounds ().thickness () / 100.0;  // cm -> m
00259     float  widthAtHalfLength=(*begin)->surface ().bounds ().widthAtHalfLength () / 100.0;
00260     
00261     det_type=1;
00262     if (subdet> 2)  det_type=2;
00263     int id = (*begin)->geographicalId().rawId();
00264     getPos(det_type,bar_fow,nlay,ring,nmod,idetector,ipart,ilayer,iring,module);
00265    if(saveTrackerdat=="true")(*output) <<ntotmod<< " "<< det_type << " " <<bar_fow<<" "<<nlay << " " << ring << " " << nmod <<" " << posx << " " <<posy << " " <<posz<<" "<<length<<" " <<width<<" "<<thickness<<" "<<widthAtHalfLength<<" "<<" "<<id<< std::endl;
00266     //output <<ntotmod<<" "<< det_type<<" "<<bar_fow<<" "<<nlay<<" "<<ring<<" "<<nmod<<" "<<idetector<<" "<<ipart<<" "<<ilayer<<" "<<iring<<" "<<module<<endl;
00267     mod = mmap[(*begin)->geographicalId().rawId()];
00268     if(mod==0){
00269       cout << "error in filling module: idetector="<<idetector <<" ipart="<<ipart<<" ilayer="<<ilayer<<" iring="<<iring<<" module="<< module <<endl;
00270       cout << " ind=" <<theTracker->getComponent(idetector)->getComponent(ipart)->getComponent(ilayer)->getComponent(iring)->getComponent(module)<< endl;
00271     }
00272     if(!mod->notInUse()){
00273       cout << "error in VisCuTkBuilder;  trying to reload module: idetector="<<idetector <<" ipart="<<ipart<<" ilayer="<<ilayer<<" iring="<<iring<<" module="<< module <<endl;
00274       cout << " ind=" <<theTracker->getComponent(idetector)->getComponent(ipart)->getComponent(ilayer)->getComponent(iring)->getComponent(module)<< endl;
00275     }   
00276     iModule++;
00277     string mod_name1 = VisCuTracker::getModuleName((*begin)->geographicalId().rawId());
00278     string mod_name = VisCuTracker::getModuleName((*begin)->geographicalId());
00279     if(mod_name!=mod_name1)std::cout << mod_name << std::endl << mod_name1 << std::endl;
00280     if(saveTrackerdat=="true")(*output) << mod_name << endl;
00281     mod->setName(mod_name); 
00282     mod->posx = posx;
00283     mod->posy = posy;
00284     mod->posz = posz;
00285     mod->length = length;
00286     mod->width = width;
00287     mod->thickness = thickness;
00288     mod->widthAtHalfLength = widthAtHalfLength;
00289     mod->value = 0;      
00290     mod->bufvalue = 0;   
00291     mod->detUnit = (*begin);
00292     VisCuTkModuleMap::moduleMap[(*begin)]=mod;
00293     mod->setVisible(false);
00294     mod->setUsed();
00295   }             
00296   mmap.clear(); //no more needed 
00297   }             
00298 void VisCuTkBuilder::buildMap(const TrackerGeometry* geom){
00299   //  ifstream* infile;
00300   infilename="CommonTools/TrackerMap/data/tracker.dat";
00301   int det_type, bar_fow, nlay,  ring, nmod;
00302   int nmods, idetector, ipart, iring, module, ilayer;
00303   unsigned int idex;
00304   float posx, posy, posz, length, width, thickness, widthAtHalfLength;
00305   int iModule=0,old_layer=0, ntotMod =0, number_modules;
00306   string name,dummys;
00307   ifstream infile(edm::FileInPath(infilename).fullPath().c_str(),ios::in);
00308   while(!infile.eof()) {
00309     infile >> nmods >> det_type >> bar_fow >> nlay  >> ring >> nmod >> posx >> posy
00310            >> posz>> length >> width >> thickness
00311            >> widthAtHalfLength >> idex ;
00312     getline(infile,dummys); //necessary to reach end of record
00313     getline(infile,name); 
00314     if(old_layer!=ilayer){old_layer=ilayer;iModule=0;}
00315     iModule++;
00316     ntotMod++;
00317     getPos(det_type,bar_fow,nlay,ring,nmod,idetector,ipart,ilayer,iring,module);
00318     mod = theTracker->getComponent(idetector)->getComponent(ipart)->getComponent(ilayer)->getComponent(iring)->getComponent(module);
00319     if(mod==0){
00320       cout << "error in module: idetector="<<idetector <<" ipart="<<ipart<<" ilayer="<<ilayer<<" iring="<<iring<<" module="<< module <<endl;
00321       cout << " ind=" <<theTracker->getComponent(idetector)->getComponent(ipart)->getComponent(ilayer)->getComponent(iring)->getComponent(module)<< endl;
00322     }
00323     else
00324       {
00325         mmap[idex]=mod;
00326       }
00327   }
00328   infile.close();
00329   number_modules = ntotMod-1;
00330   cout <<"ntotmod1="<<  number_modules << endl;
00331 }
VisCuTkBuilder.cc
