Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-02-14 12:46:07 

0001 #include "TText.h"
0002 #include "TPaveText.h"
0003 #include "TGeoManager.h"
0004 #include "TGeoMaterial.h"
0005 #include "TH1F.h"
0006 #include "TNtuple.h"
0007 #include "TRotation.h"
0008 #include "Riostream.h"
0009 #include "TStyle.h"
0010 #include "TLine.h" 
0011 #include "TCanvas.h"
0012 #include "TFile.h"
0013 #include "TTree.h"
0014 #include "TMath.h"
0015 #include "TSystem.h"
0016 #include "TError.h"
0017 #include <fstream>
0018 #include <iostream>
0019 #include <vector>
0020 #include <string>
0021 #include <algorithm>
0022 #include <assert.h>
0023 using namespace std;
0024 //using namespace ROOT::Math;
0025 
0026 // tree values
0027 int _nEntries;
0028 int _id, _level, _mid, _mlevel, _sublevel, _useid;
0029 float _xVal, _yVal, _zVal, _rVal, _phiVal, _alphaVal, _betaVal, _gammaVal;
0030 float _dxVal, _dyVal, _dzVal, _drVal, _dphiVal, _dalphaVal, _dbetaVal, _dgammaVal;
0031 float _surWidth, _surLength;
0032 int _identifiers[6];
0033 double _surRot[9];
0034 int _i;
0035 
0036 // global variables
0037 TTree* _inTree;
0038 int _sclftr;
0039 int _sclfrt;
0040 float _sclfmodulesizex;
0041 float _sclfmodulesizey;
0042 float _sclfmodulesizez;
0043 int _subdetector1;
0044 int _subdetector2;
0045 string _outputFileName;
0046 string _line1, _line2, _line3;
0047 float _piperadius;
0048 float _pipexcoord, _pipeycoord;
0049 float _linexcoord, _lineycoord;
0050 
0051 TTree *t;
0052 float zval;
0053 vector<int> vz;
0054 
0055 bool iszilessthanzj(int i, int j)
0056 {
0057     return vz[i] < vz[j];
0058 }
0059 
0060 void sortbyz(TString infile)
0061 {
0062     TFile *f = TFile::Open(infile);
0063     if(f == 0) {
0064         cout << "***Exception Thrown: null input file***" << endl;
0065         assert(0);
0066     }
0067     t = (TTree*)f->Get("alignTree");
0068     if(t == 0){
0069         cout << "***Exception Thrown: tree is null***" << endl;
0070         assert(0);
0071     }
0072     if(t->GetEntries() == 0) {
0073         cout << "***Exception Thrown: tree has no entries***" << endl;
0074         assert(0);
0075     }
0076 
0077     t->SetBranchAddress("z", &zval);
0078 
0079     TString outfile = infile.ReplaceAll(".root","_sorted.root");
0080     TFile *newf = TFile::Open(outfile, "RECREATE");
0081     TTree *newt = t->CloneTree(0);
0082 
0083     vector<int> v;
0084     vz.clear();
0085     int length = t->GetEntries();
0086     for (int i = 0; i < length; i++)
0087     {
0088         v.push_back(i);
0089 
0090         t->GetEntry(i);
0091         vz.push_back(zval);
0092     }
0093 
0094     sort(v.begin(), v.end(), iszilessthanzj);
0095 
0096     for (int i = 0; i < length; i++)
0097     {
0098         t->GetEntry(v[i]);
0099         newt->Fill();
0100     }
0101     newt->Write();
0102     delete newf;
0103 }
0104 
0105 void getBeamVisuals(TGeoManager* geom, TGeoVolume* top, float minZ, float maxZ) {
0106     TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0,0,0);
0107     TGeoMedium *Vacuum = new TGeoMedium("Vacuum",1, matVacuum);
0108     TGeoVolume *xyaxis = geom->MakeBox( "xyaxis", Vacuum, 90., 90., 40. );
0109 
0110     TGeoMaterial *matAl = new TGeoMaterial("Al", 26.98,13,2.7);
0111     TGeoMedium *Al = new TGeoMedium("Root Material",1, matAl);
0112     //TGeoVolume *line = geom->MakeTube( "BeamLine", Al, 0, .3, (maxZ - minZ) / 2 + 5);
0113     TGeoVolume *xaxis = geom->MakeTube( "XAxis", Al, 0, .1, 30.);
0114     TGeoVolume *yaxis = geom->MakeTube( "YAxis", Al, 0, .1, 30.);
0115     //TGeoVolume *pipe = geom->MakeTube( "BeamPipe", Al, _piperadius-.05, _piperadius+.05, (maxZ - minZ) / 2 + 5);
0116     //line->SetLineColor(kRed);
0117     xaxis->SetLineColor(kBlue);
0118     yaxis->SetLineColor(kBlue);
0119     //pipe->SetLineColor(kBlack);
0120 
0121     xyaxis->AddNode(xaxis, 1, new TGeoRotation( "rtyz", 0, 90, 0));
0122     xyaxis->AddNode(yaxis, 1, new TGeoRotation( "rtxz", 90, 90, 0));
0123     
0124     TGeoCombiTrans * pipecenter = new TGeoCombiTrans( *new TGeoTranslation(_pipexcoord, _pipeycoord, 0), *new TGeoRotation());
0125     //TGeoCombiTrans * linecenter = new TGeoCombiTrans( *new TGeoTranslation(_linexcoord, _lineycoord, 0), *new TGeoRotation());
0126     //top->AddNode( pipe, 1, pipecenter);
0127     //top->AddNode( line, 1, linecenter);
0128     top->AddNode( xyaxis, 1, pipecenter);
0129 }
0130 
0131 void getModule(TGeoManager* geom, TGeoVolume* top, TGeoVolume* mod){
0132 //--- define some materials
0133     TGeoMaterial *matAl = new TGeoMaterial("Al", 26.98,13,2.7);
0134 //--- define some media
0135     TGeoMedium *Al = new TGeoMedium("Root Material",1, matAl);
0136     TGeoVolume *refMod = geom->MakeBox( "refMod", Al, 0.5*_surWidth*_sclfmodulesizex, 0.5*_surLength*_sclfmodulesizey, 0.30*_sclfmodulesizez );
0137     refMod->SetLineColor( 38 );
0138     refMod->SetFillColor( 13 );
0139     TGeoVolume *curMod = geom->MakeBox( "curMod", Al, 0.5*_surWidth*_sclfmodulesizex, 0.5*_surLength*_sclfmodulesizey, 0.30*_sclfmodulesizez );
0140 
0141     if ((_yVal < 0)&&(_zVal>=0)) curMod->SetLineColor( kRed );
0142     if ((_yVal < 0)&&(_zVal<0)) curMod->SetLineColor( kGreen );
0143     if ((_yVal >= 0)&&(_zVal>=0)) curMod->SetLineColor( kBlue );
0144     if ((_yVal >= 0)&&(_zVal<0)) curMod->SetLineColor( kMagenta );
0145     refMod->SetLineColor( 14 );
0146     //curMod->SetLineColor(kBlue);
0147     //refMod->SetLineColor(kRed);
0148 
0149     const Double_t radc = 180./TMath::Pi();
0150     TGeoTranslation *tr1 = new TGeoTranslation( 0., 0., 0. );
0151     TGeoRotation *rt1 = new TGeoRotation();
0152     double rota[9];
0153     rota[0] = _surRot[0];   rota[1] = _surRot[3];   rota[2] = _surRot[6];
0154     rota[3] = _surRot[1];   rota[4] = _surRot[4];   rota[5] = _surRot[7];
0155     rota[6] = _surRot[2];   rota[7] = _surRot[5];   rota[8] = _surRot[8];
0156     rt1->SetMatrix( rota );
0157     TGeoTranslation *tr2 = new TGeoTranslation( _sclftr*_dxVal, _sclftr*_dyVal, _sclftr*_dzVal );
0158     TGeoRotation *rt2 = new TGeoRotation( "rt2", _sclfrt*_dalphaVal*radc, _sclfrt*_dbetaVal*radc, _sclfrt*_dgammaVal*radc );
0159     rt2->MultiplyBy( rt1 );
0160     TGeoCombiTrans *combi1 = new TGeoCombiTrans( *tr1, *rt1 );
0161     TGeoCombiTrans *combi2 = new TGeoCombiTrans( *tr2, *rt2 );
0162     mod->AddNode( curMod, 1, combi2 );
0163     mod->AddNode( refMod, 1, combi1 );
0164     TGeoTranslation *trG = new TGeoTranslation( _xVal - _dxVal, _yVal - _dyVal, _zVal - _dzVal);
0165     TGeoRotation *rtG = new TGeoRotation( "rtG", -1*_dalphaVal, -1*_dbetaVal, -1*_dgammaVal );
0166     TGeoCombiTrans *combi = new TGeoCombiTrans( *trG, *rtG );
0167     top->AddNode( mod, 1, combi );
0168 }
0169 
0170 bool isRightSubDet() {
0171     return (_sublevel == _subdetector1 || _sublevel == _subdetector2);
0172 }
0173 
0174 int visualizationTracker(float minZ, float maxZ, float minX, float maxX, float theta, float phi){
0175     gSystem->Load("libGeom");
0176 //++++++++++++++++++++ Set up stuff ++++++++++++++++++++//
0177     TGeoManager *geom = new TGeoManager("simple1", "Simple geometry");
0178 //--- define some materials and media
0179     TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0,0,0);
0180     TGeoMedium *Vacuum = new TGeoMedium("Vacuum",1, matVacuum);
0181 //--- make the top container volume
0182     TGeoVolume *top = geom->MakeBox("TOP", Vacuum, 500., 500., 500.);
0183 //TGeoVolume *toptop = geom->MakeBox("TOPTOP", Vacuum, 1000., 1000., 500.);
0184     geom->SetTopVolume(top);
0185 
0186     int count = 0;
0187     for (int i = 0; i < _nEntries; ++i){
0188         _inTree->GetEntry(i);
0189         if (isRightSubDet()&&(_zVal >= minZ && _zVal < maxZ)&&(_xVal >= minX && _xVal < maxX)/*&&(_rVal <= 12)&&(_rVal >=8)*/){
0190             char modName[192];
0191             sprintf(modName, "testModule%i", i);
0192             TGeoVolume* testMod = geom->MakeBox( modName, Vacuum, 90., 90., 40. );
0193             getModule( geom, top, testMod );
0194             count++;
0195         }
0196     }
0197 
0198     if(count == 0) return -1;
0199 
0200     getBeamVisuals(geom, top, minZ, maxZ);
0201 
0202 //--- close the geometry
0203     geom->CloseGeometry();
0204 // -- draw
0205     geom->SetVisLevel(4);
0206 
0207     TCanvas * c = new TCanvas();
0208     c->SetTheta(theta);
0209     c->SetPhi(phi);
0210     top->Draw();
0211 
0212 //--- putting words on canvas...
0213     bool with0T = true;
0214 
0215     //can play with these numbers
0216     double widthofeach = 0.07;
0217     double textsize = 0.05;
0218 
0219     double xmax = 2*widthofeach;
0220     if (with0T) xmax = widthofeach;
0221 
0222     TPaveText* pt = new TPaveText(0,0,xmax,1,"brNDC");
0223     pt->SetBorderSize(0);
0224     pt->SetFillStyle(0);
0225     pt->SetTextAlign(22);
0226     pt->SetTextFont(42);
0227     pt->SetTextSize(0.1);
0228     TText *text = pt->AddText(0,0,TString("#font[42]{"+_line1+"}"));
0229     text->SetTextSize(textsize);
0230     text->SetTextAngle(90);
0231     pt->Draw();
0232 
0233     TPaveText *pt2 = new TPaveText(widthofeach, 0, 2*widthofeach, 1, "brNDC");
0234     pt2->SetBorderSize(0);
0235     pt2->SetFillStyle(0);
0236     pt2->SetTextAlign(22);
0237     pt2->SetTextFont(42);
0238     pt2->SetTextSize(0.1);
0239     TText *text2 = pt2->AddText(0,0,TString("#font[42]{"+_line2+"}"));
0240     text2->SetTextSize(textsize);
0241     text2->SetTextAngle(90);
0242     pt2->Draw();
0243 
0244     TPaveText *pt3 = new TPaveText(2*widthofeach, 0, 3*widthofeach, 1, "brNDC");
0245     pt3->SetBorderSize(0);
0246     pt3->SetFillStyle(0);
0247     pt3->SetTextAlign(22);
0248     pt3->SetTextFont(42);
0249     pt3->SetTextSize(0.1);
0250     TText *text3 = pt3->AddText(0,0,TString("#font[42]{"+_line3+"}"));
0251     text3->SetTextSize(textsize);
0252     text3->SetTextAngle(90);
0253     pt3->Draw();
0254 
0255     string str = string("i") + to_string(_i) + string(".gif");
0256     c->SaveAs(TString(str));
0257     gSystem->Exec(TString("mv "+str+" images/"+str));
0258     delete c;
0259     cout << "Created image " << str << endl;
0260     return 0;
0261 }
0262 
0263 //gets minimum and maximum values of Z and Y in the specified subdetectors
0264 void getMinMax(float & minZ, float & maxZ, float & minX, float & maxX) {
0265     int i = 0;
0266     while(i < _nEntries){
0267         _inTree->GetEntry(i);
0268         if(isRightSubDet()) {
0269             _inTree->GetEntry(i);
0270             maxX = _xVal;
0271             minX = _xVal;
0272             maxZ = _zVal;
0273             minZ = _zVal;
0274             break;
0275         }
0276         ++i;
0277     }
0278     while ( i < _nEntries ) {
0279         _inTree->GetEntry(i);
0280         if (isRightSubDet()) {
0281             if( _xVal > maxX ) {
0282                 maxX = _xVal;
0283             }
0284             if( _xVal < minX ) {
0285                 minX = _xVal;
0286             }
0287             if( _zVal > maxZ ) {
0288                 maxZ = _zVal;
0289             }
0290             if( _zVal < minZ ) {
0291                 minZ = _zVal;
0292             }
0293         }
0294         ++i;
0295     }
0296     cout << minX << endl;
0297     cout << maxX << endl;
0298     cout << minZ << endl;
0299     cout << maxZ << endl;
0300 }
0301 
0302 //gets string that is a unix command that merges gifs using gifmerge (download at http://the-labs.com/GIFMerge/)
0303 string getGifMergeCommand(int start, int breakspot1, int breakspot2, int end) {
0304     string str = "";
0305     str += "./gifmerge -192,192,192 -l0 -5 ";
0306     for (int i = start; i < breakspot1; i++) {
0307         str += "images/i"+to_string(i)+".gif ";
0308     }
0309     str += "-50 ";
0310     for (int i = breakspot1; i < breakspot2; i++) {
0311         str += "images/i"+to_string(i)+".gif ";
0312     }
0313     str += "-5 ";
0314     for (int i = breakspot2; i < end-1; i++) {
0315         str += "images/i"+to_string(i)+".gif ";
0316     }
0317     str += "-100 images/i"+to_string(end-1)+".gif > "+_outputFileName+".gif";
0318     return str;
0319 }
0320 
0321 //gets string that is a unix command that merges gifs using GraphicsMagick
0322 string getConvertCommand(int start, int breakspot1, int breakspot2, int end) {
0323     string str = "";
0324     str += "gm convert -loop 0 -delay 5 ";
0325     for (int i = start; i < breakspot1; i++) {
0326         str += "images/i"+to_string(i)+".gif ";
0327     }
0328     str += "-delay 50 ";
0329     for (int i = breakspot1; i < breakspot2; i++) {
0330         str += "images/i"+to_string(i)+".gif ";
0331     }
0332     str += "-delay 5 ";
0333     for (int i = breakspot2; i < end-1; i++) {
0334         str += "images/i"+to_string(i)+".gif ";
0335     }
0336     str += "-delay 100 images/i"+to_string(end-1)+".gif   "+_outputFileName+".gif";
0337     return str;
0338 }
0339 
0340 void runVisualizer(TString input,
0341                     string output,
0342                     string line1,
0343                     string line2,
0344                     int subdetector1,
0345                     int subdetector2,
0346                     int sclftr,
0347                     int sclfrt,
0348                     float sclfmodulesizex,
0349                     float sclfmodulesizey,
0350                     float sclfmodulesizez,
0351                     float piperadius,
0352                     float pipexcoord,
0353                     float pipeycoord,
0354                     float linexcoord,
0355                     float lineycoord ) {
0356     gErrorIgnoreLevel = kError;
0357 
0358 //------Tree Read In--------
0359     TString inputFileName = input;
0360     //output file name
0361     _outputFileName = output;
0362     //title
0363     _line1 = line1;
0364     _line2 = line2;
0365     //set subdetectors to see
0366     _subdetector1 = subdetector1;
0367     _subdetector2 = subdetector2;
0368     //translation scale factor
0369     _sclftr = sclftr;
0370     //rotation scale factor
0371     _sclfrt = sclfrt;
0372     //module size scale factor
0373     _sclfmodulesizex = sclfmodulesizex;
0374     _sclfmodulesizey = sclfmodulesizey;
0375     _sclfmodulesizez = sclfmodulesizez;
0376     //beam pipe radius
0377     _piperadius = piperadius;
0378     //beam pipe xy coordinates
0379     _pipexcoord = pipexcoord;
0380     _pipeycoord = pipeycoord;
0381     //beam line xy coordinates
0382     _linexcoord = linexcoord;
0383     _lineycoord = lineycoord;
0384 
0385 
0386     sortbyz( inputFileName );
0387     TFile *fin = TFile::Open( inputFileName.ReplaceAll(".root", "_sorted.root") );
0388     _line3 = Form("Translational Scale Factor: %i",_sclftr);
0389     //++++++++++++++++++++ Read in tree ++++++++++++++++++++//
0390     _inTree = (TTree*) fin->Get("alignTree");
0391     _nEntries = _inTree->GetEntries();
0392     _inTree->SetBranchAddress( "id", &_id );
0393     _inTree->SetBranchAddress( "level", &_level);
0394     _inTree->SetBranchAddress("mid", &_mid);
0395     _inTree->SetBranchAddress("mlevel", &_mlevel);
0396     _inTree->SetBranchAddress("sublevel", &_sublevel);
0397     _inTree->SetBranchAddress("x", &_xVal);
0398     _inTree->SetBranchAddress("y", &_yVal);
0399     _inTree->SetBranchAddress("z", &_zVal);
0400     _inTree->SetBranchAddress("r", &_rVal);
0401     _inTree->SetBranchAddress("phi", &_phiVal);
0402     _inTree->SetBranchAddress("alpha", &_alphaVal);
0403     _inTree->SetBranchAddress("beta", &_betaVal);
0404     _inTree->SetBranchAddress("gamma", &_gammaVal);
0405     _inTree->SetBranchAddress("dx", &_dxVal);
0406     _inTree->SetBranchAddress("dy", &_dyVal);
0407     _inTree->SetBranchAddress("dz", &_dzVal);
0408     _inTree->SetBranchAddress("dr", &_drVal);
0409     _inTree->SetBranchAddress("dphi", &_dphiVal);
0410     _inTree->SetBranchAddress("dalpha", &_dalphaVal);
0411     _inTree->SetBranchAddress("dbeta", &_dbetaVal);
0412     _inTree->SetBranchAddress("dgamma", &_dgammaVal);
0413     _inTree->SetBranchAddress("surW", &_surWidth);
0414     _inTree->SetBranchAddress("surL", &_surLength);
0415     _inTree->SetBranchAddress("surRot", &_surRot);
0416 //_inTree->SetBranchAddress("useid", &_useid);
0417     _inTree->SetBranchAddress("identifiers", &_identifiers);
0418 
0419     float minZ, maxZ, minX, maxX;
0420     int zpos;
0421     int numincrements;
0422     getMinMax(minZ, maxZ, minX, maxX);
0423     
0424     gSystem->mkdir("images");
0425 
0426     _i = 0;
0427     for (int i = 0; i < 90; i+=1, _i++) {
0428         visualizationTracker(minZ, maxZ, minX, maxX, i, 90);
0429     }
0430     
0431     numincrements = 12;
0432     float length;
0433     int start1 = _i;
0434     length = (maxZ - minZ) / numincrements;
0435     for(int i = numincrements-1; i >= 0; i--) {
0436         zpos = minZ + i*length;
0437         if(visualizationTracker(zpos, zpos+length, minX, maxX, 90, 90) == 0) {
0438             _i++;
0439         }
0440     }
0441     
0442     int start2 = _i;
0443     for (int i = 90; i >= 0; i-=1, _i++){
0444         visualizationTracker(minZ, maxZ, (minX + maxX) / 2 - 3, (minX + maxX) / 2 + 1, i, 90);
0445     }
0446     delete fin;
0447 
0448     //gSystem->Exec(TString(getGifMergeCommand(0, start1, start2, _i)));
0449     gSystem->Exec(TString(getConvertCommand(0, start1, start2, _i)));
0450     cout << "images merged." << endl;
0451     gSystem->Exec(TString("gm convert "+_outputFileName+".gif -rotate 90 "+_outputFileName+"_rotated.gif"));
0452     cout << "images rotated." << endl;
0453 }
0454 
0455 void runVisualizer() {
0456         //------------------------------ONLY NEEDED INPUTS-------------------------------//
0457 //------Tree Read In--------
0458     TString inputFileName = "~/normal_vs_test.Comparison_commonTracker.root";
0459     //output file name
0460     string outputFileName = "animation";
0461     //title
0462     string line1 = "";
0463     string line2 = "";
0464     //set subdetectors to see
0465     int subdetector1 = 1;
0466     int subdetector2 = 2;
0467     //translation scale factor
0468     int sclftr = 50;
0469     //rotation scale factor
0470     int sclfrt = 1;
0471     //module size scale factor
0472     float sclfmodulesizex = 1;
0473     float sclfmodulesizey = 1;
0474     float sclfmodulesizez = 1;
0475     //beam pipe radius
0476     float piperadius = 2.25;
0477     //beam pipe xy coordinates
0478     float pipexcoord = 0;
0479     float pipeycoord = 0;
0480     //beam line xy coordinates
0481     float linexcoord = 0;
0482     float lineycoord = 0;
0483 //------------------------------End of ONLY NEEDED INPUTS-------------------------------//
0484     runVisualizer(inputFileName,
0485                     outputFileName,
0486                     line1,
0487                     line2,
0488                     subdetector1,
0489                     subdetector2,
0490                     sclftr,
0491                     sclfrt,
0492                     sclfmodulesizex,
0493                     sclfmodulesizey,
0494                     sclfmodulesizez,
0495                     piperadius,
0496                     pipexcoord,
0497                     pipeycoord,
0498                     linexcoord,
0499                     lineycoord );
0500 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!