Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 10:40:32 

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