Changeset 1365 in svn for trunk/external/Hector/H_OpticalElement.cc

Timestamp:

Apr 16, 2014, 3:56:14 PM (10 years ago)

Author:

Pavel Demin

Message:

switch to a more stable Hector version

File:

• trunk/external/Hector/H_OpticalElement.cc (modified) (8 diffs)

trunk/external/Hector/H_OpticalElement.cc

@@ -1 @@
-  /* * * * * * * * * * * * * * * * * * * * * * * * * * * *
- *                                                         *
-*                   --<--<--  A fast simulator --<--<--     *
-*                 / --<--<--     of particle   --<--<--     *
-*  ----HECTOR----<                                          *
-*                 \ -->-->-- transport through -->-->--     *
-*                   -->-->-- generic beamlines -->-->--     *
-*                                                           *
-* JINST 2:P09005 (2007)                                     *
-*      X Rouby, J de Favereau, K Piotrzkowski (CP3)         *
-*       http://www.fynu.ucl.ac.be/hector.html               *
-*                                                           *
-* Center for Cosmology, Particle Physics and Phenomenology  *
-*              Universite catholique de Louvain             *
-*                 Louvain-la-Neuve, Belgium                 *
- *                                                         *
-   * * * * * * * * * * * * * * * * * * * * * * * * * * * */
+/*
+---- Hector the simulator ----
+   A fast simulator of particles through generic beamlines.
+   J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be
+
+        http://www.fynu.ucl.ac.be/hector.html
+
+   Centre de Physique des Particules et de Phénoménologie (CP3)
+   Université Catholique de Louvain (UCL)
+*/
 
 /// \file H_OpticalElement.cc
@@ -25 +18 @@
 
 // ROOT #includes
-#include "TPaveLabel.h"
+//#include "TPaveLabel.h"
 
 // local #includes
-#include "H_Parameters.h"
 #include "H_TransportMatrices.h"
 #include "H_Aperture.h"
@@ -35 +27 @@
 
 /// called by the constructors
-void H_OpticalElement::init(const string& nameE, const int typeE, const double s, const double k, const double l) {
+void H_OpticalElement::init(const string nameE, const int typeE, const double s, const double k, const double l, H_Aperture* the_app) {
         // this is called by the constructors
         // must be in public section !
@@ -48 +40 @@
         element_length = l;
         type = typeE;
-        element_mat.ResizeTo(MDIM,MDIM);
-        element_mat = driftmat(l);
-
-        // do NOT use setAperture for the initialisation ! there are protections there
+        element_mat = new TMatrix(MDIM,MDIM);
+        setAperture(the_app);
         
-        if(element_length<0)  { if(VERBOSE) cout<<"<H_OpticalElement> ERROR : Interpenetration of elements !"<<endl; }
-        if(element_length==0) { if(VERBOSE) cout<<"<H_OpticalElement> WARNING : 0-length element ! (" << name << ") " << " at " << fs << endl; }
+        if(element_length<0)  { if(VERBOSE) cout<<"\t ERROR : Interpenetration of elements !"<<endl; }
+        if(element_length==0) { if(VERBOSE) cout<<"\t WARNING : 0-length element ! (" << name << ") " << " at " << fs << endl; }
         betax =0;
         betay =0;
 }
 
-H_OpticalElement::H_OpticalElement() : element_aperture(new H_Aperture()) {
-        init("",DRIFT,0.,0.,0.1);
+H_OpticalElement::H_OpticalElement() {
+        H_Aperture* the_app = new H_Aperture();
+        init("",DRIFT,0.,0.,0.1,the_app);
 }
 
-H_OpticalElement::H_OpticalElement(const string& nameE, const int typeE, const double s, const double k, const double l) : element_aperture(new H_Aperture()) {
-        init(nameE,typeE,s,k,l);
+H_OpticalElement::H_OpticalElement(const string nameE, const int typeE, const double s, const double k, const double l) {
+        H_Aperture* the_app = new H_Aperture();
+        init(nameE,typeE,s,k,l,the_app);
 }
 
-H_OpticalElement::H_OpticalElement(const string& nameE, const int typeE, const double s, const double k, const double l, H_Aperture* the_app) : element_aperture(the_app->clone()) {
-        init(nameE,typeE,s,k,l);
+H_OpticalElement::H_OpticalElement(const string nameE, const int typeE, const double s, const double k, const double l, H_Aperture* the_app) {
+        init(nameE,typeE,s,k,l,the_app);
 }
 
-H_OpticalElement::H_OpticalElement(const int typeE, const double s, const double k, const double l, H_Aperture* the_app) : element_aperture(the_app->clone()) {
-        init("",typeE,s,k,l);
+H_OpticalElement::H_OpticalElement(const int typeE, const double s, const double k, const double l, H_Aperture* the_app) {
+        init("",typeE,s,k,l,the_app);
 }
 
-H_OpticalElement::H_OpticalElement(const int typeE, const double s, const double k, const double l) : element_aperture(new H_Aperture()) {
-        init("",typeE,s,k,l);
+H_OpticalElement::H_OpticalElement(const int typeE, const double s, const double k, const double l) {
+        H_Aperture* the_app = new H_Aperture();
+        init("",typeE,s,k,l,the_app);
 }
 
@@ -92 +85 @@
         name = el.name;
         typestring = el.typestring;
-        element_mat.ResizeTo(MDIM,MDIM);
-        element_mat = el.element_mat;
-        element_aperture = el.element_aperture->clone();
+        element_mat = new TMatrix(*(el.element_mat));
+        element_aperture = new H_Aperture(*(el.element_aperture));
 }
 
@@ -104 +96 @@
         xpos = el.xpos;
         ypos = el.ypos;
-        txpos = el.txpos;
-        typos = el.typos;
+                txpos = el.txpos;
+                typos = el.typos;
         betax = el.betax;
         betay = el.betay;
@@ -111 +103 @@
         name = el.name;
         typestring = el.typestring;
-        element_mat.ResizeTo(MDIM,MDIM);
-        element_mat = el.element_mat;
-        element_aperture = el.element_aperture->clone();
+        delete element_mat;
+        delete element_aperture;
+        element_mat = new TMatrix(*(el.element_mat));
+        element_aperture = new H_Aperture(*(el.element_aperture));
         return *this;
 }
 
-void H_OpticalElement::setAperture(const H_Aperture* ap) {
-        // do NOT use setAperture in your constructor, as element_aperture is not initialized
-        // this function do not take into account ap if ap=0
-        // do nothing if element_mat = ap
-                if (!ap) {cout << "<H_OpticalElement> Trying to set an empty pointer for the aperture ! Nothing done.\n"; return;}
-                if (element_aperture != ap) {
-                        delete element_aperture;
-                        element_aperture = ap->clone();
-                }
+void H_OpticalElement::setAperture(H_Aperture * ap) {
+//      element_aperture = const_cast<H_Aperture*>(ap);
+        element_aperture = ap;
         return;
 }
-/*
-void H_OpticalElement::setAperture(H_Aperture* ap) {
-        // do NOT use setAperture in your constructor, as element_aperture is not initialized
-        // this function do not take into account ap if ap=0
-        // do nothing if element_mat = ap
-                if (!ap) {cout << "<H_OpticalElement> Trying to set an empty pointer for the aperture ! Nothing done.\n"; return;}
-                if (element_aperture != ap) {
-                        delete element_aperture;
-                        element_aperture = ap; //->clone();
-                }
-        return;
-}
-*/
 
-std::ostream& operator<< (std::ostream& os, const H_OpticalElement& el) {
-        os << el.typestring << el.name << "\t at s = " << el.fs << "\t length = "<< el.element_length;
-        if(el.fk!=0) os <<"\t strength = " << el.fk;
-        if(el.element_aperture->getType()!=NONE) {
-                os << *(el.element_aperture) << endl;
-        }
-        os<<endl;
-        if(el.element_length<0 && VERBOSE) 
-                os <<"<H_OpticalElement> ERROR : Interpenetration of elements !"<<endl; 
-        else if(el.element_length==0 && VERBOSE) 
-                os <<"<H_OpticalElement> WARNING : 0-length "<< el.typestring << " !" << endl;
-  return os;
+void H_OpticalElement::printProperties() const {
+                cout << typestring;
+                cout << name;
+                cout <<"\t at s = " << fs;
+                cout <<"\t length = "<< element_length;
+                if(fk!=0) cout <<"\t strength = " << fk;
+                if(element_aperture->getType()!=NONE) { 
+                        cout <<"\t aperture type = " << element_aperture->getTypeString();
+                        element_aperture->printProperties();                    
+                }
+        
+                cout<<endl;
+                if(element_length<0)  { if(VERBOSE) cout<<"\t ERROR : Interpenetration of elements !"<<endl; }
+                if(element_length==0) { if(VERBOSE) cout<<"\t WARNING : 0-length "<< typestring << " !" << endl; }
+
+ return;
 }
 
 void H_OpticalElement::showMatrix() const {
-        printMatrix(element_mat);
+                printMatrix(element_mat);
         return;
 }
 
 void H_OpticalElement::drawAperture() const {
-        element_aperture->draw(1);
+        element_aperture->draw();
         return;
 }
 
-TMatrix H_OpticalElement::getMatrix() {
-        setMatrix(0,MP,1);
-        return element_mat;
+TMatrix H_OpticalElement::getMatrix() const {
+        return *element_mat;
 }
 
-TMatrix H_OpticalElement::getMatrix(const float eloss, const float p_mass, const float p_charge) {
+TMatrix H_OpticalElement::getMatrix(const float eloss, const float p_mass, const float p_charge) const {
         setMatrix(eloss,p_mass,p_charge);
-        return element_mat;
+        return *element_mat;
 }
 
 void H_OpticalElement::draw(const float meight, const float height) const{
-        /// @param meight is the minimal extend of the graph
+/*      /// @param meight is the minimal extend of the graph
         /// @param height is the maximal extend of the graph
         float x1 = getS();
@@ -188 +165 @@
         cur_box->SetFillColor((int)getType());
         cur_box->Draw();
+*/
 }
 
-TVectorD H_OpticalElement::getHitPosition(const TVectorD& init_pos, const double energy_loss, const double mp, const double qp) {
-        if(!element_length) {
-                // cout<<"O-length element ("<<getName()<<"), should not appear here !"<<endl;
-                return init_pos;
-        }
-        // some declarations
-        bool inside = false;
-        double vec[MDIM] = {init_pos[INDEX_X]/URAD,  tan(init_pos[INDEX_TX]/URAD),
-                            init_pos[INDEX_Y]/URAD,  tan(init_pos[INDEX_TY]/URAD), 
-                            -energy_loss,            1};
-        TMatrixD mat_init(1,MDIM,vec);
-        TMatrixD mat_min(1,MDIM,vec);
-        TMatrixD mat_max(1,MDIM,vec);
-        TMatrixD mat_stop(1,MDIM,vec);
-        H_OpticalElement* temp_el = clone();
-        // initialasing boundaries
-        double min_pos = 0;
-        double max_pos = element_length/2.;
-        double max_old = element_length/2.;
-        // number of iterations 
-        // (idea : fix precision instead of number of iterations + add security)
-        // (idea : interpolate between max and min and give the error)
-        const int N = 10;
-        // starting search loop
-        for(int i = 0; i < N; i++) {
-                // fixing position to be investigated 
-                temp_el->setLength(max_pos);
-                // initialising the vector at the initial vector + possible shift/tilt
-                mat_max[0][0] = mat_init[0][0] - temp_el->getX();
-                mat_max[0][1] = mat_init[0][1] - tan(temp_el->getTX());
-                mat_max[0][2] = mat_init[0][2] - temp_el->getY();
-                mat_max[0][3] = mat_init[0][3] - tan(temp_el->getTY());
-                // propagating
-                mat_max *= temp_el->getMatrix(energy_loss,mp,qp);
-                // compensating for the previously stated shifts/tilts
-                mat_max[0][0] = mat_max[0][0] + temp_el->getX();
-                mat_max[0][1] = mat_max[0][1] + tan(temp_el->getTX());
-                mat_max[0][2] = mat_max[0][2] + temp_el->getY();
-                mat_max[0][3] = mat_max[0][3] + tan(temp_el->getTY());
-                // fixing new boundaries 
-                if(temp_el->isInside(mat_max.GetMatrixArray()[0]*URAD,mat_max.GetMatrixArray()[2]*URAD)) {
-                        max_old = max_pos;
-                        max_pos = max_pos + (max_pos - min_pos)/2.;
-                        min_pos = max_old;
-                        inside = true;
-                } else {
-                        max_pos = min_pos + (max_pos - min_pos)/2.;
-                        inside = false;
-                }
-                // end of loop
-        }
-        // if it passes at the last iteration, choosing the other range²
-        if(inside) min_pos = max_old;
-        // here the interpolation method : now we are sure that the intercept is between min_pos and max_pos
-        // getting vector at min_pos (for first boundary) :
-        bool precision_estimate = false;
-        if(precision_estimate) {
-                temp_el->setLength(min_pos);
-                mat_min[0][0] = mat_init[0][0] - temp_el->getX();
-                mat_min[0][1] = mat_init[0][1] - tan(temp_el->getTX());
-                mat_min[0][2] = mat_init[0][2] - temp_el->getY();
-                mat_min[0][3] = mat_init[0][3] - tan(temp_el->getTY());
-                mat_min *= temp_el->getMatrix(energy_loss,mp,qp);
-                mat_min[0][0] = mat_min[0][0] + temp_el->getX();
-                mat_min[0][1] = mat_min[0][1] + tan(temp_el->getTX());
-                mat_min[0][2] = mat_min[0][2] + temp_el->getY();
-                mat_min[0][3] = mat_min[0][3] + tan(temp_el->getTY());
-                mat_min[0][4] = min_pos + init_pos[4];
-                // getting vector at max_pos (for second boundary) :
-                temp_el->setLength(max_pos);
-                mat_max[0][0] = mat_init[0][0] - temp_el->getX();
-                mat_max[0][1] = mat_init[0][1] - tan(temp_el->getTX());
-                mat_max[0][2] = mat_init[0][2] - temp_el->getY();
-                mat_max[0][3] = mat_init[0][3] - tan(temp_el->getTY());
-                mat_max *= temp_el->getMatrix(energy_loss,mp,qp);
-                mat_max[0][0] = mat_max[0][0] + temp_el->getX();
-                mat_max[0][1] = mat_max[0][1] + tan(temp_el->getTX());
-                mat_max[0][2] = mat_max[0][2] + temp_el->getY();
-                mat_max[0][3] = mat_max[0][3] + tan(temp_el->getTY());
-                mat_max[0][4] = max_pos + init_pos[4];
-        }
-        // getting vector in the middle (for estimate) :
-        temp_el->setLength((max_pos+min_pos)/2.);
-        mat_stop[0][0] = mat_init[0][0] - temp_el->getX();
-        mat_stop[0][1] = mat_init[0][1] - tan(temp_el->getTX());
-        mat_stop[0][2] = mat_init[0][2] - temp_el->getY();
-        mat_stop[0][3] = mat_init[0][3] - tan(temp_el->getTY());
-        mat_stop *= temp_el->getMatrix(energy_loss,mp,qp);
-        mat_stop[0][0] = mat_stop[0][0] + temp_el->getX();
-        mat_stop[0][1] = mat_stop[0][1] + tan(temp_el->getTX());
-        mat_stop[0][2] = mat_stop[0][2] + temp_el->getY();
-        mat_stop[0][3] = mat_stop[0][3] + tan(temp_el->getTY());
-        mat_stop[0][4] = (max_pos+min_pos)/2. + init_pos[4];
-
-        double xys[LENGTH_VEC];
-        xys[INDEX_X]=  mat_stop[0][0]*URAD;
-        xys[INDEX_TX]= atan(mat_stop[0][1])*URAD;
-        xys[INDEX_Y]=  mat_stop[0][2]*URAD;
-        xys[INDEX_TY]= atan(mat_stop[0][3])*URAD;
-        xys[INDEX_S]=  mat_stop[0][4] ;
-        TVectorD temp_vec(LENGTH_VEC,xys);
-
-        if(precision_estimate) {
-                cout<<"--- Results and precision estimates ---"<<endl;
-                cout<<"\t Stopping element : "<<getName()<<endl;
-                cout<<"\t hit point s  : "<<mat_stop[0][4]<<" m +- "<<(mat_max[0][4]-mat_stop[0][4])*1000.<<" mm"<<endl;
-                cout<<"\t hit point x  : "<<mat_stop[0][0]*URAD;
-                cout<<" + "<<fabs(((mat_min[0][0]<mat_max[0][0])?(mat_max[0][0]-mat_stop[0][0])*URAD:(mat_min[0][0]-mat_stop[0][0])*URAD));
-                cout<<" - "<<fabs(((mat_min[0][0]<mat_max[0][0])?(mat_stop[0][0]-mat_min[0][0])*URAD:(mat_stop[0][0]-mat_max[0][0])*URAD))<<" µm"<<endl;
-                cout<<"\t hit point y  : "<<mat_stop[0][2]*URAD;
-                cout<<" + "<<fabs(((mat_min[0][0]<mat_max[0][2])?(mat_max[0][2]-mat_stop[0][2])*URAD:(mat_min[0][2]-mat_stop[0][2])*URAD));
-                cout<<" - "<<fabs(((mat_min[0][0]<mat_max[0][2])?(mat_stop[0][2]-mat_min[0][2])*URAD:(mat_stop[0][2]-mat_max[0][2])*URAD))<<" µm"<<endl;
-                cout<<"\t hit point tx : "<<mat_stop[0][1]*URAD;
-                cout<<" + "<<fabs(((mat_min[0][0]<mat_max[0][1])?(mat_max[0][1]-mat_stop[0][1])*URAD:(mat_min[0][1]-mat_stop[0][1])*URAD));
-                cout<<" - "<<fabs(((mat_min[0][0]<mat_max[0][1])?(mat_stop[0][1]-mat_min[0][1])*URAD:(mat_stop[0][1]-mat_max[0][1])*URAD))<<" µrad"<<endl;
-                cout<<"\t hit point ty : "<<mat_stop[0][3]*URAD;
-                cout<<" + "<<fabs(((mat_min[0][0]<mat_max[0][3])?(mat_max[0][3]-mat_stop[0][3])*URAD:(mat_min[0][3]-mat_stop[0][3])*URAD));
-                cout<<" - "<<fabs(((mat_min[0][0]<mat_max[0][3])?(mat_stop[0][3]-mat_min[0][3])*URAD:(mat_stop[0][3]-mat_max[0][3])*URAD))<<" µrad"<<endl;
-
-        }
-
-        delete temp_el;
-        return temp_vec;
-}