drich-dev/src/findFocalPlane.cpp

0001 // After using debug optics mode 4/simulate.py test 14, calculate
0002 // points of closest approach of reflected parallel photon beams,
0003 // to get approximate focal region.
0004
0005 // Produces text file of focal point positions and directions.
0006 // To draw these points alongside dRICH geometry, set
0007 // DRICH_debug_optics mode 5 and DRICH_FP_file to the generated
0008 // text file.
0009 #include <cstdlib>
0010 #include <iostream>
0011
0012 // ROOT
0013 #include "TSystem.h"
0014 #include "TStyle.h"
0015 #include "TRegexp.h"
0016 #include "TCanvas.h"
0017 #include "TApplication.h"
0018 #include "TBox.h"
0019 #include "ROOT/RDataFrame.hxx"
0020 #include "TVector3.h"
0021 #include "TMatrix.h"
0022 #include "TMatrixD.h"
0023 #include "TMatrixDSparse.h"
0024 // edm4hep
0025 #include "edm4hep/MCParticleCollection.h"
0026 #include "edm4hep/SimTrackerHitCollection.h"
0027
0028 using std::cout;
0029 using std::cerr;
0030 using std::endl;
0031
0032 using namespace ROOT;
0033 using namespace ROOT::VecOps;
0034 using namespace edm4hep;
0035
0036 TCanvas *CreateCanvas(TString name, Bool_t logx=0, Bool_t logy=0, Bool_t logz=0);
0037 std::vector<TVector3> d, e, fp, dirout;
0038 int nphotons = 0;
0039 int nAcc = 0;
0040 TVector3 focalPoint(std::vector<TVector3> dvec, std::vector<TVector3> endvec){
0041   TMatrixD a(3,3), b(3,1);
0042   TArrayI row(3),col(3);
0043
0044   for(int i = 0; i < 3; i++) row[i] = col[i] = i;
0045   TArrayD idarr(3); idarr.Reset(1.);
0046   TMatrixDSparse identity(3,3);
0047   identity.SetMatrixArray(3,row.GetArray(),col.GetArray(),idarr.GetArray());
0048
0049   for(int i = 0; i < dvec.size(); i++){
0050     double darr[] = {dvec[i].X(),dvec[i].Y(),dvec[i].Z()};
0051     double earr[] = {endvec[i].X(),endvec[i].Y(),endvec[i].Z()};
0052
0053     auto matd = TMatrixD(3,1,darr);
0054     auto mate = TMatrixD(3,1,earr);
0055     auto matdT = TMatrixD(matd);
0056     matdT.T();
0057
0058     a += (identity - matd*matdT);
0059     b += (identity - matd*matdT)*mate;
0060
0061   }
0062   auto c = (a.Invert()*b);
0063   TVector3 x(TMatrixDRow(c,0)(0),TMatrixDRow(c,1)(0),TMatrixDRow(c,2)(0));
0064   return x;
0065 };
0066
0067 TVector3 avgDir(std::vector<TVector3> dvec){
0068   double xavg = 0;
0069   double yavg = 0;
0070   double zavg = 0;
0071   double n = 0;
0072   for(int i = 0; i < dvec.size(); i++){
0073     xavg += dvec[i].X();
0074     yavg += dvec[i].Y();
0075     zavg += dvec[i].Z();
0076     n+=1.0;
0077   }
0078
0079   TVector3 outdir;
0080   if(n > 0) outdir.SetXYZ(xavg/n, yavg/n, zavg/n);
0081   return outdir;
0082 };
0083
0084 int main(int argc, char** argv) {
0085   // setup
0086   TString infileN="out/sim.edm4hep.root";
0087   if(argc>1) infileN = TString(argv[1]);
0088
0089   RDataFrame dfIn("events",infileN.Data());
0090   gStyle->SetOptStat(0);
0091
0092   /* lambdas
0093    * - most of these transform an `RVec<T1>` to an `RVec<T2>` using `VecOps::Map` or `VecOps::Filter`
0094    * - see NPdet/src/dd4pod/dd4hep.yaml for POD syntax
0095    */
0096   // calculate number of hits
0097   auto numHits = [](RVec<SimTrackerHitData> hits) { return hits.size(); };
0098   // filter for thrown particles
0099   auto isThrown = [](RVec<MCParticleData> parts){
0100     return Filter(parts,[](auto p){
0101         return p.generatorStatus==1;
0102         });
0103   };
0104   // get positions for each hit (units=cm)
0105   auto hitPos = [](RVec<SimTrackerHitData> hits){ return Map(hits,[](auto h){ return h.position; }); };
0106   auto hitPosX = [](RVec<Vector3d> v){ return Map(v,[](auto p){ return p.x; }); };
0107   auto hitPosY = [](RVec<Vector3d> v){ return Map(v,[](auto p){ return p.y; }); };
0108   auto hitPosZ = [](RVec<Vector3d> v){ return Map(v,[](auto p){ return p.z; }); };
0109
0110   auto findReflected = [](RVec<MCParticleData> photons){ // want to skip non-reflected photons
0111     return Map( photons,[](auto p){
0112       auto imom = p.momentum;
0113       auto fmom = p.momentumAtEndpoint;
0114       return !(imom==fmom);
0115     });
0116   };
0117   auto dirVec = [](RVec<MCParticleData> parts){
0118     return Map(parts,[](auto p){
0119       auto dMom = TVector3(p.momentumAtEndpoint.x,p.momentumAtEndpoint.y,p.momentumAtEndpoint.z);
0120       return dMom.Unit();
0121       //return dMom;
0122     });
0123   };
0124   auto endVec = [](RVec<MCParticleData> parts){
0125     return Map(parts,[](auto p){
0126       return TVector3(p.endpoint.x/10., p.endpoint.y/10., p.endpoint.z/10.);
0127     });
0128   };
0129
0130   // transformations
0131   auto df1 = dfIn
0132     .Define("reflectedPhotons",findReflected,{"MCParticles"})
0133     .Define("direcVec",dirVec,{"MCParticles"})
0134     .Define("endVec",endVec,{"MCParticles"})
0135     ;
0136
0137   df1.Foreach(
0138           [](RVec<bool> refl, RVec<TVector3> dir, RVec<TVector3> end){
0139         if( refl[0]  ){  d.push_back(dir[0]); e.push_back(end[0]); nAcc++;}
0140         nphotons++;
0141         if(nphotons==50){
0142           if(nAcc>0){
0143             fp.push_back(focalPoint(d,e));
0144             dirout.push_back(avgDir(d));
0145           }
0146           nAcc = 0;
0147           nphotons = 0;
0148           d.clear();
0149           e.clear();
0150         }
0151           },
0152           {"reflectedPhotons","direcVec","endVec"}
0153           );
0154   auto dfFinal = df1;
0155   FILE * outtxt = fopen("focalPoints.txt","w");
0156   for(int i = 0 ; i < fp.size(); i++){
0157     if( std::abs(fp[i].X()) < 1000 && std::abs(fp[i].Y()) < 1000 && std::abs(fp[i].Z()) < 1000){
0158       fprintf(outtxt, "%lf %lf %lf %lf %lf %lf \n", fp[i].X(), fp[i].Y(), fp[i].Z(), dirout[i].X(), dirout[i].Y(), dirout[i].Z());
0159     }
0160   }
0161   fclose(outtxt);
0162
0163   return 0;
0164 };
0165
0166
0167 TCanvas *CreateCanvas(TString name, Bool_t logx, Bool_t logy, Bool_t logz) {
0168   TCanvas *c = new TCanvas("canv_"+name,"canv_"+name,800,600);
0169   c->SetGrid(1,1);
0170   if(logx) c->SetLogx(1);
0171   if(logy) c->SetLogy(1);
0172   if(logz) c->SetLogz(1);
0173   return c;
0174 };
0175