obsolete/ak/xray.C

0001
0002 #define _DISTANCE_ 1185.5
0003 #define _X0_          0.0
0004 #define _Y0_          0.0
0005 //#define _SIZE_      250.0
0006 #define _XYDIM_       500
0007 //#define _SIZE_     400.0
0008 //#define _XYDIM_       300
0009 #define _SIZE_     1300.0
0010 //#define _XYDIM_       300
0011 #define _BWIDTH_ (_SIZE_/_XYDIM_)
0012
0013 // Help stepping not get stuck; fine with the pfRICH; may not work well for other detectors;
0014 #define _MAX_STEP_COUNT_              50
0015 // In [cm]; suppress fake IP file air;
0016 #define _MAX_POSSIBLE_STEP_LENGTH_  1000
0017 // Suppress vacuum entries;
0018 #define _MAX_ACCOUNTABLE_RAD_LENGTH_ 1E6
0019 // Not really interested in anything beyond this length along the scan line from the IP;
0020 #define _MAX_SCAN_DEPTH_             500
0021
0022 void xray(const char *fname)
0023 {
0024   TGeoManager::Import(fname);
0025
0026   auto hrlen = new TH2D("hrlen", "pfRICH radiation length scan [%]",
0027             _XYDIM_, _X0_ - _SIZE_/2, _X0_ + _SIZE_/2,
0028             _XYDIM_, _Y0_ - _SIZE_/2, _Y0_ + _SIZE_/2);
0029   hrlen->GetXaxis()->SetTitle("pfRICH vessel front side X, [mm]");
0030   hrlen->GetYaxis()->SetTitle("pfRICH vessel front side Y, [mm]");
0031   hrlen->GetXaxis()->SetTitleOffset(1.20);
0032   hrlen->GetYaxis()->SetTitleOffset(1.40);
0033
0034   double z = -_DISTANCE_;
0035   for(unsigned i = 0; i < _XYDIM_; i++) {
0036     double x = _X0_ - _SIZE_/2 + _BWIDTH_*(i + 0.5);
0037
0038     for(unsigned j = 0; j < _XYDIM_; j++) {
0039       double y = _Y0_ - _SIZE_/2 + _BWIDTH_*(j + 0.5);
0040       double r = sqrt(x*x + y*y + z*z);
0041       double xx[3] = {0.0, 0.0, 0.0}, nn[3] = {x, y, z};
0042       for(auto &coord: nn)
0043     coord /= r;
0044
0045       gGeoManager->SetCurrentPoint    (xx);
0046       gGeoManager->SetCurrentDirection(nn);
0047
0048       unsigned counter = 0;
0049       double accu = 0.0;//, length = 0.0;
0050       for(TGeoNode *node = gGeoManager->GetCurrentNode(); ; ) {
0051     TGeoMaterial *material = node->GetVolume()->GetMaterial();
0052     double radlen = material->GetRadLen();
0053
0054     auto xx = gGeoManager->GetCurrentPoint();
0055     //printf("%f %f %f\n", xx[0], xx[1], xx[2]);
0056     gGeoManager->FindNextBoundary();
0057     double thickness = gGeoManager->GetStep(), distance = sqrt(xx[0]*xx[0] + xx[1]*xx[1] + xx[2]*xx[2]);;
0058     //length += thickness;
0059     if (distance < _MAX_SCAN_DEPTH_ && thickness < _MAX_POSSIBLE_STEP_LENGTH_ &&
0060         radlen < _MAX_ACCOUNTABLE_RAD_LENGTH_) {
0061       //if (100 * thickness / radlen > 10)
0062       //printf("%7.2f [cm] %10.2f [cm] (D = %7.2f [cm]) --> %7.3f%%\n",
0063       //       thickness, radlen, distance, 100 * thickness / radlen);
0064       if (thickness && radlen) accu += 100 * thickness / radlen;
0065     } //if
0066
0067     // Switch to next node along {xx, nn[]} 3D line;
0068     node = gGeoManager->Step();
0069     assert(gGeoManager->IsEntering());
0070
0071     // Once out of volume, break;
0072     if (gGeoManager->IsOutside()) break;
0073
0074     if (counter++ >= _MAX_STEP_COUNT_) break;
0075       } //for inf
0076
0077       hrlen->SetBinContent(i+1, j+1, accu);
0078     } //for j
0079   } //for i
0080
0081   gStyle->SetOptStat(0);
0082   auto cv = new TCanvas("cv", "", 800, 800);
0083   hrlen->SetMinimum( 0);
0084   //hrlen->SetMaximum(47);
0085   //hrlen->SetMaximum(35);
0086   hrlen->Draw("COLZ");
0087
0088   for(unsigned ir=0; ir<2; ir++) {
0089     double r = _DISTANCE_*tan(2*atan(exp(ir ? -3.5 : -4.0)));
0090     TEllipse *el = new TEllipse(0, 0, r, r);
0091     el->SetFillStyle(kNone);
0092     el->SetLineColor(ir ? kGreen : kRed);
0093     el->SetLineStyle(9);//kDashed);
0094     el->SetLineWidth(4);
0095     //+el->Draw("SAME");
0096   } //for ir
0097 } // xray()