Changeset 1365 in svn for trunk/external
- Timestamp:
- Apr 16, 2014, 3:56:14 PM (11 years ago)
- Location:
- trunk/external/Hector
- Files:
-
- 54 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/external/Hector/H_AbstractBeamLine.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 11 18 12 19 13 /// \file H_AbstractBeamLine.cc … … 27 21 28 22 // ROOT #includes 29 #include "TPaveLabel.h"30 #include "TLine.h"31 #include "TGaxis.h"32 #include "TLegend.h"33 #include "TF1.h"34 #include "TROOT.h"23 //#include "TPaveLabel.h" 24 //#include "TLine.h" 25 //#include "TGaxis.h" 26 //#include "TLegend.h" 27 //#include "TF1.h" 28 //#include "TROOT.h" 35 29 36 30 // local #includes … … 39 33 #include "H_Drift.h" 40 34 #include "H_AbstractBeamLine.h" 35 #include "H_BeamParticle.h" 41 36 #include "H_RomanPot.h" 42 37 using namespace std; 43 38 44 39 void H_AbstractBeamLine::init(const float length) { 45 beam_mat.ResizeTo(MDIM,MDIM); 46 beam_mat = driftmat(length); 40 beam_mat = new TMatrix(MDIM,MDIM); 47 41 beam_length = length; 48 42 H_Drift * drift0 = new H_Drift("Drift0",0.,length); … … 51 45 } 52 46 53 H_AbstractBeamLine::H_AbstractBeamLine(const H_AbstractBeamLine& beamline) : 54 matrices(beamline.matrices) { 55 //elements = beamline.elements; //<-- bad ! the new vector contains the same pointers as the previous one 56 cloneElements(beamline); 57 beam_mat.ResizeTo(MDIM,MDIM); 58 beam_mat = beamline.beam_mat; 47 H_AbstractBeamLine::H_AbstractBeamLine(const H_AbstractBeamLine& beamline) { 48 elements = beamline.elements; 49 matrices = beamline.matrices; 50 beam_mat = new TMatrix(*(beamline.beam_mat)); 59 51 beam_length = beamline.beam_length; 60 52 } … … 62 54 H_AbstractBeamLine& H_AbstractBeamLine::operator=(const H_AbstractBeamLine& beamline) { 63 55 if(this== &beamline) return *this; 64 //elements = beamline.elements; //<-- bad ! the new vector contains the same pointers as the previous one 65 cloneElements(beamline); 56 elements = beamline.elements; 66 57 matrices = beamline.matrices; 67 beam_mat = beamline.beam_mat;58 beam_mat = new TMatrix(*(beamline.beam_mat)); 68 59 beam_length = beamline.beam_length; 69 60 return *this; 70 }71 72 H_AbstractBeamLine* H_AbstractBeamLine::clone() const {73 H_AbstractBeamLine* temp_beam = new H_AbstractBeamLine(beam_length);74 vector<H_OpticalElement*>::const_iterator element_i;75 for (element_i = elements.begin(); element_i<elements.end(); element_i++) {76 if((*element_i)->getType()!=DRIFT) {77 H_OpticalElement* temp_el = (*element_i)->clone();78 temp_beam->add(temp_el);79 }80 }81 temp_beam->beam_mat = beam_mat;82 temp_beam->matrices = matrices;83 return temp_beam;84 61 } 85 62 … … 91 68 elements.clear(); 92 69 matrices.clear(); 93 } 94 95 void H_AbstractBeamLine::cloneElements(const H_AbstractBeamLine& beam) { 96 vector<H_OpticalElement*>::const_iterator element_i; 97 for (element_i = beam.elements.begin(); element_i< beam.elements.end(); element_i++) { 98 H_OpticalElement* temp_el = (*element_i)->clone(); 99 elements.push_back(temp_el); 100 } 70 delete beam_mat; 101 71 } 102 72 … … 110 80 if (a > beam_length) { 111 81 beam_length = a; 112 if(VERBOSE) cout<<" <H_AbstractBeamLine>WARNING : element ("<< newElement->getName()<<") too far away. The beam length has been extended to "<< beam_length << ". "<<endl;82 if(VERBOSE) cout<<"WARNING : element ("<< newElement->getName()<<") too far away. The beam length has been extended to "<< beam_length << ". "<<endl; 113 83 } 114 84 calcSequence(); 115 85 calcMatrix(); 116 return; 117 } 118 119 const TMatrix H_AbstractBeamLine::getBeamMatrix() const { 120 return beam_mat; 121 } 122 123 const TMatrix H_AbstractBeamLine::getBeamMatrix(const float eloss,const float p_mass, const float p_charge) { 86 } 87 88 void H_AbstractBeamLine::add(H_OpticalElement & newElement) { 89 /// @param newElement is added to the beamline 90 // H_OpticalElement * el = new H_OpticalElement(newElement); 91 // elements.push_back(el); 92 elements.push_back(&newElement); 93 float a = newElement.getS()+newElement.getLength(); 94 if (a > beam_length) { 95 beam_length = a; 96 if(VERBOSE) cout<<"WARNING : element ("<< newElement.getName()<<") too far away. The beam length has been extended to "<< beam_length << ". "<<endl; 97 } 98 calcSequence(); 99 calcMatrix(); 100 } 101 102 const TMatrix * H_AbstractBeamLine::getBeamMatrix() const { 103 TMatrix * mat = new TMatrix(*beam_mat); 104 return mat; 105 } 106 107 const TMatrix * H_AbstractBeamLine::getBeamMatrix(const float eloss,const float p_mass, const float p_charge) { 124 108 125 109 vector<H_OpticalElement*>::iterator element_i; … … 134 118 calc_mat *= (*element_i)->getMatrix(eloss,p_mass,p_charge); 135 119 } 136 return calc_mat; 137 } 138 139 const TMatrix H_AbstractBeamLine::getPartialMatrix(const string& elname, const float eloss, const float p_mass, const float p_charge) { 120 const TMatrix* bmat = new TMatrix(calc_mat); 121 return bmat; 122 } 123 124 const TMatrix * H_AbstractBeamLine::getPartialMatrix(const string elname, const float eloss, const float p_mass, const float p_charge) { 140 125 141 126 vector<H_OpticalElement*>::iterator element_i; … … 147 132 calc_mat *= (*element_i)->getMatrix(eloss,p_mass,p_charge); 148 133 if(elname==(*element_i)->getName()) { 149 return calc_mat; 150 } 151 } 152 cout<<"<H_AbstractBeamLine> Element "<<elname<<" desn't exist. Returning full beam matrix"<<endl; 153 return calc_mat; 154 } 155 156 const TMatrix H_AbstractBeamLine::getPartialMatrix(const unsigned int element_position) const { 134 const TMatrix* bmat = new TMatrix(calc_mat); 135 return bmat; 136 } 137 } 138 cout<<"Element "<<elname<<" desn't exist. Returning full beam matrix"<<endl; 139 const TMatrix* bmat = new TMatrix(calc_mat); 140 return bmat; 141 } 142 143 const TMatrix * H_AbstractBeamLine::getPartialMatrix(const unsigned int element_position) const { 157 144 //const int N = (element_position<0)?0:(( (element_position)>elements.size()-1)?elements.size()-1:element_position); 158 145 const int N = (element_position>elements.size()-1)?elements.size()-1:element_position; 159 return *(matrices.begin()+N); // //for optimization of the code :same as return &matrices[N];160 } 161 162 const TMatrix H_AbstractBeamLine::getPartialMatrix(const H_OpticalElement * element) const{146 return &(*(matrices.begin()+N)); // //for optimization of the code :same as return &matrices[N]; 147 } 148 149 const TMatrix * H_AbstractBeamLine::getPartialMatrix(const H_OpticalElement * element) const{ 163 150 // returns the transport matrix to transport until the end of the specified element 164 151 // !!! 2 elements should never have the same name in "elements" !!! … … 166 153 vector<H_OpticalElement*>::const_iterator element_i; 167 154 vector<TMatrix>::const_iterator matrix_i; 168 TMatrix calc_mat(MDIM,MDIM);155 TMatrix * calc_mat = new TMatrix(MDIM,MDIM); 169 156 170 157 // parses the list of optical elements and find the searched one … … 172 159 if(element->getName() == (*element_i)->getName()) { 173 160 // element has been found 174 calc_mat = *matrix_i;161 calc_mat = const_cast<TMatrix*>( &(*matrix_i)); 175 162 } 176 163 } 177 return calc_mat;164 return (const TMatrix*) calc_mat; 178 165 } 179 166 … … 183 170 } 184 171 185 H_OpticalElement * H_AbstractBeamLine::getElement(const unsigned int element_position) const {172 const H_OpticalElement * H_AbstractBeamLine::getElement(const unsigned int element_position) const { 186 173 const unsigned int N = (element_position>elements.size())?elements.size():element_position; 187 174 return *(elements.begin()+N);//for optimization of the code :same as return &elements[N]; … … 189 176 190 177 191 H_OpticalElement * H_AbstractBeamLine::getElement(const string &el_name) {178 H_OpticalElement * H_AbstractBeamLine::getElement(const string el_name) { 192 179 for(unsigned int i=0; i < elements.size(); i++) { 193 180 if( (*(elements.begin()+i))->getName() == el_name ) 194 181 return *(elements.begin()+i); 195 182 } // if found -> return ; else : not found at all ! 196 cout<<" <H_AbstractBeamLine>Element "<<el_name<<" not found"<<endl;183 cout<<"Element "<<el_name<<" not found"<<endl; 197 184 return *(elements.begin()+1); 198 185 } 199 186 200 H_OpticalElement * H_AbstractBeamLine::getElement(const string&el_name) const {187 const H_OpticalElement * H_AbstractBeamLine::getElement(const string el_name) const { 201 188 for(unsigned int i=0; i < elements.size(); i++) { 202 189 if( (*(elements.begin()+i))->getName() == el_name) 203 190 return *(elements.begin()+i); 204 191 } // if found -> return ; else : not found at all ! 205 cout<<" <H_AbstractBeamLine>Element "<<el_name<<" not found"<<endl;192 cout<<"Element "<<el_name<<" not found"<<endl; 206 193 return *(elements.begin()+1); 207 194 } 208 195 209 std::ostream& operator<< (std::ostream& os, const H_AbstractBeamLine& be) { 210 vector<H_OpticalElement*>::const_iterator element_i; 211 os << "Beamline content" << endl; 212 for (element_i = be.elements.begin(); element_i < be.elements.end(); element_i++) { 213 os << (int)(element_i - be.elements.begin()) << "\t" << (*element_i)->getName() << "\t" << (*element_i)->getS() << endl; 214 } 215 return os; 216 } 217 218 void H_AbstractBeamLine::printProperties() const { cout << *this; return; } 196 void H_AbstractBeamLine::printProperties() const { 197 vector<H_OpticalElement*>::const_iterator element_i; 198 cout << "Pointeurs des elements du faisceau" << endl; 199 for (element_i = elements.begin(); element_i < elements.end(); element_i++) { 200 cout << (int)(element_i-elements.begin()) << "\t" << (*element_i)->getName() << "\t" << (*element_i)->getS() << endl; 201 } 202 return; 203 } 219 204 220 205 void H_AbstractBeamLine::showElements() const{ … … 254 239 temp = *matrix_i; 255 240 cout << "Matrix for transport until s=" << (*element_i)->getS() + (*element_i)->getLength() << "m (" << (*element_i)->getName() << "). " << endl; 256 printMatrix( temp);241 printMatrix(&temp); 257 242 cout << endl; 258 243 } … … 271 256 // getting rid of drifts before calculating 272 257 for(element_i = elements.begin(); element_i < elements.end(); element_i++) { 273 if((*element_i)->getType() == DRIFT) { delete (*element_i);elements.erase(element_i); }258 if((*element_i)->getType() == DRIFT) {elements.erase(element_i); } 274 259 } 275 260 … … 296 281 temp_elements.push_back(dr); 297 282 } 298 299 // cleaning : avoid some memory leaks300 283 elements.clear(); 301 302 284 for(element_i=temp_elements.begin(); element_i < temp_elements.end(); element_i++) { 303 285 elements.push_back(*element_i); … … 321 303 } 322 304 323 beam_mat.ResizeTo(MDIM,MDIM); 324 beam_mat = calc_mat; 305 *beam_mat = calc_mat; 325 306 return; 326 307 } … … 338 319 } 339 320 340 void H_AbstractBeamLine::draw( const float xmin, const float ymin, const float xmax, const float ymax) const{341 gROOT->SetStyle("Plain");342 TLegend* leg = new TLegend( xmin,ymin,xmax,ymax,"");321 void H_AbstractBeamLine::draw() const{ 322 /* gROOT->SetStyle("Plain"); 323 TLegend* leg = new TLegend(0.85,0.50,1,1,"Legend"); 343 324 leg->SetBorderSize(1); 344 leg->SetFillColor(kWhite);345 325 TBox* b1 = new TBox(); 346 326 TBox* b2 = new TBox(0,0,10,10); … … 365 345 leg->AddEntry(b7,"RCollimator"); 366 346 leg->Draw(); 347 */ 367 348 /* TLine* l1 = new TLine(0.05,0.5,0.95,0.5); 368 349 TLine* l2 = new TLine(0.1,0.1,0.1,0.9); … … 413 394 } 414 395 415 void H_AbstractBeamLine::drawX(const float a_min, const float a_max , const float scale) const{396 void H_AbstractBeamLine::drawX(const float a_min, const float a_max) const{ 416 397 /// @param a_min defines the size of the drawing 417 398 /// @param a_max defines the size of the drawing 418 /// @param scale allows to multiply the drawing, i.e. changing the units419 399 const int N = getNumberOfElements(); 420 400 for(int i=0;i<N;i++) { … … 422 402 float meight = fabs(a_min); 423 403 float size = (height>meight)?meight:height; 424 float middle = getElement(i)->getX()*URAD *scale;404 float middle = getElement(i)->getX()*URAD; 425 405 if(getElement(i)->getType()!=DRIFT) getElement(i)->draw(middle+size/2.,middle-size/2.); 426 406 } … … 440 420 } 441 421 442 void H_AbstractBeamLine::moveElement(const string &name, const float new_s) {422 void H_AbstractBeamLine::moveElement(const string name, const float new_s) { 443 423 /// @param name identifies the element to move 444 424 /// @param new_s is where to put it … … 453 433 } 454 434 455 void H_AbstractBeamLine::alignElement(const string &name, const float disp_x, const float disp_y) {435 void H_AbstractBeamLine::alignElement(const string name, const float disp_x, const float disp_y) { 456 436 /// @param name identifies the element to move 457 437 /// @param disp_x identifies the displacement to add in x [\f$ \mu m \f$] … … 465 445 } 466 446 } 467 cout<<"<H_AbstractBeamLine> WARNING : Element "<<name<<" not found."<<endl; 447 cout<<"Element "<<name<<" not found."<<endl; 448 if(VERBOSE) cout<<"Element "<<name<<" not found."<<endl; 468 449 return; 469 450 } 470 451 471 void H_AbstractBeamLine::tiltElement(const string &name, const float ang_x, const float ang_y) {452 void H_AbstractBeamLine::tiltElement(const string name, const float ang_x, const float ang_y) { 472 453 /// @param name identifies the element to move 473 454 /// @param ang_x identifies the angle to add in x … … 481 462 } 482 463 } 483 cout<<"<H_AbstractBeamLine> WARNING : Element "<<name<<" not found."<<endl;464 if(VERBOSE) cout<<"Element "<<name<<" not found."<<endl; 484 465 return; 485 466 } … … 500 481 } 501 482 483 /* 502 484 TGraph * H_AbstractBeamLine::getBetaX() const{ 503 485 const int N = elements.size(); … … 648 630 return rely; 649 631 } 650 632 */ -
trunk/external/Hector/H_AbstractBeamLine.h
r1360 r1365 2 2 #define _H_AbstractBeamLine_ 3 3 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 5 * * 6 * --<--<-- A fast simulator --<--<-- * 7 * / --<--<-- of particle --<--<-- * 8 * ----HECTOR----< * 9 * \ -->-->-- transport through -->-->-- * 10 * -->-->-- generic beamlines -->-->-- * 11 * * 12 * JINST 2:P09005 (2007) * 13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 14 * http://www.fynu.ucl.ac.be/hector.html * 15 * * 16 * Center for Cosmology, Particle Physics and Phenomenology * 17 * Universite catholique de Louvain * 18 * Louvain-la-Neuve, Belgium * 19 * * 20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 4 /* 5 ---- Hector the simulator ---- 6 A fast simulator of particles through generic beamlines. 7 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 21 8 9 http://www.fynu.ucl.ac.be/hector.html 10 11 Centre de Physique des Particules et de Phénoménologie (CP3) 12 Université Catholique de Louvain (UCL) 13 */ 22 14 23 15 /// \file H_AbstractBeamLine.h 24 16 /// \brief Class aiming at simulating the LHC beamline. 25 17 /// 26 /// Units : angles [ urad], distances [um], energies [GeV], c=[1].18 /// Units : angles [µrad], distances [µm], energies [GeV], c=[1]. 27 19 28 20 /// default length of the beam line … … 37 29 // ROOT #includes 38 30 #include "TMatrix.h" 39 #include "TGraph.h"31 ////#include "TGraph.h" 40 32 41 33 // local #includes … … 47 39 48 40 public: 41 void init(const float ); 49 42 /// Constructors, destructor and operator 50 43 //@{ … … 53 46 H_AbstractBeamLine(const H_AbstractBeamLine &); 54 47 H_AbstractBeamLine& operator=(const H_AbstractBeamLine&); 55 H_AbstractBeamLine* clone() const ;56 48 ~H_AbstractBeamLine(); 57 49 //@} 58 50 /// Adds an element to the beamline 51 //@{ 59 52 void add(H_OpticalElement *); 53 void add(H_OpticalElement &); 54 //@} 60 55 /// Returns the (float) length of the beamline 61 constfloat getLength() const { return beam_length;};56 inline float getLength() const { return beam_length;}; 62 57 /// Returns the (int) number of optics element of the beamline, including drifts 63 const unsigned int getNumberOfElements() const { returnelements.size();};58 inline int getNumberOfElements() const { return (int)elements.size();}; 64 59 /// Returns the transport matrix for the whole beam 65 const TMatrix getBeamMatrix() const;60 const TMatrix * getBeamMatrix() const; 66 61 /// Returns the transport matrix for the whole beam, for given energy loss/mass/charge 67 const TMatrix getBeamMatrix(const float , const float, const float ) ;62 const TMatrix * getBeamMatrix(const float , const float, const float ) ; 68 63 /// Returns the transport matrix for a part of the beam from the IP to a given element 69 const TMatrix getPartialMatrix(const H_OpticalElement *) const;64 const TMatrix * getPartialMatrix(const H_OpticalElement *) const; 70 65 /// Returns the transport matrix for a part of the beam from the IP to the ith element 71 const TMatrix getPartialMatrix(const unsigned int ) const;66 const TMatrix * getPartialMatrix(const unsigned int ) const; 72 67 /// Returns the transport matrix for a part of the beam from the IP to a given element, given energy loss/mass/charge 73 const TMatrix getPartialMatrix(const string&, const float, const float, const float);68 const TMatrix * getPartialMatrix(const string, const float, const float, const float); 74 69 /// Returns the ith element of the beamline 75 70 //@{ 76 71 H_OpticalElement * getElement(const unsigned int ); 77 72 const H_OpticalElement * getElement(const unsigned int ) const; 78 73 //@} 79 74 /// Returns a given element of the beamline, choosen by name 80 75 //@{ 81 H_OpticalElement * getElement(const string &);82 H_OpticalElement * getElement(const string&) const;76 H_OpticalElement * getElement(const string ); 77 const H_OpticalElement * getElement(const string ) const; 83 78 //@} 84 79 /// Print some info … … 96 91 /// Computes global transport matrix 97 92 void calcMatrix(); 98 /// Draws the legend of theelements of the beam99 void draw( const float xmin =0.85, const float ymin=0.5, const float xmax=1, const float ymax=1) const;93 /// Draws the elements of the beam 94 void draw() const; 100 95 /// Draws the elements of the beam in the (x,s) plane 101 void drawX(const float, const float , const float scale=1) const;96 void drawX(const float, const float) const; 102 97 /// Draws the elements of the beam in the (y,s) plane 103 98 void drawY(const float, const float) const; 104 99 /// Moves an element in the list, reorders the lists and recomputes the transport matrix 105 void moveElement(const string &, const float );100 void moveElement(const string, const float ); 106 101 /// Moves the given element tranversely by given amounts. 107 void alignElement(const string &, const float, const float);102 void alignElement(const string, const float, const float); 108 103 /// Tilts the given element tranversely by given angles. 109 void tiltElement(const string &, const float, const float);104 void tiltElement(const string, const float, const float); 110 105 /// Offsets all element in X pos from the start position 111 106 void offsetElements(const float start, const float offset); 112 107 /// Draws the beta functions, from MAD 113 108 //@{ 114 TGraph * getBetaX() const;115 TGraph * getBetaY() const;109 ////TGraph * getBetaX() const; 110 ////TGraph * getBetaY() const; 116 111 //@} 117 112 /// Draws the dispersion functions, from MAD 118 113 //@{ 119 TGraph * getDX() const;120 TGraph * getDY() const;114 ////TGraph * getDX() const; 115 ////TGraph * getDY() const; 121 116 //@} 122 117 /// Draws the relative position functions, from MAD 123 118 //@{ 124 TGraph * getRelX() const;125 TGraph * getRelY() const;119 ////TGraph * getRelX() const; 120 ////TGraph * getRelY() const; 126 121 //@} 127 122 128 123 129 124 private: 125 /// list of all optics elements, including drifts 126 vector<H_OpticalElement*> elements; 130 127 /// list of matrices, 1 matrix = the transport till the end of each element 131 128 vector<TMatrix> matrices; 132 129 /// transport matrix for the whole beam 133 TMatrix beam_mat; 134 /// list of all optics elements, including drifts 135 vector<H_OpticalElement*> elements; 130 TMatrix * beam_mat; 136 131 /// Orderting method for the vector of H_OpticalElement* 137 struct ordering{ bool operator()(const H_OpticalElement* el1, const H_OpticalElement* el2) const {return (*el1 < *el2);}}; 138 // private method for copying the contents of "elements" 139 void cloneElements(const H_AbstractBeamLine&); 132 struct ordering{ bool operator()(H_OpticalElement* el1, H_OpticalElement* el2) const { return (*el1 < *el2);}}; 140 133 141 134 protected: 142 void init(const float );143 135 /// total length of the beamline 144 136 float beam_length; 145 friend std::ostream& operator<< (std::ostream& os, const H_AbstractBeamLine& be);146 137 }; 147 138 -
trunk/external/Hector/H_Aperture.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 18 11 19 12 /// \file H_Aperture.cc … … 27 20 using namespace std; 28 21 29 H_Aperture::H_Aperture() : 30 type_(NONE), x1(0), x2(0), x3(0), x4(0), fx(0), fy(0), aptypestring(getApertureString()) { 22 H_Aperture::H_Aperture() { 23 type = NONE; 24 setApertureString(); 25 x1 = 0; 26 x2 = 0; 27 x3 = 0; 28 x4 = 0; 29 fx = 0; 30 fy = 0; 31 31 } 32 32 33 H_Aperture::H_Aperture(const int dtype, const float size1, const float size2, const float size3, const float size4, const float posx, const float posy) : 34 type_(dtype), x1(size1), x2(size2), x3(size3), x4(size4), fx(posx), fy(posy), aptypestring(getApertureString()) { 33 H_Aperture::H_Aperture(const int dtype, const float size1, const float size2, const float size3, const float size4, const float posx, const float posy) { 35 34 /// @param dtype defines the aperture shape 36 35 /// @param size1, size2, size3, size4 are the geometrical sizes (length/width or great/small radii) in m 37 36 /// @param posx, posy are the (x,y) coordinates of the center of the aperture [m] 37 type = dtype; 38 setApertureString(); 39 x1 = size1; 40 x2 = size2; 41 x3 = size3; 42 x4 = size4; 43 fx = posx; 44 fy = posy; 38 45 } 39 46 40 H_Aperture::H_Aperture(const H_Aperture& ap) : 41 type_(ap.type_), x1(ap.x1), x2(ap.x2), x3(ap.x3), x4(ap.x4), fx(ap.fx), fy(ap.fy), aptypestring(ap.aptypestring) { 47 H_Aperture::H_Aperture(const H_Aperture& ap) { 48 type = ap.type; 49 aptypestring = ap.aptypestring; 50 x1 = ap.x1; 51 x2 = ap.x2; 52 x3 = ap.x3; 53 x4 = ap.x4; 54 fx = ap.fx; 55 fy = ap.fy; 42 56 } 43 57 44 58 H_Aperture& H_Aperture::operator=(const H_Aperture& ap) { 45 59 if(this==&ap) return *this; 46 type_ = ap.type_; 60 type = ap.type; 61 aptypestring = ap.aptypestring; 47 62 x1 = ap.x1; 48 63 x2 = ap.x2; … … 51 66 fx = ap.fx; 52 67 fy = ap.fy; 53 aptypestring = ap.aptypestring;54 68 return *this; 55 69 } 56 70 57 std::ostream& operator<< (std::ostream& os, const H_Aperture& ap) {58 os << "Aperture shape:" << ap.aptypestring << ", parameters "<<ap.x1<<", "<<ap.x2<<", "<<ap.x3<<", "<<ap.x4<<endl;59 os << " \t Center : " << ap.fx << "," << ap.fy << endl;60 return os;61 }62 63 71 void H_Aperture::printProperties() const { 64 cout << *this; 72 cout << "Aperture shape:" << getTypeString() << ", parameters"<<x1<<", "<<x2<<", "<<x3<<", "<<x4<<endl; 73 cout << " \t Center : " << fx << "," << fy << endl; 65 74 return; 66 75 } … … 74 83 } 75 84 76 bool H_Aperture::isInside(const float , const float) const {85 bool H_Aperture::isInside(const float x, const float y) const { 77 86 /// @param x, y are the (x,y) coordinates of the proton, in [m] 78 //cout << "aperture::isInside" << endl;87 cout << "aperture::isInside" << endl; 79 88 return false; 80 89 } 81 90 82 91 83 /*void H_Aperture::setApertureString() {84 switch (type _) {92 void H_Aperture::setApertureString() { 93 switch (type) { 85 94 case NONE: aptypestring = NONENAME; break; 86 95 case RECTANGULAR: aptypestring = RECTANGULARNAME; break; … … 91 100 } 92 101 } 93 */94 95 const string H_Aperture::getApertureString() const {96 string str;97 switch (type_) {98 case NONE: str = NONENAME; break;99 case RECTANGULAR: str = RECTANGULARNAME; break;100 case ELLIPTIC: str = ELLIPTICNAME; break;101 case CIRCULAR: str = CIRCULARNAME; break;102 case RECTELLIPSE: str = RECTELLIPSENAME; break;103 default: str = NONENAME; break;104 }105 return str;106 }107 -
trunk/external/Hector/H_Aperture.h
r1360 r1365 2 2 #define _H_Aperture_ 3 3 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 5 * * 6 * --<--<-- A fast simulator --<--<-- * 7 * / --<--<-- of particle --<--<-- * 8 * ----HECTOR----< * 9 * \ -->-->-- transport through -->-->-- * 10 * -->-->-- generic beamlines -->-->-- * 11 * * 12 * JINST 2:P09005 (2007) * 13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 14 * http://www.fynu.ucl.ac.be/hector.html * 15 * * 16 * Center for Cosmology, Particle Physics and Phenomenology * 17 * Universite catholique de Louvain * 18 * Louvain-la-Neuve, Belgium * 19 * * 20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 4 /* 5 ---- Hector the simulator ---- 6 A fast simulator of particles through generic beamlines. 7 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 8 9 http://www.fynu.ucl.ac.be/hector.html 10 11 Centre de Physique des Particules et de Phénoménologie (CP3) 12 Université Catholique de Louvain (UCL) 13 */ 21 14 22 15 /// \file H_Aperture.h … … 25 18 // C++ #defines 26 19 #include <string> 27 #include <iostream>28 20 using namespace std; 29 21 30 22 // local #defines 31 enum {NONE=0, RECTANGULAR, ELLIPTIC, CIRCULAR, RECTELLIPSE}; 23 #define NONE 0 24 #define RECTANGULAR 1 25 #define ELLIPTIC 2 26 #define CIRCULAR 3 27 #define RECTELLIPSE 4 32 28 #define NONENAME "None " 33 29 #define RECTANGULARNAME "Rectangle " … … 46 42 H_Aperture(const H_Aperture&); 47 43 H_Aperture& operator=(const H_Aperture&); 48 virtual ~H_Aperture() { }; 49 virtual H_Aperture* clone() const { return new H_Aperture(type_,x1,x2,x3,x4,fx,fy); }; 44 virtual ~H_Aperture() {return;}; 50 45 //@} 51 46 52 47 /// Prints the aperture features 53 //virtual void printProperties() const; 54 void printProperties() const; 48 virtual void printProperties() const; 55 49 /// Draws the aperture shape 56 virtual void draw( const float) const {return;};50 virtual void draw() const {return;}; 57 51 /// Sets the (x,y) position in [m] 58 52 void setPosition(const float,const float); … … 60 54 virtual bool isInside(const float, const float) const; 61 55 /// Returns the (int) type of aperture 62 inline int getType() const {return type _;};56 inline int getType() const {return type;}; 63 57 /// Returns the (string) type of the aperture 64 58 inline const string getTypeString() const { return aptypestring; } … … 67 61 protected: 68 62 /// Aperture shape (either RECTANGULAR or ELLIPTIC or ...) 69 int type_; 63 int type; 64 /// Aperture shape string 65 string aptypestring; 66 70 67 /// Aperture geometrical sizes (length/width or great/small radii) [m] 71 68 //@{ 72 69 float x1, x2, x3, x4; 73 70 //@} 71 74 72 /// Horizontal coordinate of the aperture center [m] (from the nominal beam position). 75 73 //@{ 76 74 float fx, fy; 77 75 //@} 78 /// Aperture shape string 79 string aptypestring; 80 // Sets the name of the aperture from its type. 81 //void setApertureString(); 82 /// Gets the name of the aperture from its type. 83 const string getApertureString() const; 84 85 friend std::ostream& operator<< (std::ostream& os, const H_Aperture& ap); 76 /// Sets the name of the aperture from its type. 77 void setApertureString(); 86 78 }; 87 79 88 89 80 #endif -
trunk/external/Hector/H_Beam.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 18 11 19 12 /// \file H_Beam.cc … … 22 15 23 16 // ROOT #includes 24 #include "TGraph.h"17 //#include "TGraph.h" 25 18 #include "TRandom.h" 26 19 … … 62 55 }; 63 56 64 void H_Beam::createBeamParticles(const unsigned int Number_of_particles, const double p_mass, const double p_charge, TRandom* r) { 57 void H_Beam::createBeamParticles(const unsigned int Number_of_particles) { 58 createBeamParticles(Number_of_particles,MP,QP); 59 } 60 61 void H_Beam::createBeamParticles(const unsigned int Number_of_particles, const double p_mass, const double p_charge) { 65 62 beamParticles.clear(); 66 63 Nparticles = (Number_of_particles<1) ? 1 : Number_of_particles; … … 69 66 p.setPosition(fx_ini,fy_ini,tx_ini,ty_ini,fs_ini); 70 67 p.setE(fe_ini); 71 p.smearPos(x_disp,y_disp ,r);72 p.smearAng(tx_disp,ty_disp ,r);73 p.smearE(e_disp ,r);74 p.smearS(s_disp ,r);68 p.smearPos(x_disp,y_disp); 69 p.smearAng(tx_disp,ty_disp); 70 p.smearE(e_disp); 71 p.smearS(s_disp); 75 72 if (VERBOSE) {if (i==0) cout << " x_ini , tx_ini " << p.getX() << " " << p.getTX() << endl;} 76 73 beamParticles.push_back(p); 77 74 } 78 75 } 79 80 //void H_Beam::particleGun(const unsigned int Number_of_particles, const float E_min=BE, const float E_max=BE, const float fs_min=0, const float fs_max=0, const float fx_min=0, const float fx_max=0, const float fy_min=0, const float fy_max=0, const float tx_min=-PI/2., const float tx_max=PI/2., const float ty_min=-PI/2., const float ty_max=PI/2., const float p_mass=MP, const double p_charge=QP) {81 void H_Beam::particleGun(const unsigned int Number_of_particles, const float E_min, const float E_max, const float fs_min, const float fs_max, const float fx_min, const float fx_max, const float fy_min, const float fy_max, const float tx_min, const float tx_max, const float ty_min, const float ty_max, const float p_mass, const double p_charge, const bool flat, TRandom* r) {82 beamParticles.clear();83 Nparticles = (Number_of_particles<2) ? 2 : Number_of_particles;84 float gx,gy,gs,gtx,gty,gE;85 for (unsigned int i=0; i<Nparticles; i++) {86 H_BeamParticle p(p_mass,p_charge);87 if (flat) {88 gx = r->Uniform(fx_min,fx_max);89 gy = r->Uniform(fy_min,fy_max);90 gs = r->Uniform(fs_min,fs_max);91 gtx = r->Uniform(tx_min,tx_max);92 gty = r->Uniform(ty_min,ty_max);93 gE = r->Uniform(E_min,E_max);94 } else {95 gx = r->Gaus((fx_min+fx_max)/2,(-fx_min+fx_max)/2);96 gy = r->Gaus((fy_min+fy_max)/2,(-fy_min+fy_max)/2);97 gs = r->Gaus((fs_min+fs_max)/2,(-fs_min+fs_max)/2);98 gtx = r->Gaus((tx_min+tx_max)/2,(-tx_min+tx_max)/2);99 gty = r->Gaus((ty_min+ty_max)/2,(-ty_min+ty_max)/2);100 gE = r->Gaus ((E_min+E_max)/2,(-E_min+E_max)/2);101 }102 p.setPosition(gx,gy,gtx,gty,gs);103 p.setE(gE);104 beamParticles.push_back(p);105 }106 return;107 }108 109 76 110 77 void H_Beam::createXScanningBeamParticles(const unsigned int Number_of_particles, const float fx_max) { … … 182 149 } 183 150 151 void H_Beam::computePath(const H_AbstractBeamLine * beamline) { 152 computePath(beamline,false); 153 } 154 184 155 /// Propagates the beam until a given s 185 156 void H_Beam::propagate(const float position) { … … 188 159 particle_i->propagate(position); 189 160 } 161 } 162 163 void H_Beam::emitGamma(const double gee, const double gq2) { 164 /// @param gee = \f$ E_{\gamma} \f$ is the photon energy 165 /// @param gq2 = \f$ Q^2 < 0 \f$ is virtuality of photon \f$ Q^{2} = E^{2}-\vec{k}^{2} \f$ 166 emitGamma(gee,gq2,0,2*PI); 190 167 } 191 168 … … 237 214 return EY2; 238 215 } 239 216 /* 240 217 TGraphErrors * H_Beam::getBetaX(const float length, const unsigned int number_of_points) { 241 218 /// @param length [m] … … 277 254 return betay; 278 255 } 256 */ 279 257 280 258 float H_Beam::getX(const float s, float& error_on_posx) { … … 313 291 } 314 292 315 vector<TVectorD> H_Beam::getStoppingElements(const H_AbstractBeamLine * beamline, vector<H_OpticalElement>& list, vector<int>& numb) { 316 317 vector<TVectorD> stop_positions; 293 void H_Beam::getStoppingElements(const H_AbstractBeamLine * beamline, vector<H_OpticalElement>& list, vector<int>& numb) { 318 294 vector<H_BeamParticle>::iterator particle_i; 319 295 vector<H_OpticalElement>::iterator element_i; … … 331 307 if(particle_i->stopped(beamline)) { 332 308 temp_el = *(particle_i->getStoppingElement()); 333 stop_positions.push_back(*(particle_i->getStopPosition()));334 309 if(list.size()==0) { 335 310 number=1; 336 311 list.push_back(temp_el); 337 312 numb.push_back(number); 338 339 313 } else { 340 314 for (element_i = list.begin(), n_i = numb.begin(); element_i < list.end(); element_i++, n_i++) { … … 356 330 } // if particle_i->stopped 357 331 }// for particle_i 358 return stop_positions;359 332 } // H_Beam::getStoppingElements 360 333 … … 379 352 } 380 353 381 std::ostream& operator<< (std::ostream& os, const H_Beam& be){354 void H_Beam::printProperties() const { 382 355 vector<H_BeamParticle>::const_iterator particle_i; 383 cout << "There are " << be.Nparticles << " in the beam." << endl; 384 for (particle_i = be.beamParticles.begin(); particle_i < be.beamParticles.end(); particle_i++) { 385 cout << *particle_i; 386 } 387 return os; 356 cout << "There are " << Nparticles << " in the beam." << endl; 357 for (particle_i = beamParticles.begin();particle_i < beamParticles.end(); particle_i++) { 358 particle_i->printProperties(); 359 } 388 360 } 389 361 … … 400 372 } 401 373 } // H_Beam::printStoppingElements 402 403 TH2F * H_Beam::drawAngleProfile(const float s) { 404 /// not a const method because does a propagate to s! 405 char title[50]; 406 sprintf(title,"Beam profile at %.2f m",s); 407 vector<H_BeamParticle>::iterator particle_i; 408 float xmax, xmin, ymax, ymin; 409 float xx, yy, xborder, yborder; 410 411 particle_i=beamParticles.begin(); 412 xmin = particle_i->getTX(); 413 xmax = particle_i->getTX(); 414 ymin = particle_i->getTY(); 415 ymax = particle_i->getTY(); 416 417 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) { 418 particle_i->propagate(s); 419 xx = particle_i->getTX(); 420 yy = particle_i->getTY(); 421 422 xmax = xx>xmax ? xx : xmax; 423 ymax = yy>ymax ? yy : ymax; 424 xmin = xx<xmin ? xx : xmin; 425 ymin = yy<ymin ? yy : ymin; 426 } 427 428 // in order to avoid some drawing problems, when the beam divergence is null 429 if(!(xmax || xmin)) xmax +=0.1; 430 if(!(ymax || ymin)) xmax +=0.1; 431 432 if(xmax == xmin) xmax *= 1.1; 433 if(ymax == ymin) ymax *= 1.1; 434 435 xborder = (xmax-xmin)*0.2; 436 yborder = (ymax-ymin)*0.2; 437 438 xmax += xborder; 439 xmin -= xborder; 440 ymax += yborder; 441 ymin -= yborder; 442 443 TH2F * profile = new TH2F("profile",title,10000,xmin,xmax,1000,ymin,ymax); 444 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) { 445 profile->Fill(particle_i->getTX(), particle_i->getTY()); 446 } 447 return profile; 448 } 449 450 374 /* 451 375 TH2F * H_Beam::drawProfile(const float s) { 452 376 /// not a const method because does a propagate to s! … … 475 399 476 400 // in order to avoid some drawing problems, when the beam divergence is null 477 if(!(xmax || xmin)) xmax +=0.1;478 if(!(ymax || ymin)) xmax +=0.1;479 480 401 if(xmax == xmin) xmax += 0.1; 481 402 if(ymax == ymin) ymax += 0.1; … … 494 415 } 495 416 return profile; 496 } 497 417 }*/ 418 /* 498 419 TMultiGraph * H_Beam::drawBeamX(const int color) const { 499 420 int mycolor = color; … … 519 440 return beam_profile_y; 520 441 } 442 */ -
trunk/external/Hector/H_Beam.h
r1360 r1365 2 2 #define _H_Beam_ 3 3 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 5 * * 6 * --<--<-- A fast simulator --<--<-- * 7 * / --<--<-- of particle --<--<-- * 8 * ----HECTOR----< * 9 * \ -->-->-- transport through -->-->-- * 10 * -->-->-- generic beamlines -->-->-- * 11 * * 12 * JINST 2:P09005 (2007) * 13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 14 * http://www.fynu.ucl.ac.be/hector.html * 15 * * 16 * Center for Cosmology, Particle Physics and Phenomenology * 17 * Universite catholique de Louvain * 18 * Louvain-la-Neuve, Belgium * 19 * * 20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 4 /* 5 ---- Hector the simulator ---- 6 A fast simulator of particles through generic beamlines. 7 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 8 9 http://www.fynu.ucl.ac.be/hector.html 10 11 Centre de Physique des Particules et de Phénoménologie (CP3) 12 Université Catholique de Louvain (UCL) 13 */ 21 14 22 15 /// \file H_Beam.h … … 31 24 32 25 // ROOT #includes 33 #include "TH2F.h" 34 #include "TGraphErrors.h" 35 #include "TMultiGraph.h" 36 #include "TMath.h" 26 ////#include "TH2F.h" 27 ////#include "TGraphErrors.h" 28 ////#include "TMultiGraph.h" 37 29 38 30 // local #includes … … 63 55 ~H_Beam(); 64 56 //@} 65 /// Fills the beam with particles in given position/angle/energy intervals (flat distribution)66 void particleGun(const unsigned int Number_of_particles, const float E_min=BE, const float E_max=BE, const float fs_min=0, const float fs_max=0, const float fx_min=0, const float fx_max=0, const float fy_min=0, const float fy_max=0, const float tx_min=-TMath::Pi()/2., const float tx_max=TMath::Pi()/2., const float ty_min=-TMath::Pi()/2., const float ty_max=TMath::Pi()/2., const float p_mass=MP, const double p_charge=QP, const bool flat = true, TRandom* r=gRandom);67 68 57 /// Fills the beam with particles 69 void createBeamParticles(const unsigned int Number_of_particles, const double p_mass=MP, const double p_charge=QP, TRandom* r=gRandom); 58 //@{ 59 void createBeamParticles(const unsigned int , const double , const double ); 60 void createBeamParticles(const unsigned int); 61 //@} 70 62 /// Fills the beam with particles with incremental offset and no initial transverse momentum 71 63 //@{ … … 86 78 void add(const H_BeamParticle&); 87 79 /// Compute the position of each particle in the beamline 88 void computePath(const H_AbstractBeamLine * beamline, const bool NonLinear=false); 80 //@{ 81 void computePath(const H_AbstractBeamLine *, const bool); 82 void computePath(const H_AbstractBeamLine *); 83 //@} 89 84 // Photon emission by the particle 90 void emitGamma(const double gee, const double gq2, const double phimin=0, const double phimax=2*TMath::Pi()); 85 // @{ 86 void emitGamma(const double, const double, const double, const double); 87 void emitGamma(const double, const double); 88 //@} 91 89 // Propagates the beam until a given s 92 90 void propagate(const float ); … … 98 96 /// Draws the \f$ \beta \f$ function of the beam 99 97 //@{ 100 101 TGraphErrors * getBetaY(const float, const unsigned int);98 //// TGraphErrors * getBetaX(const float, const unsigned int); 99 //// TGraphErrors * getBetaY(const float, const unsigned int); 102 100 //@} 103 101 /// Returns the position of the beam … … 108 106 /// Returns the emittance \f$ \epsilon \f$ of the beam in x and y 109 107 //@{ 110 inline const float getEmittanceX() const {111 108 inline const float getEmittanceX() const { 109 if(!x_disp*tx_disp) cout<<"Warning : Degenerate Beam : x-emittance = 0"<<endl; 112 110 return x_disp * tan(tx_disp/URAD)/URAD; 113 111 } … … 140 138 unsigned int getStoppedNumber(const H_AbstractBeamLine *); 141 139 /// Returns the list of the stopping elements in the beamline 142 v ector<TVectorD>getStoppingElements(const H_AbstractBeamLine *, vector<H_OpticalElement>&, vector<int>&);140 void getStoppingElements(const H_AbstractBeamLine *, vector<H_OpticalElement>&, vector<int>&); 143 141 /// Prints the initial parameters 144 142 void printInitialState() const; 145 143 /// Prints the properties for each particle 146 void printProperties() const {cout << *this; return;}144 void printProperties() const; 147 145 /// Prints the list of the stopping elements in the beamline 148 146 void printStoppingElements(const vector<H_OpticalElement>&, const vector<int>&) const; … … 150 148 const int getNumberOfBeamParticles() const {return Nparticles;} 151 149 /// Draws the beam profile at a given s 152 TH2F * drawProfile(const float); 153 TH2F * drawAngleProfile(const float); 150 //// TH2F * drawProfile(const float); 154 151 /// Draws the beam width and height 155 152 //@{ 156 TMultiGraph * drawBeamX(const int) const;157 TMultiGraph * drawBeamY(const int) const ;153 //// TMultiGraph * drawBeamX(const int) const; 154 //// TMultiGraph * drawBeamY(const int) const ; 158 155 //@} 159 156 … … 183 180 /// Number of particles in this beam 184 181 unsigned int Nparticles; 185 186 friend std::ostream& operator<< (std::ostream& os, const H_Beam& be);187 182 }; 188 183 -
trunk/external/Hector/H_BeamLine.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 18 11 19 12 /// \file H_BeamLine.cc … … 24 17 #include <fstream> 25 18 #include <sstream> 26 #include <cstdlib>27 19 28 20 // local #includes … … 47 39 // are opposite to Wille's. See fill() for more details. 48 40 49 H_BeamLine::H_BeamLine(): H_AbstractBeamLine(), 50 direction(1), ips(0), ipx(0), ipy(0), iptx(0), ipty(0) { 51 } 52 53 H_BeamLine::H_BeamLine(const int si, const float length) : H_AbstractBeamLine(length), 54 direction((si >= abs(si)) ? 1 : -1), ips(0), ipx(0), ipy(0), iptx(0), ipty(0) { 55 } 56 57 H_BeamLine::H_BeamLine(const H_BeamLine& beam) : H_AbstractBeamLine(beam), 58 direction(beam.direction), ips(beam.ips), ipx(beam.ipx), ipy(beam.ipy), iptx(beam.iptx), ipty(beam.ipty) { 59 } 60 61 H_BeamLine& H_BeamLine::operator=(const H_BeamLine& beam) { 62 if(this==&beam) return *this; 63 H_AbstractBeamLine::operator=(beam); // call the mother's operator= 41 H_BeamLine::H_BeamLine(const int si, const float length) : H_AbstractBeamLine(length){ 42 direction = (si >= abs(si)) ? 1 : -1; 43 ips=0; 44 ipx=0; 45 ipy=0; 46 iptx=0; 47 ipty=0; 48 } 49 50 H_BeamLine::H_BeamLine(const H_BeamLine& beam) : H_AbstractBeamLine(beam) { 64 51 direction = beam.direction; 65 52 ips = beam.ips; … … 68 55 iptx = beam.iptx; 69 56 ipty = beam.ipty; 57 } 58 59 H_BeamLine& H_BeamLine::operator=(const H_BeamLine& beam) { 60 if(this==&beam) return *this; 61 direction = beam.direction; 62 ips = beam.ips; 63 ipx = beam.ipx; 64 ipy = beam.ipy; 65 iptx = beam.iptx; 66 ipty = beam.ipty; 70 67 return *this; 71 68 } 72 69 73 void H_BeamLine::findIP(const string& filename, const string& ipname) { 70 void H_BeamLine::findIP(const string filename) { 71 findIP(filename,"IP5"); 72 return; 73 } 74 75 void H_BeamLine::findIP(const string filename, const string ipname) { 74 76 // searches for the IP position in the extended table. 75 77 ifstream tabfile(filename.c_str()); 76 if (! tabfile.is_open()) cout << " <H_BeamLine>ERROR: I Can't open \"" << filename << "\"" << endl;78 if (! tabfile.is_open()) cout << "\t ERROR: I Can't open \"" << filename << "\"" << endl; 77 79 bool found = false; 78 80 int N_col=0; … … 112 114 } // while (!eof) 113 115 114 if (!found) cout << " <H_BeamLine>ERROR ! IP not found." << endl;116 if (!found) cout << "\t ERROR ! IP not found." << endl; 115 117 tabfile.close(); 116 118 } … … 125 127 } 126 128 127 128 void H_BeamLine::fill(const string& filename, const int dir, const string& ipname) { 129 void H_BeamLine::fill(const string filename) { 130 fill(filename,1,"IP5"); 131 return; 132 } 133 134 135 void H_BeamLine::fill(const string filename, const int dir, const string ipname) { 129 136 string headers[40]; // table of all column headers 130 137 int header_type[40]; // table of all column types 131 138 findIP(filename,ipname); 132 139 ifstream tabfile(filename.c_str()); 133 if (! tabfile.is_open()) cout << " <H_BeamLine>ERROR: I Can't open \"" << filename << "\"" << endl;140 if (! tabfile.is_open()) cout << "\t ERROR: I Can't open \"" << filename << "\"" << endl; 134 141 if(VERBOSE) { 135 142 cout<<"Using file : "<< filename <<" in the "<<((direction>0)?"positive":"negative")<<" direction."<<endl; … … 187 194 } 188 195 196 ap = 0; //init 197 if(e.aper_1 !=0 || e.aper_2 !=0 || e.aper_3 !=0 || e.aper_4 != 0) { 198 e.aper_1 *= URAD; 199 e.aper_2 *= URAD; 200 e.aper_3 *= URAD; 201 e.aper_4 *= URAD; // in [m] in the tables ! 202 203 if(strstr(e.apertype.c_str(),"RECTELLIPSE")) ap = new H_RectEllipticAperture(e.aper_1,e.aper_2,e.aper_3,e.aper_4,0,0); 204 else if(strstr(e.apertype.c_str(),"CIRCLE")) ap = new H_CircularAperture(e.aper_1,0,0); 205 else if(strstr(e.apertype.c_str(),"RECTANGLE")) ap = new H_RectangularAperture(e.aper_1,e.aper_2,0,0); 206 else if(strstr(e.apertype.c_str(),"ELLIPSE")) ap = new H_EllipticAperture(e.aper_1,e.aper_2,0,0); 207 } 208 189 209 if(el!=0) { 190 ap = 0; //init 191 if(e.aper_1 !=0 || e.aper_2 !=0 || e.aper_3 !=0 || e.aper_4 != 0) { 192 e.aper_1 *= URAD; 193 e.aper_2 *= URAD; 194 e.aper_3 *= URAD; 195 e.aper_4 *= URAD; // in [m] in the tables ! 196 197 if(strstr(e.apertype.c_str(),"RECTELLIPSE")) 198 ap = new H_RectEllipticAperture(e.aper_1,e.aper_2,e.aper_3,e.aper_4,0,0); 199 else if(strstr(e.apertype.c_str(),"CIRCLE")) 200 ap = new H_CircularAperture(e.aper_1,0,0); 201 else if(strstr(e.apertype.c_str(),"RECTANGLE")) 202 ap = new H_RectangularAperture(e.aper_1,e.aper_2,0,0); 203 else if(strstr(e.apertype.c_str(),"ELLIPSE")) 204 ap = new H_EllipticAperture(e.aper_1,e.aper_2,0,0); 210 if (direction<0) { 211 el->setBetaX(e.betx); 212 el->setBetaY(e.bety); 213 el->setDX(e.dx); 214 el->setDY(e.dy); 215 el->setRelX(e.x); 216 el->setRelY(e.y); 217 } else { 218 el->setBetaX(previous_betax); 219 el->setBetaY(previous_betay); 220 el->setDX(previous_dx); 221 el->setDY(previous_dy); 222 el->setRelX(previous_x); 223 el->setRelY(previous_y); 224 } 225 if(ap!=0) { 226 el->setAperture(ap); 227 // delete ap; // ap deleted in H_AbstractBeamLine::~H_AbstractBeamLine 205 228 } 206 229 207 if (direction<0) { 208 el->setBetaX(e.betx); el->setBetaY(e.bety); 209 el->setDX(e.dx); el->setDY(e.dy); 210 el->setRelX(e.x); el->setRelY(e.y); 211 } else { 212 el->setBetaX(previous_betax); el->setBetaY(previous_betay); 213 el->setDX(previous_dx); el->setDY(previous_dy); 214 el->setRelX(previous_x); el->setRelY(previous_y); 215 } 216 if(ap!=0) { 217 el->setAperture(ap); 218 delete ap; 219 // if memory is allocated, i.e. if ap!=0, 220 // it wont be deallocated in H_OpticalElement::~H_OpticalElement 221 // ap should then be deleted here 222 } // memory leak-free! 223 224 /// Parses all the elements from the beginning of the file, 225 /// but only keeps the ones from the IP till the desired length 226 if(e.s>=0 && e.s<beam_length) add(el); 227 else { delete el;} 228 // NB : if "el" is added to the beamline, it will be 229 // deleted by H_AbstractBeamLine::~H_AbstractBeamLine 230 // Otherwise, it should be deleted explicitly 231 230 /// Parses all the elements, but only keeps the ones from the IP till the desired length 231 if(e.s>=0 && e.s<beam_length) add(*el); 232 233 // delete el; // el deleted in H_AbstractBeamLine::~H_AbstractBeamLine 232 234 } // if(el!=0) 233 235 } // if (found) -
trunk/external/Hector/H_BeamLine.h
r1360 r1365 2 2 #define _H_BeamLine_ 3 3 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 5 * * 6 * --<--<-- A fast simulator --<--<-- * 7 * / --<--<-- of particle --<--<-- * 8 * ----HECTOR----< * 9 * \ -->-->-- transport through -->-->-- * 10 * -->-->-- generic beamlines -->-->-- * 11 * * 12 * JINST 2:P09005 (2007) * 13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 14 * http://www.fynu.ucl.ac.be/hector.html * 15 * * 16 * Center for Cosmology, Particle Physics and Phenomenology * 17 * Universite catholique de Louvain * 18 * Louvain-la-Neuve, Belgium * 19 * * 20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 4 /* 5 ---- Hector the simulator ---- 6 A fast simulator of particles through generic beamlines. 7 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 8 9 http://www.fynu.ucl.ac.be/hector.html 10 11 Centre de Physique des Particules et de Phénoménologie (CP3) 12 Université Catholique de Louvain (UCL) 13 */ 21 14 22 15 /// \file H_BeamLine.h … … 39 32 /// Constructors and destructor 40 33 //@{ 41 H_BeamLine() ;34 H_BeamLine():H_AbstractBeamLine() {direction=1; ips=0; }; 42 35 H_BeamLine(const H_BeamLine& ); 43 36 H_BeamLine(const int, const float); 44 37 H_BeamLine& operator=(const H_BeamLine& ); 45 ~H_BeamLine() { };38 ~H_BeamLine() {return;}; 46 39 //@ 47 40 /// Finds the IP position (s) from the MAD table. Should be "IP5" or "IP1". 48 void findIP(const string& filename, const string& ipname="IP5"); 41 //@{ 42 void findIP(const string); 43 void findIP(const string, const string); 44 //@} 49 45 /// Reader for the external MAD table 50 void fill(const string& filename, const int dir=1, const string& ipname="IP5"); 46 //@{ 47 void fill(const string); 48 void fill(const string, const int, const string ); 49 //@} 51 50 /// Returns the IP position (s) 52 double getIP() const{return ips;};51 double getIP() {return ips;}; 53 52 /// Returns positions and angles of beam at IP 54 53 double* getIPProperties(); -
trunk/external/Hector/H_BeamLineParser.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 18 11 19 12 /// \file H_BeamLineParser.cc 20 13 /// \brief Reader for madx tables 21 14 /// 22 /// Notes : V erifier que tous les SBEND sont toujours appeles MB. Et seulement comme ca!23 /// V erifier qu'il n'y a pas de problemesd'inversion H/V QUADRUPOLES15 /// Notes : Vï¿œifier que tous les SBEND sont toujours appelï¿œ MB. Et seulement comme ï¿œ ! 16 /// Vï¿œifier qu'il n'y a pas de problï¿œe d'inversion H/V QUADRUPOLES 24 17 /// no distinction between H and Vkickers ? 25 18 /// The identification of the element is based on the values of k1l, k2l, hkick, vkick and on their name. … … 38 31 39 32 /// Identifies the column content from its header 40 int column_identification(const string &header) {33 int column_identification(const string header) { 41 34 // identifies the column type from its name 42 35 -
trunk/external/Hector/H_BeamLineParser.h
r1360 r1365 2 2 #define _H_BeamLineParser_ 3 3 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 5 * * 6 * --<--<-- A fast simulator --<--<-- * 7 * / --<--<-- of particle --<--<-- * 8 * ----HECTOR----< * 9 * \ -->-->-- transport through -->-->-- * 10 * -->-->-- generic beamlines -->-->-- * 11 * * 12 * JINST 2:P09005 (2007) * 13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 14 * http://www.fynu.ucl.ac.be/hector.html * 15 * * 16 * Center for Cosmology, Particle Physics and Phenomenology * 17 * Universite catholique de Louvain * 18 * Louvain-la-Neuve, Belgium * 19 * * 20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 4 /* 5 ---- Hector the simulator ---- 6 A fast simulator of particles through generic beamlines. 7 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 8 9 http://www.fynu.ucl.ac.be/hector.html 10 11 Centre de Physique des Particules et de Phénoménologie (CP3) 12 Université Catholique de Louvain (UCL) 13 */ 21 14 22 15 /// \file H_BeamLineParser.h … … 78 71 */ 79 72 80 extern int column_identification(const string &);73 extern int column_identification(const string ); 81 74 82 75 /// \brief Reader for madx tables to use in H_BeamLine … … 89 82 //@{ 90 83 H_BeamLineParser() {init();} 91 ~H_BeamLineParser() { }84 ~H_BeamLineParser() {return;} 92 85 //@} 93 86 void init(); -
trunk/external/Hector/H_BeamParticle.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 18 11 19 12 /// \file H_BeamParticle.cc … … 26 19 // c++ #includes 27 20 #include <iostream> 28 #include <fstream>29 21 #include <iomanip> 30 22 #include <vector> … … 35 27 #include "H_Parameters.h" 36 28 #include "TRandom.h" 37 #include "TVectorD.h" 29 //#include "TView.h" 30 //#include "TPolyLine3D.h" 38 31 #ifdef _include_pythia_ 39 32 #include "TPythia6.h" 33 #include "TRandom.h" 40 34 #endif 41 42 35 // local #includes 43 36 #include "H_OpticalElement.h" … … 47 40 using namespace std; 48 41 49 void H_BeamParticle::init() { // for more efficiency, put the objects away from init!42 void H_BeamParticle::init() { 50 43 mp = MP; 51 44 qp = QP; … … 65 58 } 66 59 67 H_BeamParticle::H_BeamParticle() { 60 H_BeamParticle::H_BeamParticle() { 68 61 init(); 69 62 } 70 63 71 H_BeamParticle::H_BeamParticle(const H_BeamParticle& p): 72 mp(p.mp), qp(p.qp), fs(p.fs), fx(p.fx), fy(p.fy), thx(p.thx), thy(p.thy), 73 energy(p.energy), hasstopped(p.hasstopped), hasemitted(p.hasemitted), 74 isphysical(p.isphysical), stop_position(new TVectorD(*(p.stop_position))), 75 stop_element(0), positions(p.positions) { 64 H_BeamParticle::H_BeamParticle(const H_BeamParticle& p) { 65 mp = p.mp; 66 qp = p.qp; 67 fx = p.fx; 68 fy = p.fy; 69 thx = p.thx; 70 thy = p.thy; 71 fs = p.fs; 72 energy = p.energy; 73 hasstopped = p.hasstopped; 74 hasemitted = p.hasemitted; 75 isphysical = p.isphysical; 76 stop_position = new TVectorD(*(p.stop_position)); 76 77 if(p.hasstopped) stop_element = new H_OpticalElement(*(p.stop_element)); 78 positions = p.positions; 77 79 } 78 80 … … 99 101 stop_position = new TVectorD(*(p.stop_position)); 100 102 if(p.hasstopped) stop_element = new H_OpticalElement(*(p.stop_element)); 101 else stop_element = 0;102 103 positions = p.positions; 103 104 return *this; 104 105 } 105 106 106 constbool H_BeamParticle::stopped(const H_AbstractBeamLine * beamline) {107 bool H_BeamParticle::stopped(const H_AbstractBeamLine * beamline) { 107 108 vector<TVectorD>::const_iterator position_i; 108 109 for(position_i = positions.begin(); position_i < positions.end()-1; position_i++) { … … 110 111 if(beamline->getElement(pos)->getAperture()->getType()!=NONE) { 111 112 bool has_passed_entrance = beamline->getElement(pos)->isInside((*position_i)[INDEX_X],(*position_i)[INDEX_Y]); 112 bool has_passed_exit = beamline->getElement(pos)->isInside((*(position_i+1))[INDEX_X],(*(position_i+1))[INDEX_Y]); 113 // here we should distinguish between particles passing the input or not. 114 // - particles not passing the input are logically stopped at the input position. period. 115 // - particles passing the input but not the output are stopped somewhere in the element 116 // - particles passing both the input and the output could have been stopped somewhere inside too (less likely) 117 if(!has_passed_entrance) { 118 if(VERBOSE) cout<<"Stopped at the entrance of "<<beamline->getElement(pos)->getName(); 119 hasstopped=true; 120 stop_element = const_cast<H_OpticalElement*>(beamline->getElement(pos)); 121 *stop_position = *position_i; 122 if(VERBOSE) cout<<" at s = "<<(*stop_position)[4]<<endl; 123 return hasstopped; 124 } else if (!has_passed_exit) { 125 if(VERBOSE) cout<<"Stopped inside "<<beamline->getElement(pos)->getName(); 126 hasstopped=true; 127 stop_element = const_cast<H_OpticalElement*>(beamline->getElement(pos)); 128 // this should be computed using the element-based method "H_OpticalElement::getHitPosition" (caution : always nonlinear). 129 *stop_position = beamline->getElement(pos)->getHitPosition(*position_i,BE-energy,mp,qp); 130 if(VERBOSE) cout<<" at s = "<<(*stop_position)[4]<<endl; 131 return hasstopped; 132 } 133 /* 113 bool has_passed_exit = beamline->getElement(pos)->isInside((*(position_i+1))[INDEX_X],(*(position_i+1))[INDEX_Y]); 134 114 if(!(has_passed_entrance && has_passed_exit)) { 115 // cout << "p =" << (*position_i)[INDEX_X] << "; el=" << beamline->getElement(pos)->getX()<<endl; 116 // beamline->getElement(position_i-positions.begin())->printProperties(); 135 117 if(VERBOSE) cout<<"particle stopped at "<<(beamline->getElement(pos))->getName(); 136 118 if(VERBOSE) cout<<" (s = "<<(*position_i)[4] << ")" << endl; 119 if(VERBOSE && !has_passed_exit) cout << "Particle stopped inside the element" << endl; 137 120 hasstopped=true; 138 121 stop_element = const_cast<H_OpticalElement*>(beamline->getElement(pos)); … … 140 123 return hasstopped; 141 124 } // if 142 */ 125 // else cout << "outside aperture " << endl; 143 126 } // if 144 127 } // for … … 159 142 } 160 143 161 162 163 144 void H_BeamParticle::smearPos(const double dx,const double dy, TRandom* r) { 164 145 // the beam is centered on (fx,fy) at IP … … 191 172 } 192 173 193 194 174 void H_BeamParticle::set4Momentum(const double px, const double py, const double pz, const double ene) { 195 175 /// @param px, py, pz, ene is \f$ (\vec p , E) [GeV]\f$ … … 207 187 addPosition(fx,thx,fy,thy,fs); 208 188 return; 209 }210 211 void H_BeamParticle::set4Momentum(const TLorentzVector& pmu) {212 /// @param pmu is the particle 4-momentum \f$ p^\mu \f$213 ///214 /// Clears the H_BeamParticle::positions vector.215 set4Momentum(pmu.Px(), pmu.Py(), pmu.Pz(), pmu.E());216 189 } 217 190 … … 242 215 } 243 216 217 void H_BeamParticle::emitGamma(const double gee, const double gq2) { 218 emitGamma(gee,gq2,0,2*PI); 219 return; 220 } 221 244 222 void H_BeamParticle::emitGamma(const double gee, const double gq2, const double phimin, const double phimax) { 245 223 /// @param gee = \f$ E_{\gamma} \f$ is the photon energy … … 271 249 272 250 if( (gq2>q2max) || (gq2<q2min)) { 273 if(VERBOSE) cout<<" <H_BeamParticle> WARNING :Non physical particle ! Q2 (" << q2 << ") and E ("<<gee << ") are not compatible." << endl;251 if(VERBOSE) cout<<"Non physical particle ! Q2 (" << q2 << ") and E ("<<gee << ") are not compatible." << endl; 274 252 isphysical = false; 275 253 } 276 254 277 if(hasemitted) { cout<<" <H_BeamParticle> WARNING :particle has already emitted at least one gamma !"<<endl;}255 if(hasemitted) { cout<<"particle has already emitted at least one gamma !"<<endl;} 278 256 hasemitted = true; 279 257 energy = energy - gee; … … 330 308 } 331 309 332 std::ostream& operator<< (std::ostream& os, const H_BeamParticle& p){333 os << " M = " << p.getM() << "GeV ";334 os << " Q = " << p.getQ() << "e";335 os << " fx = " << p.getX() << "m ";336 os << " fy = " << p.getY() << "m ";337 os << " thx = " << p.getTX() << "rad ";338 os << " thy = " << p.getTY() << "rad ";339 os<< endl;340 return os;310 void H_BeamParticle::printProperties() const { 311 cout << " M = " << getM() << "GeV "; 312 cout << " Q = " << getQ() << "e"; 313 cout << " fx = " << getX() << "m "; 314 cout << " fy = " << getY() << "m "; 315 cout << " thx = " << getTX() << "rad "; 316 cout << " thy = " << getTY() << "rad "; 317 cout << endl; 318 return; 341 319 } 342 320 … … 378 356 if((*position_i)[INDEX_S]>=position) { 379 357 if(position_i==positions.begin()) { 380 if(VERBOSE) cout<<" <H_BeamParticle>ERROR : non reachable value"<<endl;358 if(VERBOSE) cout<<"ERROR : non reachable value"<<endl; 381 359 return; 382 360 } 383 361 l = (*position_i)[INDEX_S] - (*(position_i-1))[INDEX_S]; 384 362 if(l==0) { 385 if(VERBOSE) cout<<" <H_BeamParticle>WARNING : no luck in choosing position, no propagation done"<<endl;363 if(VERBOSE) cout<<"WARNING : no luck in choosing position, no propagation done"<<endl; 386 364 return; 387 365 } … … 396 374 position_i = positions.begin(); 397 375 cout << "Desired position is : " << position << " & positions.begin() is " << (*position_i)[INDEX_S] << endl; 398 cout<<" <H_BeamParticle>ERROR : position not reachable"<<endl;376 cout<<"ERROR : position not reachable"<<endl; 399 377 return; 400 378 } 401 379 } 402 380 403 /// Caution : do not use this method (obsolete)!!!381 /// Caution : do not use this method !!! 404 382 void H_BeamParticle::propagate(const H_AbstractBeamLine * beam, const H_OpticalElement * element) { 405 383 TMatrixD X(*getV()); 406 X *= beam->getPartialMatrix(element);384 X *= *(beam->getPartialMatrix(element)); 407 385 fx = URAD*(X.GetMatrixArray())[0]; 408 386 thx = URAD*atan((X.GetMatrixArray())[1]); … … 412 390 } 413 391 414 void H_BeamParticle::propagate(const H_AbstractBeamLine * beam, const string &el_name) {415 propagate(beam->getElement(el_name)->getS() +beam->getElement(el_name)->getLength());392 void H_BeamParticle::propagate(const H_AbstractBeamLine * beam, const string el_name) { 393 propagate(beam->getElement(el_name)->getS()); 416 394 return; 417 395 } … … 419 397 void H_BeamParticle::propagate(const H_AbstractBeamLine * beam) { 420 398 TMatrixD X(*getV()); 421 X *= beam->getBeamMatrix();399 X *= *(beam->getBeamMatrix()); 422 400 fx = URAD*(X.GetMatrixArray())[0]; 423 401 thx = URAD*atan((X.GetMatrixArray())[1]); … … 436 414 return ; 437 415 } 438 416 /* 439 417 TGraph * H_BeamParticle::getPath(const int x_or_y, const int color) const{ 440 418 /// @param x_or_y = 0(1) draws the x(y) component; 441 419 442 420 const int N = (int) positions.size(); 443 444 if(N<2) cout<<" <H_BeamParticle> WARNING :particle positions not calculated : please run computePath"<<endl;421 int mycolor = color; 422 if(N<2) cout<<"particle positions not calculated : please run computePath"<<endl; 445 423 double * s = new double[N], * graph = new double[N]; 446 424 … … 461 439 } 462 440 463 void H_BeamParticle::getPath(const int x_or_y, const string& filename) const{ 464 /// @param x_or_y = 0(1) draws the x(y) component; 465 ofstream outfile(filename.c_str()); 466 467 int index; 468 if(x_or_y==0) {index = INDEX_X;} else {index = INDEX_Y;} 441 TPolyLine3D * H_BeamParticle::getPath3D(const H_AbstractBeamLine * beam, const bool isfirst, const int color, const int side) const{ 442 const int N = (int) positions.size(); 443 int mycolor = color; 444 if(N<2) cout<<"WARNING : particle positions not calculated. Run computePath"<<endl; 445 double * s = new double[N], * graphx = new double[N], * graphy = new double[N]; 446 int direction = (side<0)?-1:1; 447 469 448 470 449 vector<TVectorD>::const_iterator position_i; 471 450 for(position_i = positions.begin(); position_i < positions.end(); position_i++) { 472 outfile << (*position_i)[index] << "\t" << (*position_i)[INDEX_S] << endl; 451 graphx[(int)(position_i-positions.begin())] = (*position_i)[INDEX_X]; 452 graphy[(int)(position_i-positions.begin())] = (*position_i)[INDEX_Y]; 453 s[(int)(position_i-positions.begin())] = (*position_i)[INDEX_S]*1000*direction; 473 454 } 474 outfile.close(); 475 } 476 455 456 float coi[3] = {beam->getLength()*(-1000),-10000,-5000}; 457 float cof[3] = {beam->getLength()*1000,10000,5000}; 458 TView *view = TView::CreateView(11); 459 view->SetRange(coi[0],coi[1],coi[2],cof[0],cof[1],cof[2]); 460 TPolyLine3D* ppath = new TPolyLine3D(N,s,graphx,graphy); 461 ppath->SetLineColor(mycolor); 462 463 if(isfirst) { 464 ppath->Draw(); 465 view->ShowAxis(); 466 } else { 467 ppath->Draw("same"); 468 } 469 470 delete [] s; 471 delete [] graphx; 472 delete [] graphy; 473 return ppath; 474 } // getPath3D 475 */ 476 void H_BeamParticle::computePath(const H_AbstractBeamLine * beam) { 477 computePath(beam,true); 478 } 477 479 478 480 // should be removed later, to keep only computePath(const H_AbstractBeamLine & , const bool) … … 505 507 for (int i=0; i<N; i++) { 506 508 const unsigned pos = i; 507 if(fs < beam->getElement(pos)->getS() && fs > beam->getElement(pos-1)->getS()) { 508 if(beam->getElement(pos-1)->getType()!=DRIFT) { 509 cout<<"Path starts inside element "<<beam->getElement(pos-1)->getName()<<endl; 510 } else { 511 cout<<"Path starts inside unnamed drift "<<endl; 512 } 513 H_OpticalElement* temp_el = beam->getElement(pos-1)->clone(); 514 temp_el->setS(fs); 515 temp_el->setLength(beam->getElement(pos)->getS() - fs); 516 mat[0][0] = mat[0][0] - temp_el->getX(); 517 mat[0][1] = mat[0][1] - tan(temp_el->getTX()); 518 mat[0][2] = mat[0][2] - temp_el->getY(); 519 mat[0][3] = mat[0][3] - tan(temp_el->getTY()); 520 mat *= temp_el->getMatrix(energy_loss,mp,qp); 521 mat[0][0] = mat[0][0] + temp_el->getX(); 522 mat[0][1] = mat[0][1] + tan(temp_el->getTX()); 523 mat[0][2] = mat[0][2] + temp_el->getY(); 524 mat[0][3] = mat[0][3] + tan(temp_el->getTY()); 525 xys[0] = mat.GetMatrixArray()[0]*URAD; 526 xys[1] = atan(mat.GetMatrixArray()[1])*URAD; 527 xys[2] = mat.GetMatrixArray()[2]*URAD; 528 xys[3] = atan(mat.GetMatrixArray()[3])*URAD; 529 xys[4] = temp_el->getS()+temp_el->getLength(); 530 addPosition(xys[0],xys[1],xys[2],xys[3],xys[4]); 531 if (temp_el) delete temp_el; 532 } 533 if(fs <= beam->getElement(pos)->getS()) { 534 mat[0][0] = mat[0][0] - beam->getElement(pos)->getX(); 535 mat[0][1] = mat[0][1] - tan(beam->getElement(pos)->getTX()/URAD)*URAD; 536 mat[0][2] = mat[0][2] - beam->getElement(pos)->getY(); 537 mat[0][3] = mat[0][3] - tan(beam->getElement(pos)->getTY()/URAD)*URAD; 538 mat *= beam->getElement(pos)->getMatrix(energy_loss,mp,qp); 539 mat[0][0] = mat[0][0] + beam->getElement(pos)->getX(); 540 mat[0][1] = mat[0][1] + tan(beam->getElement(pos)->getTX()/URAD)*URAD; 541 mat[0][2] = mat[0][2] + beam->getElement(pos)->getY(); 542 mat[0][3] = mat[0][3] + tan(beam->getElement(pos)->getTY()/URAD)*URAD; 543 xys[0] = mat.GetMatrixArray()[0]*URAD; 544 xys[1] = atan(mat.GetMatrixArray()[1])*URAD; 545 xys[2] = mat.GetMatrixArray()[2]*URAD; 546 xys[3] = atan(mat.GetMatrixArray()[3])*URAD; 547 xys[4] = beam->getElement(pos)->getS()+beam->getElement(pos)->getLength(); 548 addPosition(xys[0],xys[1],xys[2],xys[3],xys[4]); 549 } 550 } 551 } 552 553 // part about non-ip particle is not ready yet. use the above method in the meantime 509 mat[0][0] = mat[0][0] - beam->getElement(pos)->getX(); 510 mat[0][1] = mat[0][1] - tan(beam->getElement(pos)->getTX()/URAD)*URAD; 511 mat[0][2] = mat[0][2] - beam->getElement(pos)->getY(); 512 mat[0][3] = mat[0][3] - tan(beam->getElement(pos)->getTY()/URAD)*URAD; 513 mat *= beam->getElement(pos)->getMatrix(energy_loss,mp,qp); 514 mat[0][0] = mat[0][0] + beam->getElement(pos)->getX(); 515 mat[0][1] = mat[0][1] + tan(beam->getElement(pos)->getTX()/URAD)*URAD; 516 mat[0][2] = mat[0][2] + beam->getElement(pos)->getY(); 517 mat[0][3] = mat[0][3] + tan(beam->getElement(pos)->getTY()/URAD)*URAD; 518 xys[0] = mat.GetMatrixArray()[0]*URAD; 519 xys[1] = atan(mat.GetMatrixArray()[1])*URAD; 520 xys[2] = mat.GetMatrixArray()[2]*URAD; 521 xys[3] = atan(mat.GetMatrixArray()[3])*URAD; 522 xys[4] = beam->getElement(pos)->getS()+beam->getElement(pos)->getLength(); 523 addPosition(xys[0],xys[1],xys[2],xys[3],xys[4]); 524 fx = xys[0]; 525 fy = xys[2]; 526 thx = xys[1]; 527 thy = xys[3]; 528 } 529 } 530 554 531 void H_BeamParticle::computePath(const H_AbstractBeamLine & beam, const bool NonLinear) { 555 532 TMatrixD temp_mat(MDIM,MDIM); … … 578 555 double energy_loss = NonLinear?BE-energy:0; 579 556 580 // modify here to allow starting at non-IP positions581 // s is distance to IP582 // initial position is already in positions vector ?583 557 for (int i=0; i<N; i++) { 584 558 const unsigned pos = i; 585 // if we are inside an element, we should start by adding the action 586 // of the rest of this element 587 if(pos > 0 && fs < beam.getElement(pos)->getS() && fs > beam.getElement(pos-1)->getS()) { 588 cout<<"Path starts inside element "<<beam.getElement(pos-1)->getName()<<endl; 589 H_OpticalElement* temp_el = new H_OpticalElement(*(beam.getElement(pos-1))); 590 temp_el->setS(fs); 591 temp_el->setLength(beam.getElement(pos)->getS() - fs); 592 mat[0][0] = mat[0][0] - temp_el->getX(); 593 mat[0][1] = mat[0][1] - tan(temp_el->getTX()); 594 mat[0][2] = mat[0][2] - temp_el->getY(); 595 mat[0][3] = mat[0][3] - tan(temp_el->getTY()); 596 mat *= temp_el->getMatrix(energy_loss,mp,qp); 597 mat[0][0] = mat[0][0] + temp_el->getX(); 598 mat[0][1] = mat[0][1] + tan(temp_el->getTX()); 599 mat[0][2] = mat[0][2] + temp_el->getY(); 600 mat[0][3] = mat[0][3] + tan(temp_el->getTY()); 601 } else if(fs >= beam.getElement(pos)->getS()) { 602 mat[0][0] = mat[0][0] - beam.getElement(pos)->getX(); 603 mat[0][1] = mat[0][1] - tan(beam.getElement(pos)->getTX()); 604 mat[0][2] = mat[0][2] - beam.getElement(pos)->getY(); 605 mat[0][3] = mat[0][3] - tan(beam.getElement(pos)->getTY()); 606 mat *= beam.getElement(pos)->getMatrix(energy_loss,mp,qp); 607 mat[0][0] = mat[0][0] + beam.getElement(pos)->getX(); 608 mat[0][1] = mat[0][1] + tan(beam.getElement(pos)->getTX()); 609 mat[0][2] = mat[0][2] + beam.getElement(pos)->getY(); 610 mat[0][3] = mat[0][3] + tan(beam.getElement(pos)->getTY()); 611 } 612 xys[0] = mat.GetMatrixArray()[0]*URAD; 613 xys[1] = atan(mat.GetMatrixArray()[1])*URAD; 614 xys[2] = mat.GetMatrixArray()[2]*URAD; 615 xys[3] = atan(mat.GetMatrixArray()[3])*URAD; 616 xys[4] = beam.getElement(pos)->getS()+beam.getElement(pos)->getLength(); 617 addPosition(xys[0],xys[1],xys[2],xys[3],xys[4]); 618 fx = xys[0]; 619 fy = xys[2]; 620 thx = xys[1]; 621 thy = xys[3]; 622 } 559 mat[0][0] = mat[0][0] - beam.getElement(pos)->getX(); 560 mat[0][1] = mat[0][1] - tan(beam.getElement(pos)->getTX()); 561 mat[0][2] = mat[0][2] - beam.getElement(pos)->getY(); 562 mat[0][3] = mat[0][3] - tan(beam.getElement(pos)->getTY()); 563 mat *= beam.getElement(pos)->getMatrix(energy_loss,mp,qp); 564 mat[0][0] = mat[0][0] + beam.getElement(pos)->getX(); 565 mat[0][1] = mat[0][1] + tan(beam.getElement(pos)->getTX()); 566 mat[0][2] = mat[0][2] + beam.getElement(pos)->getY(); 567 mat[0][3] = mat[0][3] + tan(beam.getElement(pos)->getTY()); 568 xys[0] = mat.GetMatrixArray()[0]*URAD; 569 xys[1] = atan(mat.GetMatrixArray()[1])*URAD; 570 xys[2] = mat.GetMatrixArray()[2]*URAD; 571 xys[3] = atan(mat.GetMatrixArray()[3])*URAD; 572 xys[4] = beam.getElement(pos)->getS()+beam.getElement(pos)->getLength(); 573 addPosition(xys[0],xys[1],xys[2],xys[3],xys[4]); 574 fx = xys[0]; 575 fy = xys[2]; 576 thx = xys[1]; 577 thy = xys[3]; 578 } 623 579 } 624 580 -
trunk/external/Hector/H_BeamParticle.h
r1360 r1365 2 2 #define _H_BeamParticle_ 3 3 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 5 * * 6 * --<--<-- A fast simulator --<--<-- * 7 * / --<--<-- of particle --<--<-- * 8 * ----HECTOR----< * 9 * \ -->-->-- transport through -->-->-- * 10 * -->-->-- generic beamlines -->-->-- * 11 * * 12 * JINST 2:P09005 (2007) * 13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 14 * http://www.fynu.ucl.ac.be/hector.html * 15 * * 16 * Center for Cosmology, Particle Physics and Phenomenology * 17 * Universite catholique de Louvain * 18 * Louvain-la-Neuve, Belgium * 19 * * 20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 4 /* 5 ---- Hector the simulator ---- 6 A fast simulator of particles through generic beamlines. 7 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 8 9 http://www.fynu.ucl.ac.be/hector.html 10 11 Centre de Physique des Particules et de Phénoménologie (CP3) 12 Université Catholique de Louvain (UCL) 13 */ 21 14 22 15 /// \file H_BeamParticle.h … … 33 26 #include "TMatrixD.h" 34 27 #include "TVectorD.h" 35 #include "TLorentzVector.h"36 #include "TRandom 3.h"28 //#include "TPolyLine3D.h" 29 #include "TRandom.h" 37 30 38 31 // local #includes … … 43 36 using namespace std; 44 37 38 // local defines 39 #define LENGTH_VEC 5 40 #define INDEX_X 0 41 #define INDEX_TX 1 42 #define INDEX_Y 2 43 #define INDEX_TY 3 44 #define INDEX_S 4 45 // (x,theta_x,y,theta_y,s) 46 45 47 /// Defines a particle from the beam and its transport through the beamline 46 48 class H_BeamParticle { 47 49 48 50 public: 51 void init(); 49 52 /// Constructors and Destructor 50 53 //@{ … … 53 56 H_BeamParticle(const double, const double); 54 57 H_BeamParticle& operator=(const H_BeamParticle&); 55 ~H_BeamParticle() { if(stop_position) delete stop_position; if(!stop_element) delete stop_element; positions.clear();}58 ~H_BeamParticle() {delete stop_position; if(!stop_element) delete stop_element; positions.clear(); return; } 56 59 //@} 57 60 /// Smears the (x,y) coordinates of the particle [\f$ \mu m \f$] … … 67 70 /// Sets the particle 4-momentum \f$ P^\mu \f$ 68 71 void set4Momentum(const double, const double, const double, const double); 69 /// Sets the particle 4-momentum \f$ P^\mu \f$70 void set4Momentum(const TLorentzVector&);71 72 /// Clears H_BeamParticle::positions and sets the initial one. 72 73 void setPosition(const double , const double , const double , const double , const double ); 73 74 /// Returns the particle mass [GeV] 74 constdouble getM() const {return mp;};75 double getM() const {return mp;}; 75 76 /// Returns the particle charge [e] 76 constdouble getQ() const {return qp;};77 double getQ() const {return qp;}; 77 78 /// Returns the current x coordinate [\f$ \mu \f$m] 78 constdouble getX() const {return fx;};79 double getX() const {return fx;}; 79 80 /// Returns the current y coordinate [\f$ \mu \f$m] 80 constdouble getY() const {return fy;};81 double getY() const {return fy;}; 81 82 /// Returns the current s coordinate [m] 82 constdouble getS() const {return fs;};83 inline double getS() const {return fs;}; 83 84 /// Returns the current \f$ \theta_x \f$ angular coordinate [\f$ \mu \f$rad] 84 constdouble getTX() const {return thx;};85 inline double getTX() const {return thx;}; 85 86 /// Returns the current \f$ \theta_y \f$ angular coordinate [\f$ \mu \f$rad] 86 constdouble getTY() const {return thy;};87 inline double getTY() const {return thy;}; 87 88 /// Returns the current particle energy [GeV] 88 constdouble getE() const {return energy;};89 inline double getE() const {return energy;}; 89 90 /// Returns all the positions 90 91 vector<TVectorD> getPositions() const {return positions;}; 91 constbool isPhysical() const {return isphysical;};92 bool isPhysical() const {return isphysical;}; 92 93 /// \brief Simulates the emission of a photon in a random direction 93 94 /// … … 98 99 /// \f$ Q^{2}_{max} = -2 * \big( \frac{M_{p} E_{\gamma}}{p_{1}+p_{2}} \big) \big[ 1 + \frac{E^{2}_{1} + E^{2}_{2} - M^{2}_{p} }{ E_{1} E_{2} + p_{1} p_{2}} \big] \f$ 99 100 //@{ 100 void emitGamma(const double gee, const double gq2, const double phimin=0, const double phimax=2*TMath::Pi()); 101 void emitGamma(const double, const double, const double, const double); 102 void emitGamma(const double, const double); 101 103 //@} 102 104 /// uses Pythia to generate some inelastic pp->pX collision as background … … 110 112 void propagate(const H_AbstractBeamLine *, const H_OpticalElement *); 111 113 /// Propagates the particle accross the beamline until a given element 112 void propagate(const H_AbstractBeamLine *, const string &);114 void propagate(const H_AbstractBeamLine *, const string); 113 115 /// Propagates the particle until the end of the beamline 114 116 void propagate(const H_AbstractBeamLine *); … … 118 120 const TVectorD * getPosition(const int ) const; 119 121 /// Prints the properties of the particle 120 void printProperties() const {cout << *this; return;};122 void printProperties() const; 121 123 /// Prints the phase vector of the particle 122 124 void printV() const; … … 124 126 const H_OpticalElement * getStoppingElement() const; 125 127 /// Checks if the particle has been stopped in any element of the beamline 126 constbool stopped(const H_AbstractBeamLine *);128 bool stopped(const H_AbstractBeamLine *); 127 129 /// Returns the StopPosition vector 128 130 inline const TVectorD * getStopPosition() const { return stop_position; }; … … 131 133 void showPositions() const; 132 134 /// Returns the particle path in the beamline 133 TGraph * getPath(const int , const int ) const; 134 void getPath(const int x_or_y, const string& filename) const; 135 ////TGraph * getPath(const int , const int ) const; 136 /// Draws the particle path in the beamline in 3D 137 ////TPolyLine3D * getPath3D(const H_AbstractBeamLine *, const bool, const int, const int) const; 138 /// Computes the position of the particle at the end of each element of the beam, without non linear effects 139 void computePath(const H_AbstractBeamLine *); 135 140 /// Computes the position of the particle at the end of each element of the beam. 136 void computePath(const H_AbstractBeamLine * beam, const bool NonLinear=true);141 void computePath(const H_AbstractBeamLine *, const bool); 137 142 /// Computes the position of the particle at the end of each element of the beam. 138 void computePath(const H_AbstractBeamLine & beam, const bool NonLinear=true);143 void computePath(const H_AbstractBeamLine &, const bool); 139 144 /// Clears H_BeamParticle::positions but keeps the initial vector. 140 145 void resetPath(); 141 146 142 147 private: 143 void init();144 148 /// Particle mass [GeV] 145 149 double mp; … … 172 176 /// Adds a new vector (x,tx,y,ty,s) at the end of H_BeamParticle::positions 173 177 void addPosition(const double , const double , const double , const double , const double ); 174 175 friend std::ostream& operator<< (std::ostream& os, const H_BeamParticle& p);176 178 }; 177 179 #endif -
trunk/external/Hector/H_CircularAperture.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 18 11 19 12 /// \file H_CircularAperture.cc … … 30 23 using namespace std; 31 24 32 H_CircularAperture* H_CircularAperture::clone() const { 33 return new H_CircularAperture(x1,fx,fy); 25 /// Circular apertures 26 H_CircularAperture::H_CircularAperture(const float r, const float posx, const float posy) :H_EllipticAperture(r,r,posx,posy) { 27 /// @param r is the radius of the circular shape 28 /// @param posx, posy are the (x,y) coordinates of the center of the circle 29 type= CIRCULAR; 34 30 } 35 31 36 std::ostream& operator<< (std::ostream& os, const H_CircularAperture& ap){37 os << "Aperture shape:" << ap.aptypestring << ", aperture radius : " << ap.x1 << endl;38 os << " \t Center : " << ap.fx << "," << ap.fy << endl;39 return os;32 void H_CircularAperture::printProperties() const { 33 cout << "Aperture shape:" << getTypeString() << ", aperture radius : " << x1 << endl; 34 cout << " \t Center : " << fx << "," << fy << endl; 35 return; 40 36 } -
trunk/external/Hector/H_CircularAperture.h
r1360 r1365 2 2 #define _H_CircularAperture_ 3 3 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 5 * * 6 * --<--<-- A fast simulator --<--<-- * 7 * / --<--<-- of particle --<--<-- * 8 * ----HECTOR----< * 9 * \ -->-->-- transport through -->-->-- * 10 * -->-->-- generic beamlines -->-->-- * 11 * * 12 * JINST 2:P09005 (2007) * 13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 14 * http://www.fynu.ucl.ac.be/hector.html * 15 * * 16 * Center for Cosmology, Particle Physics and Phenomenology * 17 * Universite catholique de Louvain * 18 * Louvain-la-Neuve, Belgium * 19 * * 20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 4 /* 5 ---- Hector the simulator ---- 6 A fast simulator of particles through generic beamlines. 7 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 8 9 http://www.fynu.ucl.ac.be/hector.html 10 11 Centre de Physique des Particules et de Phénoménologie (CP3) 12 Université Catholique de Louvain (UCL) 13 */ 21 14 22 15 /// \file H_CircularAperture.h … … 33 26 /// Constructors and destructor 34 27 //@{ 35 H_CircularAperture():H_EllipticAperture(CIRCULAR,0,0,0,0) {}; 36 H_CircularAperture(const float r, const float posx, const float posy) : H_EllipticAperture(CIRCULAR,r,r,posx,posy) {}; 37 /// @param r is the radius of the circular shape 38 /// @param posx, posy are the (x,y) coordinates of the center of the circle 39 ~H_CircularAperture() {}; 40 H_CircularAperture* clone() const; 28 H_CircularAperture():H_EllipticAperture(0,0,0,0) {type = CIRCULAR; setApertureString();} 29 H_CircularAperture(const float, const float, const float); 30 ~H_CircularAperture() {return;}; 41 31 //@} 42 friend std::ostream& operator<< (std::ostream& os, const H_CircularAperture& ap);32 virtual void printProperties() const; 43 33 }; 44 34 -
trunk/external/Hector/H_Dipole.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 18 11 19 12 /// \file H_Dipole.cc … … 26 19 // needed for call from R- and S-Dipoles constructor 27 20 // must be in public section 28 element_mat.ResizeTo(MDIM,MDIM);29 21 setTypeString(); 30 22 if (fk !=0 ) setMatrix(0,MP,QP); … … 32 24 } 33 25 34 std::ostream& operator<< (std::ostream& os, const H_Dipole& el) { 35 os << el.typestring << el.name <<"\t at s = "<< el.fs <<"\t length = "<< el.element_length <<"\t k0 = "<<el.fk << endl; 36 if(el.element_aperture->getType()!=NONE) { 37 os << *(el.element_aperture) << endl; 26 void H_Dipole::printProperties() const { 27 cout << typestring; 28 cout << name; 29 cout<<"\t at s = "<< fs; 30 cout<<"\t length = "<< element_length; 31 cout<<"\t k0 = "<<fk; 32 cout<<endl; 33 if(element_aperture->getType()!=NONE) { 34 cout <<"\t aperture type = " << element_aperture->getTypeString(); 35 element_aperture->printProperties(); 38 36 } 39 37 40 if(el.element_length<0) { if(VERBOSE) os <<"<H_Dipole> ERROR : Interpenetration of elements !"<<endl; } 41 if(el.element_length==0) { if(VERBOSE) os <<"<H_Dipole> WARNING : 0-length "<< el.name << " !" << endl; } 42 return os; 38 if(element_length<0) { if(VERBOSE) cout<<"\t ERROR : Interpenetration of elements !"<<endl; } 39 if(element_length==0) { if(VERBOSE) cout<<"\t WARNING : 0-length "<< name << " !" << endl; } 43 40 } -
trunk/external/Hector/H_Dipole.h
r1360 r1365 2 2 /// \brief Class aiming at simulating LHC beam dipoles. 3 3 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 5 * * 6 * --<--<-- A fast simulator --<--<-- * 7 * / --<--<-- of particle --<--<-- * 8 * ----HECTOR----< * 9 * \ -->-->-- transport through -->-->-- * 10 * -->-->-- generic beamlines -->-->-- * 11 * * 12 * JINST 2:P09005 (2007) * 13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 14 * http://www.fynu.ucl.ac.be/hector.html * 15 * * 16 * Center for Cosmology, Particle Physics and Phenomenology * 17 * Universite catholique de Louvain * 18 * Louvain-la-Neuve, Belgium * 19 * * 20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 4 /* 5 ---- Hector the simulator ---- 6 A fast simulator of particles through generic beamlines. 7 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 8 9 http://www.fynu.ucl.ac.be/hector.html 10 11 Centre de Physique des Particules et de Phénoménologie (CP3) 12 Université Catholique de Louvain (UCL) 13 */ 21 14 22 15 #ifndef _H_Dipole_ … … 33 26 H_Dipole():H_OpticalElement() {} 34 27 H_Dipole(const int dtype, const double s, const double k, const double l):H_OpticalElement(dtype,s,k,l){} 35 H_Dipole(const string &nameE, const int dtype, const double s, const double k, const double l):H_OpticalElement(nameE,dtype,s,k,l){}36 virtual ~H_Dipole() { };28 H_Dipole(const string nameE, const int dtype, const double s, const double k, const double l):H_OpticalElement(nameE,dtype,s,k,l){} 29 virtual ~H_Dipole() {return;}; 37 30 //@} 38 31 /// Prints the properties of the element 39 virtual void printProperties() const { cout << *this; return;}; 40 virtual H_Dipole* clone() const =0; 41 void init(); 32 virtual void printProperties() const; 33 void init(); 42 34 43 35 protected: 44 36 virtual void setTypeString() =0; 45 virtual void setMatrix(const float, const float, const float) = 0; 46 47 friend std::ostream& operator<< (std::ostream& os, const H_Dipole& el); 37 virtual void setMatrix(const float, const float, const float) const =0; 38 48 39 }; 49 40 -
trunk/external/Hector/H_Drift.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 18 11 19 12 /// \file H_Drift.cc … … 24 17 #include "H_TransportMatrices.h" 25 18 26 27 19 void H_Drift::init() { 28 20 // must be in public section 29 element_mat.ResizeTo(MDIM,MDIM);30 21 setTypeString(); 31 22 setMatrix(0,MP,QP); … … 33 24 } 34 25 35 std::ostream& operator<< (std::ostream& os, const H_Drift& el) { 36 os << el.typestring << el.name << "\t\t at s = " << el.fs << "\t length = " << el.element_length <<endl; 37 if(el.element_aperture->getType()!=NONE) { 38 os << *(el.element_aperture) << endl; 26 void H_Drift::printProperties() const { 27 cout << typestring << name; 28 cout << "\t\t at s = " << fs; 29 cout << "\t length = " << element_length; 30 cout<<endl; 31 if(element_aperture->getType()!=NONE) { 32 cout <<"\t aperture type = " << element_aperture->getTypeString(); 33 element_aperture->printProperties(); 39 34 } 40 if(el.element_length<0 && VERBOSE) os <<"<H_Drift> ERROR : Interpenetration of elements !"<<endl; 41 else if(el.element_length==0 && VERBOSE) os <<"<H_Drift> WARNING : 0-length "<< el.name << " !" << endl; 42 return os; 35 if(element_length<0) { if(VERBOSE) cout<<"\t ERROR : Interpenetration of elements !"<<endl; } 36 if(element_length==0) { if(VERBOSE) cout<<"\t WARNING : 0-length "<< name << " !" << endl; } 43 37 } 44 38 45 //void H_Drift::setMatrix(const float eloss, const float p_mass, const float p_charge) { 46 void H_Drift::setMatrix(const float , const float , const float ) { 47 element_mat = driftmat(element_length); 39 void H_Drift::setMatrix(const float eloss, const float p_mass, const float p_charge) const { 40 *element_mat = driftmat(element_length); 48 41 return ; 49 42 } 50 51 H_Drift* H_Drift::clone() const {52 H_Drift* temp_drift = new H_Drift(name,fs,element_length);53 temp_drift->setX(xpos);54 temp_drift->setY(ypos);55 temp_drift->setTX(txpos);56 temp_drift->setTY(typos);57 temp_drift->setBetaX(betax);58 temp_drift->setBetaY(betay);59 return temp_drift;60 }61 -
trunk/external/Hector/H_Drift.h
r1360 r1365 2 2 #define _H_Drift_ 3 3 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 5 * * 6 * --<--<-- A fast simulator --<--<-- * 7 * / --<--<-- of particle --<--<-- * 8 * ----HECTOR----< * 9 * \ -->-->-- transport through -->-->-- * 10 * -->-->-- generic beamlines -->-->-- * 11 * * 12 * JINST 2:P09005 (2007) * 13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 14 * http://www.fynu.ucl.ac.be/hector.html * 15 * * 16 * Center for Cosmology, Particle Physics and Phenomenology * 17 * Universite catholique de Louvain * 18 * Louvain-la-Neuve, Belgium * 19 * * 20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 4 /* 5 ---- Hector the simulator ---- 6 A fast simulator of particles through generic beamlines. 7 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 8 9 http://www.fynu.ucl.ac.be/hector.html 10 11 Centre de Physique des Particules et de Phénoménologie (CP3) 12 Université Catholique de Louvain (UCL) 13 */ 21 14 22 15 /// \file H_Drift.h … … 34 27 H_Drift():H_OpticalElement(DRIFT,0.,0.,0.) {init();} 35 28 H_Drift(const double s, const double l):H_OpticalElement(DRIFT,s,0.,l){init();} 36 H_Drift(const string &nameE, const double s, const double l):H_OpticalElement(nameE,DRIFT,s,0.,l){init();}37 ~H_Drift() { };29 H_Drift(const string nameE, const double s, const double l):H_OpticalElement(nameE,DRIFT,s,0.,l){init();} 30 ~H_Drift() { return; }; 38 31 //@} 32 virtual void printProperties() const; 39 33 void init(); 40 virtual void printProperties() const { cout << *this; return;};41 H_Drift* clone() const;42 34 43 pr otected:35 private: 44 36 virtual void setTypeString() {typestring = DRIFTNAME;}; 45 virtual void setMatrix(const float, const float, const float) ; 46 friend std::ostream& operator<< (std::ostream& os, const H_Drift& el); 37 virtual void setMatrix(const float, const float, const float) const ; 47 38 }; 48 39 -
trunk/external/Hector/H_EllipticAperture.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 18 11 19 12 /// \file H_EllipticAperture.cc … … 24 17 25 18 // ROOT #includes 26 #include "TEllipse.h"19 //#include "TEllipse.h" 27 20 28 21 // local #includes … … 32 25 using namespace std; 33 26 34 H_EllipticAperture::H_EllipticAperture(const int type, const float l, const float h, const float posx, const float posy) : 35 H_Aperture(type,l,h,0,0,posx,posy) { 36 /// @param type is the aperture type (ELLIPTIC or CIRCULAR) 37 /// @param l, h are the length and height of the elliptic shape 27 H_EllipticAperture::H_EllipticAperture(const float x1, const float x2, const float posx, const float posy) :H_Aperture(ELLIPTIC,x1,x2,0,0,posx,posy) { 28 /// @param x1, x2 are the length and height of the elliptic shape 38 29 /// @param posx, posy are the (x,y) coordinates of the center of the ellipse 39 40 if (type!= ELLIPTIC && type != CIRCULAR) {41 cout << "Warning: trying to define an EllipticalAperture which is neither elliptical nor circular." << endl;42 cout << "'Elliptical' type forced\n";43 type_ = ELLIPTIC;44 aptypestring = getApertureString();45 }46 }47 48 H_EllipticAperture::H_EllipticAperture(const float l, const float h, const float posx, const float posy) :49 H_Aperture(ELLIPTIC,l,h,0,0,posx,posy) {50 /// @param l, h are the length and height of the elliptic shape51 /// @param posx, posy are the (x,y) coordinates of the center of the ellipse52 }53 54 H_EllipticAperture* H_EllipticAperture::clone() const {55 return new H_EllipticAperture(x1,x2,fx,fy);56 30 } 57 31 … … 61 35 } 62 36 63 void H_EllipticAperture::draw( const float scale) const {64 TEllipse* te = new TEllipse(fx*scale,fy*scale,x1*scale,x2*scale);37 void H_EllipticAperture::draw() const { 38 /* Ellipse* te = new TEllipse(fx,fy,x1,x2); 65 39 te->SetFillStyle(3003); 66 te->SetLineColor(39); 67 te->SetFillColor(39); 68 te->Draw("f"); 40 te->SetLineColor(2); 41 te->SetFillColor(2); 42 te->Draw(); 43 return; 44 */ 45 } 46 47 void H_EllipticAperture::printProperties() const { 48 cout<< "Aperture shape:" << getTypeString() << ", ellipse axes : "<<x1<<", "<<x2<<endl; 49 cout << " \t Center : " << fx << "," << fy << endl; 69 50 return; 70 51 } 71 72 std::ostream& operator<< (std::ostream& os, const H_EllipticAperture& ap) {73 os<< "Aperture shape:" << ap.aptypestring << ", ellipse axes : "<< ap.x1 <<", " << ap.x2 << endl;74 os << " \t Center : " << ap.fx << "," << ap.fy << endl;75 return os;76 } -
trunk/external/Hector/H_EllipticAperture.h
r1360 r1365 2 2 #define _H_EllipticAperture_ 3 3 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 5 * * 6 * --<--<-- A fast simulator --<--<-- * 7 * / --<--<-- of particle --<--<-- * 8 * ----HECTOR----< * 9 * \ -->-->-- transport through -->-->-- * 10 * -->-->-- generic beamlines -->-->-- * 11 * * 12 * JINST 2:P09005 (2007) * 13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 14 * http://www.fynu.ucl.ac.be/hector.html * 15 * * 16 * Center for Cosmology, Particle Physics and Phenomenology * 17 * Universite catholique de Louvain * 18 * Louvain-la-Neuve, Belgium * 19 * * 20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 4 /* 5 ---- Hector the simulator ---- 6 A fast simulator of particles through generic beamlines. 7 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 8 9 http://www.fynu.ucl.ac.be/hector.html 10 11 Centre de Physique des Particules et de Phénoménologie (CP3) 12 Université Catholique de Louvain (UCL) 13 */ 21 14 22 15 /// \file H_EllipticAperture.h … … 33 26 //@{ 34 27 H_EllipticAperture():H_Aperture(ELLIPTIC,0,0,0,0,0,0) {} 35 H_EllipticAperture(const int, const float, const float, const float, const float);36 28 H_EllipticAperture(const float, const float, const float, const float); 37 ~H_EllipticAperture() {}; 38 virtual H_EllipticAperture* clone() const; 29 ~H_EllipticAperture() {return;}; 39 30 //@} 40 31 /// Checks whether the point is inside the aperture or not 41 32 virtual bool isInside(const float, const float) const; 42 33 /// Draws the aperture shape. 43 virtual void draw( const float scale=1) const;44 friend std::ostream& operator<< (std::ostream& os, const H_EllipticAperture& ap);34 virtual void draw() const; 35 virtual void printProperties() const; 45 36 }; 46 37 -
trunk/external/Hector/H_HorizontalKicker.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 18 11 19 12 /// \file H_HorizontalKicker.cc … … 28 21 int kickers_on = 0; 29 22 30 void H_HorizontalKicker::setMatrix(const float eloss, const float p_mass, const float p_charge) {23 void H_HorizontalKicker::setMatrix(const float eloss, const float p_mass, const float p_charge) const { 31 24 extern int kickers_on; 32 25 if(kickers_on) { 33 element_mat = hkickmat(element_length,fk,eloss,p_mass,p_charge);26 *element_mat = hkickmat(element_length,fk,eloss,p_mass,p_charge); 34 27 } else { 35 element_mat = driftmat(element_length);28 *element_mat = driftmat(element_length); 36 29 } 37 30 return ; 38 31 } 39 40 H_HorizontalKicker* H_HorizontalKicker::clone() const {41 H_HorizontalKicker* temp_kick = new H_HorizontalKicker(name,fs,fk,element_length);42 temp_kick->setAperture(element_aperture);43 temp_kick->setX(xpos);44 temp_kick->setY(ypos);45 temp_kick->setTX(txpos);46 temp_kick->setTY(typos);47 temp_kick->setBetaX(betax);48 temp_kick->setBetaY(betay);49 return temp_kick;50 }51 -
trunk/external/Hector/H_HorizontalKicker.h
r1360 r1365 1 #ifndef _H_HorizontalKicker_2 #define _H_HorizontalKicker_3 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * *5 * *6 * --<--<-- A fast simulator --<--<-- *7 * / --<--<-- of particle --<--<-- *8 * ----HECTOR----< *9 * \ -->-->-- transport through -->-->-- *10 * -->-->-- generic beamlines -->-->-- *11 * *12 * JINST 2:P09005 (2007) *13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) *14 * http://www.fynu.ucl.ac.be/hector.html *15 * *16 * Center for Cosmology, Particle Physics and Phenomenology *17 * Universite catholique de Louvain *18 * Louvain-la-Neuve, Belgium *19 * *20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */21 22 1 /// \file H_HorizontalKicker.h 23 2 /// \brief Classes aiming at simulating horizontal kickers in beamline. 24 3 /// 25 4 /// fk [rad] for kickers !!!! 5 6 /* 7 ---- Hector the simulator ---- 8 A fast simulator of particles through generic beamlines. 9 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 10 11 http://www.fynu.ucl.ac.be/hector.html 12 13 Centre de Physique des Particules et de Phénoménologie (CP3) 14 Université Catholique de Louvain (UCL) 15 */ 16 17 #ifndef _H_HorizontalKicker_ 18 #define _H_HorizontalKicker_ 26 19 27 20 // local #includes … … 36 29 H_HorizontalKicker():H_Kicker(HKICKER,0.,0.,0.) {init();} 37 30 H_HorizontalKicker(const double s, const double k, const double l) :H_Kicker(HKICKER,s,k,l){init();} 38 H_HorizontalKicker(const string &nameE, const double s, const double k, const double l) :H_Kicker(nameE,HKICKER,s,k,l){init();}39 ~H_HorizontalKicker() { };31 H_HorizontalKicker(const string nameE, const double s, const double k, const double l) :H_Kicker(nameE,HKICKER,s,k,l){init();} 32 ~H_HorizontalKicker() {return;}; 40 33 //@} 41 H_HorizontalKicker* clone() const;42 34 private: 43 35 virtual void setTypeString() {typestring=HKICKERNAME;}; 44 virtual void setMatrix(const float, const float, const float) ;36 virtual void setMatrix(const float, const float, const float) const ; 45 37 }; 46 38 -
trunk/external/Hector/H_HorizontalQuadrupole.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 18 11 19 12 /// \file H_HorizontalQuadrupole.cc 20 /// \brief Classes aiming at simulating horizontal quadrupoles in beamline13 /// \brief Classes aiming at simulating horizontal kickers in beamline 21 14 22 15 // local #includes 23 16 #include "H_HorizontalQuadrupole.h" 24 17 25 void H_HorizontalQuadrupole::setMatrix(const float eloss, const float p_mass, const float p_charge) {26 if (fk>0) { if(VERBOSE) cout<<" <H_HorizontalQuadrupole> ERROR : k1 should be < 0(" << name << ")!"<<endl; }27 if (fk !=0 ) element_mat = hquadmat(element_length,fk,eloss, p_mass, p_charge);18 void H_HorizontalQuadrupole::setMatrix(const float eloss, const float p_mass, const float p_charge) const { 19 if (fk>0) { if(VERBOSE) cout<<"\t ERROR : k1 should be < 0 for H_HorizontalQuadrupole (" << name << ")!"<<endl; } 20 if (fk !=0 ) *element_mat = hquadmat(element_length,fk,eloss, p_mass, p_charge); 28 21 else { 29 element_mat = driftmat(element_length);30 if(VERBOSE) cout<<" <H_HorizontalQuadrupole>WARNING : k1= 0 ; drift-like quadrupole (" << name << ") !" << endl;22 *element_mat = driftmat(element_length); 23 if(VERBOSE) cout<<"\t WARNING : k1= 0 ; drift-like quadrupole (" << name << ") !" << endl; 31 24 } 32 25 return ; 33 26 } 34 35 H_HorizontalQuadrupole* H_HorizontalQuadrupole::clone() const {36 H_HorizontalQuadrupole* temp_quad = new H_HorizontalQuadrupole(name,fs,fk,element_length);37 temp_quad->setAperture(element_aperture);38 temp_quad->setX(xpos);39 temp_quad->setY(ypos);40 temp_quad->setTX(txpos);41 temp_quad->setTY(typos);42 temp_quad->setBetaX(betax);43 temp_quad->setBetaY(betay);44 return temp_quad;45 } -
trunk/external/Hector/H_HorizontalQuadrupole.h
r1360 r1365 2 2 #define _H_HorizontalQuadrupole_ 3 3 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 5 * * 6 * --<--<-- A fast simulator --<--<-- * 7 * / --<--<-- of particle --<--<-- * 8 * ----HECTOR----< * 9 * \ -->-->-- transport through -->-->-- * 10 * -->-->-- generic beamlines -->-->-- * 11 * * 12 * JINST 2:P09005 (2007) * 13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 14 * http://www.fynu.ucl.ac.be/hector.html * 15 * * 16 * Center for Cosmology, Particle Physics and Phenomenology * 17 * Universite catholique de Louvain * 18 * Louvain-la-Neuve, Belgium * 19 * * 20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 4 /* 5 ---- Hector the simulator ---- 6 A fast simulator of particles through generic beamlines. 7 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 8 9 http://www.fynu.ucl.ac.be/hector.html 10 11 Centre de Physique des Particules et de Phénoménologie (CP3) 12 Université Catholique de Louvain (UCL) 13 */ 21 14 22 15 /// \file H_HorizontalQuadrupole.h … … 34 27 H_HorizontalQuadrupole():H_Quadrupole(HQUADRUPOLE,0.,0.,0.) {init();} 35 28 H_HorizontalQuadrupole(const double s, const double k, const double l) : H_Quadrupole(HQUADRUPOLE,s,k,l){init();} 36 H_HorizontalQuadrupole(const string &nameE, const double s, const double k, const double l) : H_Quadrupole(nameE,HQUADRUPOLE,s,k,l){init();}37 ~H_HorizontalQuadrupole() { };29 H_HorizontalQuadrupole(const string nameE, const double s, const double k, const double l) : H_Quadrupole(nameE,HQUADRUPOLE,s,k,l){init();} 30 ~H_HorizontalQuadrupole() {return;}; 38 31 //@} 39 H_HorizontalQuadrupole* clone() const ;40 32 private: 41 33 virtual void setTypeString() {typestring = HQUADRUPOLENAME;}; 42 virtual void setMatrix(const float, const float, const float) ;34 virtual void setMatrix(const float, const float, const float) const ; 43 35 }; 44 36 -
trunk/external/Hector/H_Kicker.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 18 11 19 12 /// \file H_Kicker.cc … … 34 27 } 35 28 36 std::ostream& operator<< (std::ostream& os, const H_Kicker& el) { 37 os << el.typestring << el.name <<"\t at s = "<< el.fs <<"\t length = "<< el.element_length; 38 os<<"\t k0 = "<< el.fk <<endl; 39 if(el.element_aperture->getType()!=NONE) { 40 os << *(el.element_aperture) << endl; 29 void H_Kicker::printProperties() const { 30 cout << typestring; 31 cout << name; 32 cout<<"\t at s = "<< fs; 33 cout<<"\t length = "<< element_length; 34 cout<<"\t k0 = "<<fk; 35 cout<<endl; 36 if(element_aperture->getType()!=NONE) { 37 cout <<"\t aperture type = " << element_aperture->getTypeString(); 38 element_aperture->printProperties(); 41 39 } 42 40 43 if(el.element_length<0 && VERBOSE) os<<"<H_Kicker> ERROR : Interpenetration of elements !"<<endl; 44 else if(el.element_length==0 && VERBOSE) os<<"<H_Kicker> WARNING : 0-length "<< el.name << " !" << endl; 45 return os; 41 if(element_length<0) { if(VERBOSE) cout<<"\t ERROR : Interpenetration of elements !"<<endl; } 42 if(element_length==0) { if(VERBOSE) cout<<"\t WARNING : 0-length "<< name << " !" << endl; } 46 43 } -
trunk/external/Hector/H_Kicker.h
r1360 r1365 1 #ifndef _H_Kicker_2 #define _H_Kicker_3 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * *5 * *6 * --<--<-- A fast simulator --<--<-- *7 * / --<--<-- of particle --<--<-- *8 * ----HECTOR----< *9 * \ -->-->-- transport through -->-->-- *10 * -->-->-- generic beamlines -->-->-- *11 * *12 * JINST 2:P09005 (2007) *13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) *14 * http://www.fynu.ucl.ac.be/hector.html *15 * *16 * Center for Cosmology, Particle Physics and Phenomenology *17 * Universite catholique de Louvain *18 * Louvain-la-Neuve, Belgium *19 * *20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */21 22 1 /// \file H_Kicker.h 23 2 /// \brief Classes aiming at simulating kickers in LHC beamline. 24 3 /// fk [rad] for kickers !!!! 4 5 /* 6 ---- Hector the simulator ---- 7 A fast simulator of particles through generic beamlines. 8 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 9 10 http://www.fynu.ucl.ac.be/hector.html 11 12 Centre de Physique des Particules et de Phénoménologie (CP3) 13 Université Catholique de Louvain (UCL) 14 */ 15 16 #ifndef _H_Kicker_ 17 #define _H_Kicker_ 25 18 26 19 // local #includes … … 35 28 H_Kicker():H_OpticalElement() {} 36 29 H_Kicker(const int dtype, const double s, const double k, const double l):H_OpticalElement(dtype,s,k,l){} 37 H_Kicker(const string &nameE, const int dtype, const double s, const double k, const double l):H_OpticalElement(nameE,dtype,s,k,l){}38 virtual ~H_Kicker() { };30 H_Kicker(const string nameE, const int dtype, const double s, const double k, const double l):H_OpticalElement(nameE,dtype,s,k,l){} 31 virtual ~H_Kicker() {return;}; 39 32 //@} 40 virtual H_Kicker* clone() const = 0;41 virtual void printProperties() const { cout << *this; return;};33 /// prints the kicker properties 34 virtual void printProperties() const; 42 35 void init(); 43 36 44 37 protected: 45 virtual void setTypeString() = 0; 46 virtual void setMatrix(const float, const float, const float) = 0; 47 friend std::ostream& operator<< (std::ostream& os, const H_Kicker& el); 38 virtual void setTypeString() =0; 39 virtual void setMatrix(const float, const float, const float) const =0; 48 40 }; 49 41 -
trunk/external/Hector/H_Marker.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 18 11 19 12 /// \file H_Marker.cc … … 31 24 } 32 25 33 std::ostream& operator<< (std::ostream& os, const H_Marker& el) { 34 os << el.typestring << el.name << "\t\t at s = " << el.fs; 35 if(el.element_aperture->getType()!=NONE) { 36 os << *(el.element_aperture) << endl; 26 void H_Marker::printProperties() const { 27 cout << typestring << name; 28 cout << "\t\t at s = " << fs << endl; 29 if(element_aperture->getType()!=NONE) { 30 cout <<"\t aperture type = " << element_aperture->getTypeString(); 31 element_aperture->printProperties(); 37 32 } 38 return os;39 33 } 40 34 41 //void H_Marker::setMatrix(const float eloss, const float p_mass, const float p_charge) { 42 void H_Marker::setMatrix(const float , const float , const float ) { 43 element_mat = driftmat(0); 35 void H_Marker::setMatrix(const float eloss, const float p_mass, const float p_charge) const { 36 *element_mat = driftmat(0); 44 37 return ; 45 38 } 46 47 H_Marker* H_Marker::clone() const {48 H_Marker* temp_mkr = new H_Marker(name,fs);49 temp_mkr->setAperture(element_aperture);50 temp_mkr->setX(xpos);51 temp_mkr->setY(ypos);52 temp_mkr->setTX(txpos);53 temp_mkr->setTY(typos);54 temp_mkr->setBetaX(betax);55 temp_mkr->setBetaY(betay);56 return temp_mkr;57 } -
trunk/external/Hector/H_Marker.h
r1360 r1365 2 2 #define _H_Marker_ 3 3 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 5 * * 6 * --<--<-- A fast simulator --<--<-- * 7 * / --<--<-- of particle --<--<-- * 8 * ----HECTOR----< * 9 * \ -->-->-- transport through -->-->-- * 10 * -->-->-- generic beamlines -->-->-- * 11 * * 12 * JINST 2:P09005 (2007) * 13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 14 * http://www.fynu.ucl.ac.be/hector.html * 15 * * 16 * Center for Cosmology, Particle Physics and Phenomenology * 17 * Universite catholique de Louvain * 18 * Louvain-la-Neuve, Belgium * 19 * * 20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 4 /* 5 ---- Hector the simulator ---- 6 A fast simulator of particles through generic beamlines. 7 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 8 9 http://www.fynu.ucl.ac.be/hector.html 10 11 Centre de Physique des Particules et de Phénoménologie (CP3) 12 Université Catholique de Louvain (UCL) 13 */ 21 14 22 15 /// \file H_Marker.h … … 24 17 25 18 // local includes 26 #include "H_ Drift.h"19 #include "H_OpticalElement.h" 27 20 28 21 /// \brief Class defining a marker in the beamline (e.g. for interaction point) 29 class H_Marker : public H_ Drift {22 class H_Marker : public H_OpticalElement { 30 23 31 24 public: 32 25 /// Constructors and destructor 33 26 //@{ 34 H_Marker():H_ Drift() { type = MARKER;init();}35 H_Marker(const double s):H_ Drift(s,0.) { type =MARKER;init();}36 H_Marker(const string & nameE, const double s):H_Drift(nameE,s,0.) { type=MARKER;init();}37 ~H_Marker() { };27 H_Marker():H_OpticalElement(MARKER,0.,0.,0.) {init();} 28 H_Marker(const double s):H_OpticalElement(MARKER,s,0.,0.){init();} 29 H_Marker(const string nameE, const double s):H_OpticalElement(nameE,MARKER,s,0.,0.){init();} 30 ~H_Marker() { return; }; 38 31 //@} 39 H_Marker* clone() const;32 virtual void printProperties() const; 40 33 void init(); 41 34 42 35 private: 43 36 virtual void setTypeString() {typestring = MARKERNAME;}; 44 virtual void setMatrix(const float , const float, const float) ; 45 friend std::ostream& operator<< (std::ostream& os, const H_Marker& el); 37 virtual void setMatrix(const float , const float, const float) const ; 46 38 }; 47 39 -
trunk/external/Hector/H_OpticalElement.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 18 11 19 12 /// \file H_OpticalElement.cc … … 25 18 26 19 // ROOT #includes 27 #include "TPaveLabel.h"20 //#include "TPaveLabel.h" 28 21 29 22 // local #includes 30 #include "H_Parameters.h"31 23 #include "H_TransportMatrices.h" 32 24 #include "H_Aperture.h" … … 35 27 36 28 /// called by the constructors 37 void H_OpticalElement::init(const string & nameE, const int typeE, const double s, const double k, const double l) {29 void H_OpticalElement::init(const string nameE, const int typeE, const double s, const double k, const double l, H_Aperture* the_app) { 38 30 // this is called by the constructors 39 31 // must be in public section ! … … 48 40 element_length = l; 49 41 type = typeE; 50 element_mat.ResizeTo(MDIM,MDIM); 51 element_mat = driftmat(l); 52 53 // do NOT use setAperture for the initialisation ! there are protections there 42 element_mat = new TMatrix(MDIM,MDIM); 43 setAperture(the_app); 54 44 55 if(element_length<0) { if(VERBOSE) cout<<" <H_OpticalElement>ERROR : Interpenetration of elements !"<<endl; }56 if(element_length==0) { if(VERBOSE) cout<<" <H_OpticalElement>WARNING : 0-length element ! (" << name << ") " << " at " << fs << endl; }45 if(element_length<0) { if(VERBOSE) cout<<"\t ERROR : Interpenetration of elements !"<<endl; } 46 if(element_length==0) { if(VERBOSE) cout<<"\t WARNING : 0-length element ! (" << name << ") " << " at " << fs << endl; } 57 47 betax =0; 58 48 betay =0; 59 49 } 60 50 61 H_OpticalElement::H_OpticalElement() : element_aperture(new H_Aperture()) { 62 init("",DRIFT,0.,0.,0.1); 51 H_OpticalElement::H_OpticalElement() { 52 H_Aperture* the_app = new H_Aperture(); 53 init("",DRIFT,0.,0.,0.1,the_app); 63 54 } 64 55 65 H_OpticalElement::H_OpticalElement(const string& nameE, const int typeE, const double s, const double k, const double l) : element_aperture(new H_Aperture()) { 66 init(nameE,typeE,s,k,l); 56 H_OpticalElement::H_OpticalElement(const string nameE, const int typeE, const double s, const double k, const double l) { 57 H_Aperture* the_app = new H_Aperture(); 58 init(nameE,typeE,s,k,l,the_app); 67 59 } 68 60 69 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()) {70 init(nameE,typeE,s,k,l );61 H_OpticalElement::H_OpticalElement(const string nameE, const int typeE, const double s, const double k, const double l, H_Aperture* the_app) { 62 init(nameE,typeE,s,k,l,the_app); 71 63 } 72 64 73 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()){74 init("",typeE,s,k,l );65 H_OpticalElement::H_OpticalElement(const int typeE, const double s, const double k, const double l, H_Aperture* the_app) { 66 init("",typeE,s,k,l,the_app); 75 67 } 76 68 77 H_OpticalElement::H_OpticalElement(const int typeE, const double s, const double k, const double l) : element_aperture(new H_Aperture()) { 78 init("",typeE,s,k,l); 69 H_OpticalElement::H_OpticalElement(const int typeE, const double s, const double k, const double l) { 70 H_Aperture* the_app = new H_Aperture(); 71 init("",typeE,s,k,l,the_app); 79 72 } 80 73 … … 92 85 name = el.name; 93 86 typestring = el.typestring; 94 element_mat.ResizeTo(MDIM,MDIM); 95 element_mat = el.element_mat; 96 element_aperture = el.element_aperture->clone(); 87 element_mat = new TMatrix(*(el.element_mat)); 88 element_aperture = new H_Aperture(*(el.element_aperture)); 97 89 } 98 90 … … 104 96 xpos = el.xpos; 105 97 ypos = el.ypos; 106 txpos = el.txpos;107 typos = el.typos;98 txpos = el.txpos; 99 typos = el.typos; 108 100 betax = el.betax; 109 101 betay = el.betay; … … 111 103 name = el.name; 112 104 typestring = el.typestring; 113 element_mat.ResizeTo(MDIM,MDIM); 114 element_mat = el.element_mat; 115 element_aperture = el.element_aperture->clone(); 105 delete element_mat; 106 delete element_aperture; 107 element_mat = new TMatrix(*(el.element_mat)); 108 element_aperture = new H_Aperture(*(el.element_aperture)); 116 109 return *this; 117 110 } 118 111 119 void H_OpticalElement::setAperture(const H_Aperture* ap) { 120 // do NOT use setAperture in your constructor, as element_aperture is not initialized 121 // this function do not take into account ap if ap=0 122 // do nothing if element_mat = ap 123 if (!ap) {cout << "<H_OpticalElement> Trying to set an empty pointer for the aperture ! Nothing done.\n"; return;} 124 if (element_aperture != ap) { 125 delete element_aperture; 126 element_aperture = ap->clone(); 127 } 112 void H_OpticalElement::setAperture(H_Aperture * ap) { 113 // element_aperture = const_cast<H_Aperture*>(ap); 114 element_aperture = ap; 128 115 return; 129 116 } 130 /*131 void H_OpticalElement::setAperture(H_Aperture* ap) {132 // do NOT use setAperture in your constructor, as element_aperture is not initialized133 // this function do not take into account ap if ap=0134 // do nothing if element_mat = ap135 if (!ap) {cout << "<H_OpticalElement> Trying to set an empty pointer for the aperture ! Nothing done.\n"; return;}136 if (element_aperture != ap) {137 delete element_aperture;138 element_aperture = ap; //->clone();139 }140 return;141 }142 */143 117 144 std::ostream& operator<< (std::ostream& os, const H_OpticalElement& el) { 145 os << el.typestring << el.name << "\t at s = " << el.fs << "\t length = "<< el.element_length; 146 if(el.fk!=0) os <<"\t strength = " << el.fk; 147 if(el.element_aperture->getType()!=NONE) { 148 os << *(el.element_aperture) << endl; 149 } 150 os<<endl; 151 if(el.element_length<0 && VERBOSE) 152 os <<"<H_OpticalElement> ERROR : Interpenetration of elements !"<<endl; 153 else if(el.element_length==0 && VERBOSE) 154 os <<"<H_OpticalElement> WARNING : 0-length "<< el.typestring << " !" << endl; 155 return os; 118 void H_OpticalElement::printProperties() const { 119 cout << typestring; 120 cout << name; 121 cout <<"\t at s = " << fs; 122 cout <<"\t length = "<< element_length; 123 if(fk!=0) cout <<"\t strength = " << fk; 124 if(element_aperture->getType()!=NONE) { 125 cout <<"\t aperture type = " << element_aperture->getTypeString(); 126 element_aperture->printProperties(); 127 } 128 129 cout<<endl; 130 if(element_length<0) { if(VERBOSE) cout<<"\t ERROR : Interpenetration of elements !"<<endl; } 131 if(element_length==0) { if(VERBOSE) cout<<"\t WARNING : 0-length "<< typestring << " !" << endl; } 132 133 return; 156 134 } 157 135 158 136 void H_OpticalElement::showMatrix() const { 159 printMatrix(element_mat);137 printMatrix(element_mat); 160 138 return; 161 139 } 162 140 163 141 void H_OpticalElement::drawAperture() const { 164 element_aperture->draw( 1);142 element_aperture->draw(); 165 143 return; 166 144 } 167 145 168 TMatrix H_OpticalElement::getMatrix() { 169 setMatrix(0,MP,1); 170 return element_mat; 146 TMatrix H_OpticalElement::getMatrix() const { 147 return *element_mat; 171 148 } 172 149 173 TMatrix H_OpticalElement::getMatrix(const float eloss, const float p_mass, const float p_charge) {150 TMatrix H_OpticalElement::getMatrix(const float eloss, const float p_mass, const float p_charge) const { 174 151 setMatrix(eloss,p_mass,p_charge); 175 return element_mat;152 return *element_mat; 176 153 } 177 154 178 155 void H_OpticalElement::draw(const float meight, const float height) const{ 179 /// @param meight is the minimal extend of the graph156 /* /// @param meight is the minimal extend of the graph 180 157 /// @param height is the maximal extend of the graph 181 158 float x1 = getS(); … … 188 165 cur_box->SetFillColor((int)getType()); 189 166 cur_box->Draw(); 167 */ 190 168 } 191 169 192 TVectorD H_OpticalElement::getHitPosition(const TVectorD& init_pos, const double energy_loss, const double mp, const double qp) {193 if(!element_length) {194 // cout<<"O-length element ("<<getName()<<"), should not appear here !"<<endl;195 return init_pos;196 }197 // some declarations198 bool inside = false;199 double vec[MDIM] = {init_pos[INDEX_X]/URAD, tan(init_pos[INDEX_TX]/URAD),200 init_pos[INDEX_Y]/URAD, tan(init_pos[INDEX_TY]/URAD),201 -energy_loss, 1};202 TMatrixD mat_init(1,MDIM,vec);203 TMatrixD mat_min(1,MDIM,vec);204 TMatrixD mat_max(1,MDIM,vec);205 TMatrixD mat_stop(1,MDIM,vec);206 H_OpticalElement* temp_el = clone();207 // initialasing boundaries208 double min_pos = 0;209 double max_pos = element_length/2.;210 double max_old = element_length/2.;211 // number of iterations212 // (idea : fix precision instead of number of iterations + add security)213 // (idea : interpolate between max and min and give the error)214 const int N = 10;215 // starting search loop216 for(int i = 0; i < N; i++) {217 // fixing position to be investigated218 temp_el->setLength(max_pos);219 // initialising the vector at the initial vector + possible shift/tilt220 mat_max[0][0] = mat_init[0][0] - temp_el->getX();221 mat_max[0][1] = mat_init[0][1] - tan(temp_el->getTX());222 mat_max[0][2] = mat_init[0][2] - temp_el->getY();223 mat_max[0][3] = mat_init[0][3] - tan(temp_el->getTY());224 // propagating225 mat_max *= temp_el->getMatrix(energy_loss,mp,qp);226 // compensating for the previously stated shifts/tilts227 mat_max[0][0] = mat_max[0][0] + temp_el->getX();228 mat_max[0][1] = mat_max[0][1] + tan(temp_el->getTX());229 mat_max[0][2] = mat_max[0][2] + temp_el->getY();230 mat_max[0][3] = mat_max[0][3] + tan(temp_el->getTY());231 // fixing new boundaries232 if(temp_el->isInside(mat_max.GetMatrixArray()[0]*URAD,mat_max.GetMatrixArray()[2]*URAD)) {233 max_old = max_pos;234 max_pos = max_pos + (max_pos - min_pos)/2.;235 min_pos = max_old;236 inside = true;237 } else {238 max_pos = min_pos + (max_pos - min_pos)/2.;239 inside = false;240 }241 // end of loop242 }243 // if it passes at the last iteration, choosing the other range²244 if(inside) min_pos = max_old;245 // here the interpolation method : now we are sure that the intercept is between min_pos and max_pos246 // getting vector at min_pos (for first boundary) :247 bool precision_estimate = false;248 if(precision_estimate) {249 temp_el->setLength(min_pos);250 mat_min[0][0] = mat_init[0][0] - temp_el->getX();251 mat_min[0][1] = mat_init[0][1] - tan(temp_el->getTX());252 mat_min[0][2] = mat_init[0][2] - temp_el->getY();253 mat_min[0][3] = mat_init[0][3] - tan(temp_el->getTY());254 mat_min *= temp_el->getMatrix(energy_loss,mp,qp);255 mat_min[0][0] = mat_min[0][0] + temp_el->getX();256 mat_min[0][1] = mat_min[0][1] + tan(temp_el->getTX());257 mat_min[0][2] = mat_min[0][2] + temp_el->getY();258 mat_min[0][3] = mat_min[0][3] + tan(temp_el->getTY());259 mat_min[0][4] = min_pos + init_pos[4];260 // getting vector at max_pos (for second boundary) :261 temp_el->setLength(max_pos);262 mat_max[0][0] = mat_init[0][0] - temp_el->getX();263 mat_max[0][1] = mat_init[0][1] - tan(temp_el->getTX());264 mat_max[0][2] = mat_init[0][2] - temp_el->getY();265 mat_max[0][3] = mat_init[0][3] - tan(temp_el->getTY());266 mat_max *= temp_el->getMatrix(energy_loss,mp,qp);267 mat_max[0][0] = mat_max[0][0] + temp_el->getX();268 mat_max[0][1] = mat_max[0][1] + tan(temp_el->getTX());269 mat_max[0][2] = mat_max[0][2] + temp_el->getY();270 mat_max[0][3] = mat_max[0][3] + tan(temp_el->getTY());271 mat_max[0][4] = max_pos + init_pos[4];272 }273 // getting vector in the middle (for estimate) :274 temp_el->setLength((max_pos+min_pos)/2.);275 mat_stop[0][0] = mat_init[0][0] - temp_el->getX();276 mat_stop[0][1] = mat_init[0][1] - tan(temp_el->getTX());277 mat_stop[0][2] = mat_init[0][2] - temp_el->getY();278 mat_stop[0][3] = mat_init[0][3] - tan(temp_el->getTY());279 mat_stop *= temp_el->getMatrix(energy_loss,mp,qp);280 mat_stop[0][0] = mat_stop[0][0] + temp_el->getX();281 mat_stop[0][1] = mat_stop[0][1] + tan(temp_el->getTX());282 mat_stop[0][2] = mat_stop[0][2] + temp_el->getY();283 mat_stop[0][3] = mat_stop[0][3] + tan(temp_el->getTY());284 mat_stop[0][4] = (max_pos+min_pos)/2. + init_pos[4];285 286 double xys[LENGTH_VEC];287 xys[INDEX_X]= mat_stop[0][0]*URAD;288 xys[INDEX_TX]= atan(mat_stop[0][1])*URAD;289 xys[INDEX_Y]= mat_stop[0][2]*URAD;290 xys[INDEX_TY]= atan(mat_stop[0][3])*URAD;291 xys[INDEX_S]= mat_stop[0][4] ;292 TVectorD temp_vec(LENGTH_VEC,xys);293 294 if(precision_estimate) {295 cout<<"--- Results and precision estimates ---"<<endl;296 cout<<"\t Stopping element : "<<getName()<<endl;297 cout<<"\t hit point s : "<<mat_stop[0][4]<<" m +- "<<(mat_max[0][4]-mat_stop[0][4])*1000.<<" mm"<<endl;298 cout<<"\t hit point x : "<<mat_stop[0][0]*URAD;299 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));300 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;301 cout<<"\t hit point y : "<<mat_stop[0][2]*URAD;302 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));303 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;304 cout<<"\t hit point tx : "<<mat_stop[0][1]*URAD;305 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));306 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;307 cout<<"\t hit point ty : "<<mat_stop[0][3]*URAD;308 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));309 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;310 311 }312 313 delete temp_el;314 return temp_vec;315 } -
trunk/external/Hector/H_OpticalElement.h
r1360 r1365 2 2 #define _H_OpticalElement_ 3 3 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 5 * * 6 * --<--<-- A fast simulator --<--<-- * 7 * / --<--<-- of particle --<--<-- * 8 * ----HECTOR----< * 9 * \ -->-->-- transport through -->-->-- * 10 * -->-->-- generic beamlines -->-->-- * 11 * * 12 * JINST 2:P09005 (2007) * 13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 14 * http://www.fynu.ucl.ac.be/hector.html * 15 * * 16 * Center for Cosmology, Particle Physics and Phenomenology * 17 * Universite catholique de Louvain * 18 * Louvain-la-Neuve, Belgium * 19 * * 20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 4 /* 5 ---- Hector the simulator ---- 6 A fast simulator of particles through generic beamlines. 7 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 8 9 http://www.fynu.ucl.ac.be/hector.html 10 11 Centre de Physique des Particules et de Phénoménologie (CP3) 12 Université Catholique de Louvain (UCL) 13 */ 21 14 22 15 /// \file H_OpticalElement.h … … 33 26 // ROOT #includes 34 27 #include "TMatrix.h" 35 #include "TVectorD.h"36 28 37 29 // local #includes … … 43 35 44 36 // type #defines 45 enum {DRIFT=1, RDIPOLE, SDIPOLE, VQUADRUPOLE, HQUADRUPOLE, VKICKER, HKICKER, RCOLLIMATOR, ECOLLIMATOR, CCOLLIMATOR, RP, IP, MARKER}; 37 #define DRIFT 1 38 #define RDIPOLE 2 39 #define SDIPOLE 3 40 #define VQUADRUPOLE 4 41 #define HQUADRUPOLE 5 42 #define VKICKER 6 43 #define HKICKER 7 44 #define RCOLLIMATOR 8 45 #define ECOLLIMATOR 9 46 #define CCOLLIMATOR 10 47 #define RP 11 48 #define IP 12 49 #define MARKER 13 46 50 47 51 // typestring[30] #defines … … 65 69 public: 66 70 /// init method for constructors 67 void init(const string &, const int , const double , const double , const double);71 void init(const string, const int , const double , const double , const double, H_Aperture*); 68 72 /// Constructors and destructor 69 73 //@{ 70 H_OpticalElement(const string &, const int, const double, const double, const double, H_Aperture*);74 H_OpticalElement(const string, const int, const double, const double, const double, H_Aperture*); 71 75 H_OpticalElement(const int, const double, const double, const double, H_Aperture*); 72 H_OpticalElement(const string &, const int, const double, const double, const double);76 H_OpticalElement(const string, const int, const double, const double, const double); 73 77 H_OpticalElement(const int, const double, const double, const double); 74 78 H_OpticalElement(); 75 79 H_OpticalElement(const H_OpticalElement&); 76 virtual ~H_OpticalElement() {delete element_aperture;};80 virtual ~H_OpticalElement() {delete element_mat; delete element_aperture;}; 77 81 //@} 78 82 /// Prints the element features 79 virtual void printProperties() const { cout << *this; return;};83 virtual void printProperties() const ; 80 84 /// Shows the element transport matrix 81 85 void showMatrix() const ; … … 83 87 void drawAperture() const ; 84 88 /// Sets the aperture of the element 85 void setAperture( const H_Aperture*);89 void setAperture(H_Aperture *); 86 90 /// Ordering operator acting on the s coordinate 87 inline bool operator>(const H_OpticalElement & tocomp) const {return ( fs>tocomp.getS() );};91 inline bool operator>(const H_OpticalElement tocomp) const {if(fs>tocomp.getS()) { return true; } else { return false; }}; 88 92 /// Ordering operator acting on the s coordinate 89 inline bool operator<(const H_OpticalElement & tocomp) const {return ( fs<tocomp.getS() );};93 inline bool operator<(const H_OpticalElement tocomp) const {if(fs<tocomp.getS()) { return true; } else { return false; }}; 90 94 /// Copy operator 91 95 H_OpticalElement& operator=(const H_OpticalElement&); … … 93 97 //@{ 94 98 inline void setS(const double new_s) {fs=new_s;}; 95 inline void setLength(const double new_l) { element_length = new_l;};96 99 inline void setX(const double new_pos) { 97 100 /// @param new_pos in [m] … … 136 139 /// Returns the element transport matrix 137 140 //@{ 138 TMatrix getMatrix() ;139 TMatrix getMatrix(const float, const float, const float) ;141 TMatrix getMatrix() const; 142 TMatrix getMatrix(const float, const float, const float) const; 140 143 //@} 141 144 /// Returns the element aperture 142 H_Aperture * getAperture() const {return element_aperture;};145 H_Aperture * getAperture() const {return element_aperture;}; 143 146 /// Sets the beta functions 144 147 //@{ … … 172 175 inline double getRelY() const {return rely;}; 173 176 //@} 174 virtual H_OpticalElement* clone() const { return new H_OpticalElement();};175 TVectorD getHitPosition(const TVectorD& , const double, const double, const double);176 177 177 178 … … 205 206 virtual void setTypeString() {return;}; 206 207 /// Optical element transport matrix. 207 virtual void setMatrix(const float, const float, const float) { cout<<"dummy setmatrix"<<endl;return;};208 virtual void setMatrix(const float, const float, const float) const {return;}; 208 209 /// Optical element transport matrix. 209 TMatrix element_mat; 210 /// Optical element aperture. 211 H_Aperture* element_aperture; 212 213 friend std::ostream& operator<< (std::ostream& os, const H_OpticalElement& el); 210 TMatrix * element_mat; 211 /// Optical element aperture. 212 H_Aperture * element_aperture; 214 213 }; 215 214 -
trunk/external/Hector/H_Parameters.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 18 11 19 12 ///\file H_Parameters.cc -
trunk/external/Hector/H_Parameters.h
r1360 r1365 2 2 #define _Hector_parameters_ 3 3 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 5 * * 6 * --<--<-- A fast simulator --<--<-- * 7 * / --<--<-- of particle --<--<-- * 8 * ----HECTOR----< * 9 * \ -->-->-- transport through -->-->-- * 10 * -->-->-- generic beamlines -->-->-- * 11 * * 12 * JINST 2:P09005 (2007) * 13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 14 * http://www.fynu.ucl.ac.be/hector.html * 15 * * 16 * Center for Cosmology, Particle Physics and Phenomenology * 17 * Universite catholique de Louvain * 18 * Louvain-la-Neuve, Belgium * 19 * * 20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 4 /* 5 ---- Hector the simulator ---- 6 A fast simulator of particles through generic beamlines. 7 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 8 9 http://www.fynu.ucl.ac.be/hector.html 10 11 Centre de Physique des Particules et de Phénoménologie (CP3) 12 Université Catholique de Louvain (UCL) 13 */ 21 14 22 15 /// \file H_Parameters.h … … 25 18 /// Units : angles [\f$ \mu \f$rad], distances [\f$ \mu \f$m], energies [GeV], c=[1]. 26 19 27 #include <cmath>28 29 20 /* from physics and maths */ 30 21 /// proton mass [GeV] 31 22 const double MP=0.93827; 32 23 /// proton charge [e] 33 const double QP=1.; 24 const double QP=1; 25 /// pi 26 #define PI 3.14159265359 34 27 /// conversion factor for \f$\mu\f$rad <-> rad 35 const double URAD=1000000.; 28 #define URAD 1000000. 29 /// in thx = thx* / 3 <==> thx* = THX thx 30 #define THX 3.86 31 /// in dy = 10 thy* <==> thy* = dy / THY 32 #define THY 51.4 36 33 37 34 /* beam parameters */ 38 35 /// beam nominal energy in GeV 36 //#define BE 7000. 39 37 const double BE=7000.; 40 38 /// beam energy divergence, in GeV 41 const double SBE=0.79; 39 #define SBE 0.79 42 40 /// beam nominal energy in TeV 43 const double BETEV=7.; 41 #define BETEV 7. 44 42 /// beam S @ IP 45 const double PS=0.; 43 #define PS 0. 46 44 /// beam X @ IP 47 const double PX=-500.;45 #define PX -500. 48 46 /// beam Y @ IP 49 const double PY=0.; 47 #define PY 0. 50 48 /// beam longitudinal dispersion 51 const double SS=0.; 49 #define SS 0. 52 50 /// beam lateral width SX @ IP 53 const double SX=16.63; 51 #define SX 16.63 52 // #define SX 0. 54 53 /// beam lateral width SY @ IP 55 const double SY=16.63; 54 #define SY 16.63 55 // #define SY 0. 56 56 /// beam transverse direction angle TX @ IP 57 const double TX=0.; 57 #define TX 0. 58 58 /// beam transverse direction angle TY @ IP 59 const double TY=0.; 59 #define TY 0. 60 60 /// beam angular divergence STX @ IP 61 const double STX=30.23; 61 //#define STX 0. 62 #define STX 30.23 62 63 /// beam angular divergence STY @ IP 63 const double STY=30.23; 64 //#define STY 0. 65 #define STY 30.23 66 /// beam dispersion 67 //#define D 120000. 68 const double D=120000.; 64 69 /// half crossing angle at IP [\f$ \mu \f$RAD] 65 const double CRANG=142.5; 70 #define CRANG 142.5 66 71 67 // local defines, used in H_BeamParticle & H_OpticalElements 68 enum {INDEX_X=0, INDEX_TX, INDEX_Y, INDEX_TY, INDEX_S, LENGTH_VEC}; 69 // (x,theta_x,y,theta_y,s) 70 71 /// include Pythia libraries ? (not included on some ROOT installations) 72 //#define _include_pythia_ 73 74 const unsigned int TM = 0; // not used anymore. left for backward compatibility 75 const unsigned int AM = 1; // not used anymore. left for backward compatibility 76 72 /* roman pots parameters */ 73 /// granularity in X position 74 #define GRANPOSX 5. 75 /// granularity in Y position 76 #define GRANPOSY 5. 77 /// granularity in X angle 78 #define GRANANGX 10. 79 /// granularity in Y angle 80 #define GRANANGY 10. 81 /// Distance between rp's 82 #define DISTRP 4000000. 83 /// RP resolution in X, for smearing 84 #define RESX 10. 85 /// RP resolution in Y, for smearing 86 #define RESY 10. 87 /// Radius of the hole in the RP 88 #define RADIUS 1000. 77 89 78 90 /* display parameters */ 79 91 /// Verbose mode ? 80 const bool VERBOSE=false; 92 #define VERBOSE 0 93 94 /// include Pythia libraries ? (not included on some ROOT installations) 95 //#define _include_pythia_ 81 96 82 97 #endif -
trunk/external/Hector/H_Quadrupole.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 18 11 19 12 /// \file H_Quadrupole.cc … … 25 18 // needed for call from H- and V-Quadrupoles constructor 26 19 // must be in public section 27 element_mat.ResizeTo(MDIM,MDIM);28 20 setTypeString(); 29 21 setMatrix(0,MP,QP); … … 31 23 } 32 24 33 std::ostream& operator<< (std::ostream& os, const H_Quadrupole& el) { 34 os << el.typestring << el.name << "\t at s = " << el.fs << "\t length = " << el.element_length << "\t k1 = " << el.fk <<endl; 35 if(el.element_aperture->getType()!=NONE) { 36 os << *(el.element_aperture) << endl; 25 void H_Quadrupole::printProperties() const{ 26 cout << typestring; 27 cout << name; 28 cout << "\t at s = " << fs; 29 cout << "\t length = " << element_length; 30 cout << "\t k1 = " << fk; 31 cout<<endl; 32 if(element_aperture->getType()!=NONE) { 33 cout <<"\t aperture type = " << element_aperture->getTypeString(); 34 element_aperture->printProperties(); 37 35 } 38 36 39 if(el.element_length<0 && VERBOSE) os<<"<H_Quadrupole> ERROR : Interpenetration of elements !"<<endl; 40 else if(el.element_length==0 && VERBOSE) os<<"<H_Quadrupole> WARNING : 0-length "<< el.typestring << " !" << endl; 41 return os; 37 if(element_length<0) { if(VERBOSE) cout<<"\t ERROR : Interpenetration of elements !"<<endl; } 38 if(element_length==0) { if(VERBOSE) cout<<"\t WARNING : 0-length "<< typestring << " !" << endl; } 42 39 } -
trunk/external/Hector/H_Quadrupole.h
r1360 r1365 2 2 #define _H_Quadrupole_ 3 3 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 5 * * 6 * --<--<-- A fast simulator --<--<-- * 7 * / --<--<-- of particle --<--<-- * 8 * ----HECTOR----< * 9 * \ -->-->-- transport through -->-->-- * 10 * -->-->-- generic beamlines -->-->-- * 11 * * 12 * JINST 2:P09005 (2007) * 13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 14 * http://www.fynu.ucl.ac.be/hector.html * 15 * * 16 * Center for Cosmology, Particle Physics and Phenomenology * 17 * Universite catholique de Louvain * 18 * Louvain-la-Neuve, Belgium * 19 * * 20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 4 /* 5 ---- Hector the simulator ---- 6 A fast simulator of particles through generic beamlines. 7 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 8 9 http://www.fynu.ucl.ac.be/hector.html 10 11 Centre de Physique des Particules et de Phénoménologie (CP3) 12 Université Catholique de Louvain (UCL) 13 */ 21 14 22 15 /// \file H_Quadrupole.h … … 34 27 H_Quadrupole():H_OpticalElement() {} 35 28 H_Quadrupole(const int dtype, const double s, const double k, const double l) : H_OpticalElement(dtype,s,k,l) {} 36 H_Quadrupole(const string &nameE, const int dtype, const double s, const double k, const double l) : H_OpticalElement(nameE,dtype,s,k,l) {}37 virtual ~H_Quadrupole() { };29 H_Quadrupole(const string nameE, const int dtype, const double s, const double k, const double l) : H_OpticalElement(nameE,dtype,s,k,l) {} 30 virtual ~H_Quadrupole() {return;}; 38 31 //@} 39 virtual H_Quadrupole* clone() const =0 ; 40 virtual void printProperties() const { cout << *this; return;}; 32 virtual void printProperties() const; 41 33 void init(); 42 34 43 35 protected: 44 36 virtual void setTypeString() =0; 45 virtual void setMatrix(const float, const float, const float) =0;37 virtual void setMatrix(const float, const float, const float) const =0; 46 38 47 friend std::ostream& operator<< (std::ostream& os, const H_Quadrupole& el);48 39 }; 49 40 -
trunk/external/Hector/H_RecRPObject.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 18 11 19 12 /// \file H_RecRPObject.cc 20 /// \brief Class performing the reconstruction based on forward detector measurements 21 /// 22 /// Units : angles [rad], distances [m], energies [GeV], c=[1]. 13 /// \brief Classes aiming at reconstruction particle properties 14 15 // C++ #includes 16 #include <iostream> 17 #include <iomanip> 18 19 // ROOT #includes 20 //#include "TGraph.h" 21 //#include "TF1.h" 22 //#include "TCanvas.h" 23 23 24 24 // local #includes … … 28 28 using namespace std; 29 29 30 // reconstruction class for forward detectors. 31 // Featuring the brand-new reco method from the 32 // louvain group ! 33 34 #define MEGA 1000000. 35 36 H_RecRPObject::H_RecRPObject(): emin(0), emax(-1), x1(0), x2(0), y1(0), y2(0), s1(0), s2(0), 37 txip(NOT_YET_COMPUTED), tyip(NOT_YET_COMPUTED), energy(NOT_YET_COMPUTED), q2(NOT_YET_COMPUTED), pt(NOT_YET_COMPUTED), 38 thebeam(new H_AbstractBeamLine()), 39 f_1(new TF1("f_1","[0] + [1]*x + [2]*x*x ",emin,emax)), 40 f_2(new TF1("f_2","[0] + [1]*x + [2]*x*x ",emin,emax)), 41 g_1(new TF1("g_1","[0] + [1]*x + [2]*x*x ",emin,emax)), 42 g_2(new TF1("g_2","[0] + [1]*x + [2]*x*x ",emin,emax)), 43 d_1(new TF1("d_1","[0] + [1]*x + [2]*x*x ",emin,emax)), 44 d_2(new TF1("d_2","[0] + [1]*x + [2]*x*x ",emin,emax)), 45 k_1(new TF1("k_1","[0] + [1]*x + [2]*x*x ",emin,emax)), 46 k_2(new TF1("k_2","[0] + [1]*x + [2]*x*x ",emin,emax)), 47 l_1(new TF1("l_1","[0] + [1]*x + [2]*x*x ",emin,emax)), 48 l_2(new TF1("l_2","[0] + [1]*x + [2]*x*x ",emin,emax)) 49 {} 50 51 H_RecRPObject::H_RecRPObject(const float ss1, const float ss2, const H_AbstractBeamLine* beam) : emin(0), emax(-1), x1(0), x2(0), y1(0), y2(0), s1(ss1), s2(ss2), 52 txip(NOT_YET_COMPUTED), tyip(NOT_YET_COMPUTED), energy(NOT_YET_COMPUTED), q2(NOT_YET_COMPUTED), pt(NOT_YET_COMPUTED), 53 thebeam(beam->clone()), 54 f_1(new TF1("f_1","[0] + [1]*x + [2]*x*x ",emin,emax)), 55 f_2(new TF1("f_2","[0] + [1]*x + [2]*x*x ",emin,emax)), 56 g_1(new TF1("g_1","[0] + [1]*x + [2]*x*x ",emin,emax)), 57 g_2(new TF1("g_2","[0] + [1]*x + [2]*x*x ",emin,emax)), 58 d_1(new TF1("d_1","[0] + [1]*x + [2]*x*x ",emin,emax)), 59 d_2(new TF1("d_2","[0] + [1]*x + [2]*x*x ",emin,emax)), 60 k_1(new TF1("k_1","[0] + [1]*x + [2]*x*x ",emin,emax)), 61 k_2(new TF1("k_2","[0] + [1]*x + [2]*x*x ",emin,emax)), 62 l_1(new TF1("l_1","[0] + [1]*x + [2]*x*x ",emin,emax)), 63 l_2(new TF1("l_2","[0] + [1]*x + [2]*x*x ",emin,emax)) 64 {if(ss1==ss2) cout<<"<H_RecRPObject> WARNING : detectors are on same position"<<endl; 65 } 66 67 68 H_RecRPObject::H_RecRPObject(const float ss1, const float ss2, const H_AbstractBeamLine& beam) : emin(0), emax(-1), x1(0), x2(0), y1(0), y2(0), s1(ss1), s2(ss2), 69 txip(NOT_YET_COMPUTED), tyip(NOT_YET_COMPUTED), energy(NOT_YET_COMPUTED), q2(NOT_YET_COMPUTED), pt(NOT_YET_COMPUTED), 70 thebeam(beam.clone()), 71 f_1(new TF1("f_1","[0] + [1]*x + [2]*x*x ",emin,emax)), 72 f_2(new TF1("f_2","[0] + [1]*x + [2]*x*x ",emin,emax)), 73 g_1(new TF1("g_1","[0] + [1]*x + [2]*x*x ",emin,emax)), 74 g_2(new TF1("g_2","[0] + [1]*x + [2]*x*x ",emin,emax)), 75 d_1(new TF1("d_1","[0] + [1]*x + [2]*x*x ",emin,emax)), 76 d_2(new TF1("d_2","[0] + [1]*x + [2]*x*x ",emin,emax)), 77 k_1(new TF1("k_1","[0] + [1]*x + [2]*x*x ",emin,emax)), 78 k_2(new TF1("k_2","[0] + [1]*x + [2]*x*x ",emin,emax)), 79 l_1(new TF1("l_1","[0] + [1]*x + [2]*x*x ",emin,emax)), 80 l_2(new TF1("l_2","[0] + [1]*x + [2]*x*x ",emin,emax)) 81 {if(ss1==ss2) cout<<"<H_RecRPObject> WARNING : detectors are on same position"<<endl; 82 } 83 84 85 H_RecRPObject::H_RecRPObject(const H_RecRPObject& r): 86 emin(r.emin), emax(r.emax), x1(r.x1), x2(r.x2), y1(r.y1), y2(r.y2), s1(r.s1), s2(r.s2), 87 txip(r.txip), tyip(r.tyip), energy(r.energy), q2(r.q2), pt(r.pt), 88 //thebeam(r.thebeam->clone()), 89 thebeam(new H_AbstractBeamLine(*(r.thebeam))), 90 f_1(new TF1(*(r.f_1))), f_2(new TF1(*(r.f_2))), g_1(new TF1(*(r.g_1))), g_2(new TF1(*(r.g_2))), 91 d_1(new TF1(*(r.d_1))), d_2(new TF1(*(r.d_2))), k_1(new TF1(*(r.k_1))), k_2(new TF1(*(r.k_2))), 92 l_1(new TF1(*(r.l_1))), l_2(new TF1(*(r.l_2))) 93 {} 30 H_RecRPObject::H_RecRPObject() { 31 x1 = 0.; 32 x2 = 0.; 33 y1 = 0.; 34 y2 = 0.; 35 s1 = 0.; 36 s2 = 0.; 37 corr1_TM = 0; 38 corr2_TM = 0; 39 corr1_AM = 0; 40 corr2_AM = 0; 41 thx = NOT_YET_COMPUTED; 42 thy = NOT_YET_COMPUTED; 43 x0 = NOT_YET_COMPUTED; 44 y0 = NOT_YET_COMPUTED; 45 energy = NOT_YET_COMPUTED; 46 virtuality = NOT_YET_COMPUTED; 47 matrp1 = new TMatrix(MDIM,MDIM); 48 matrp2 = new TMatrix(MDIM,MDIM); 49 thebeam = new H_AbstractBeamLine(); 50 } 51 52 H_RecRPObject::H_RecRPObject(const float S1, const float S2, const H_AbstractBeamLine& beamline) { 53 x1 = 0; 54 x2 = 0; 55 y1 = 0; 56 y2 = 0; 57 s1 = S1; 58 s2 = S2; 59 thx = NOT_YET_COMPUTED; 60 thy = NOT_YET_COMPUTED; 61 x0 = NOT_YET_COMPUTED; 62 y0 = NOT_YET_COMPUTED; 63 energy = NOT_YET_COMPUTED; 64 virtuality = NOT_YET_COMPUTED; 65 // matrp1 = new TMatrix(MDIM,MDIM); 66 // matrp2 = new TMatrix(MDIM,MDIM); 67 thebeam = new H_AbstractBeamLine(beamline); 68 H_RomanPot * rp1 = new H_RomanPot("rp1",s1,0); 69 thebeam->add(rp1); 70 H_RomanPot * rp2 = new H_RomanPot("rp2",s2,0); 71 thebeam->add(rp2); 72 matrp1 = new TMatrix(*(thebeam->getPartialMatrix("rp1",0.,MP,QP))); 73 matrp2 = new TMatrix(*(thebeam->getPartialMatrix("rp2",0.,MP,QP))); 74 75 corr1_TM = getECorrectionFactor(0,TM); 76 corr2_TM = getECorrectionFactor(1,TM); 77 corr1_AM = getECorrectionFactor(0,AM); 78 corr2_AM = getECorrectionFactor(1,AM); 79 // cout << corr1_TM << " " << corr2_TM << endl; 80 // cout << corr1_AM << " " << corr2_AM << endl; 81 } 82 83 H_RecRPObject::H_RecRPObject(const H_RecRPObject& r) { 84 x1 = r.x1; 85 x2 = r.x2; 86 y1 = r.y1; 87 y2 = r.y2; 88 s1 = r.s1; 89 s2 = r.s2; 90 x0 = r.x0; 91 y0 = r.y0; 92 thx = r.thx; 93 thy = r.thy; 94 energy = r.energy; 95 virtuality = r.virtuality; 96 matrp1 = new TMatrix(*(r.matrp1)); 97 matrp2 = new TMatrix(*(r.matrp2)); 98 corr1_TM = r.corr1_TM; 99 corr2_TM = r.corr2_TM; 100 corr1_AM = r.corr1_AM; 101 corr2_AM = r.corr2_AM; 102 thebeam = new H_AbstractBeamLine(*(r.thebeam)); 103 } 94 104 95 105 H_RecRPObject& H_RecRPObject::operator=(const H_RecRPObject& r) { 96 if (this == &r) return *this; 97 emin = r.emin, emax = r.emax; 98 x1 = r.x1; x2 = r.x2; 99 y1 = r.y1; y2 = r.y2; 100 s1 = r.s1; s2 = r.s2; 101 txip= r.txip; tyip=r.tyip; 102 energy= r.energy; q2= r.q2; pt= r.pt; 103 //thebeam = r.thebeam->clone(); 104 thebeam = new H_AbstractBeamLine(*(r.thebeam)); 105 f_1 = new TF1(*(r.f_1)); 106 f_2 = new TF1(*(r.f_2)); 107 g_1 = new TF1(*(r.g_1)); 108 g_2 = new TF1(*(r.g_2)); 109 d_1 = new TF1(*(r.d_1)); 110 d_2 = new TF1(*(r.d_2)); 111 k_1 = new TF1(*(r.k_1)); 112 k_2 = new TF1(*(r.k_2)); 113 l_1 = new TF1(*(r.l_1)); 114 l_2 = new TF1(*(r.l_2)); 115 return *this; 116 } 117 118 void H_RecRPObject::initialize() { 119 // this method sets the functions that will be used for reco later 120 // it should be used only once per beamline after the energy range was fixed. 121 // copying beamline and adding detectors 106 if(this==&r) return *this; 107 x1 = r.x1; 108 x2 = r.x2; 109 y1 = r.y1; 110 y2 = r.y2; 111 s1 = r.s1; 112 s2 = r.s2; 113 x0 = r.x0; 114 y0 = r.y0; 115 thx = r.thx; 116 thy = r.thy; 117 energy = r.energy; 118 virtuality = r.virtuality; 119 matrp1 = new TMatrix(*(r.matrp1)); 120 matrp2 = new TMatrix(*(r.matrp2)); 121 corr1_TM = r.corr1_TM; 122 corr2_TM = r.corr2_TM; 123 corr1_AM = r.corr1_AM; 124 corr2_AM = r.corr2_AM; 125 thebeam = new H_AbstractBeamLine(*(r.thebeam)); 126 return *this; 127 } 128 129 float H_RecRPObject::getECorrectionFactor(const unsigned int facn, const unsigned int method) { 130 /* 131 * commented out because CMSSW does not want any TGraph/TCanvas/TFit 132 * 133 * to be fixed ! 134 * X.R. 07/05/2009 135 * 136 float beta1 = ((thebeam->getPartialMatrix("rp1",0,MP,QP))->GetMatrixArray())[1*MDIM]; 137 float beta2 = ((thebeam->getPartialMatrix("rp2",0,MP,QP))->GetMatrixArray())[1*MDIM]; 138 float disp1 = ((thebeam->getPartialMatrix("rp1",0,MP,QP))->GetMatrixArray())[4*MDIM]*URAD; 139 float disp2 = ((thebeam->getPartialMatrix("rp2",0,MP,QP))->GetMatrixArray())[4*MDIM]*URAD; 140 const int n = 20; //using 20 points to get a good quadratic fit 141 float ee[n], rece[n]; 142 143 for(int i = 0; i < n; i++) { 144 ee[i] = 10 + i*200./((float)n-1); 145 H_BeamParticle p1; 146 p1.emitGamma(ee[i],0.); 147 p1.computePath(thebeam,1); 148 p1.propagate(s1); 149 float x1 = p1.getX(); 150 p1.propagate(s2); 151 float x2 = p1.getX(); 152 switch (method) { 153 case TM: { rece[i] = -x1/disp1; }; break; 154 case AM: { rece[i] = -(beta2*x1-beta1*x2)/(beta2*disp1-beta1*disp2); }; break; 155 case PM: { rece[i] = -x1/disp1; cout<<"this method has not been implemented, using trivial reconstruction"<<endl; } break; 156 default: { rece[i] = -(beta2*x1-beta1*x2)/(beta2*disp1-beta1*disp2); }; break; 157 } 158 ee[i] = ee[i] - rece[i]; 159 } 160 char mytitle[50]; 161 sprintf(mytitle,"c_%d",method); 162 // TCanvas*c = new TCanvas(); 163 // c->SetTitle(mytitle); 164 TGraph* g1 = new TGraph(n,rece,ee); 165 TF1* fit1 = new TF1("fit1","[0]*x + [1]*x*x",0,100); 166 g1->Fit("fit1","Q"); 167 float xfact = fit1->GetParameter(facn); 168 // g1->Draw("APL"); 169 delete g1; 170 delete fit1; 171 172 return xfact; 173 */ 174 return 1.; 175 } 176 177 178 void H_RecRPObject::setPositions(const float X1, const float Y1, const float X2, const float Y2) { 179 thx = NOT_YET_COMPUTED; 180 thy = NOT_YET_COMPUTED; 181 x0 = NOT_YET_COMPUTED; 182 y0 = NOT_YET_COMPUTED; 183 energy = NOT_YET_COMPUTED; 184 virtuality = NOT_YET_COMPUTED; 185 x1 = X1; 186 x2 = X2; 187 y1 = Y1; 188 y2 = Y2; 122 189 123 if(emax<0) { 124 cout<<"<H_RecRPObject> ERROR : energy range has to be set first !"<<endl; 125 cout<<"<H_RecRPObject> Please run setERange() or computeERange()"<<endl; 126 cout<<"<H_RecRPObject> initialization aborted"<<endl; 127 return; 128 } 129 130 if(emax<emin) { 131 cout<<"<H_RecRPObject> ERROR : maximum energy lower than minimum !"<<endl; 132 cout<<"<H_RecRPObject> Please (re-)do setERange()"<<endl; 133 cout<<"<H_RecRPObject> initialization aborted"<<endl; 134 return; 135 } 136 137 H_AbstractBeamLine * b1 = thebeam->clone(); 138 H_RomanPot * rp1 = new H_RomanPot("rp1",s1,0); 139 H_RomanPot * rp2 = new H_RomanPot("rp2",s2,0); 140 b1->add(rp1); 141 b1->add(rp2); 142 143 // fitting parameters 144 const int N = 20; 145 double e_i[N]; 146 double f_1i[N], f_2i[N], g_1i[N], g_2i[N], d_1i[N], d_2i[N]; 147 double k_1i[N], k_2i[N], l_1i[N], l_2i[N]; 190 return; 191 } 192 193 void H_RecRPObject::printProperties() const { 194 cout << "Roman pot variables :" << endl; 195 cout << "\t pot 1 : (x,y,s) = (" << x1 << " , " << y1 << " , " << s1 << " )" << endl; 196 cout << "\t pot 2 : (x,y,s) = (" << x2 << " , " << y2 << " , " << s2 << " )" << endl; 197 cout << endl << "Reconstructed variables :" << endl; 198 cout << "\t IP : (x,y) = (" << x0 << " , " << y0 << ") and (theta_x, theta_y) = (" << thx << " , " << thy << " )" << endl; 199 if (energy==NOT_YET_COMPUTED) cout << "\t Energy not yet computed" << endl; 200 else cout << "\t Energy = " << energy << " GeV" << endl; 201 if (virtuality==NOT_YET_COMPUTED) cout << "\t Virtuality not yet computed" << endl; 202 else cout << "\t Virtuality = " << virtuality << " GeV^2" << endl; 203 cout << endl; 204 } 205 206 float H_RecRPObject::getE() { 207 if(energy==NOT_YET_COMPUTED) { 208 cout<<"Please first compute energy using your favourite method"<<endl; 209 return NOT_YET_COMPUTED; 210 } 211 return energy; 212 } 213 214 float H_RecRPObject::getE(const unsigned int method) { 215 switch (method) { 216 case TM: {energy = computeE_TM();} break; 217 case AM: {energy = computeE_AM();} break; 218 case PM: {energy = computeE_PM();} break; 219 default: {energy = computeE_AM();} break; 220 }; 221 return energy; 222 } 223 224 float H_RecRPObject::computeX0() { 225 if(energy==NOT_YET_COMPUTED) { 226 cout<<"Please first compute energy using your favourite method"<<endl; 227 return NOT_YET_COMPUTED; 228 } 229 float alpha1 = (matrp1->GetMatrixArray())[0]; 230 float alpha2 = (matrp2->GetMatrixArray())[0]; 231 float disp1 = (matrp1->GetMatrixArray())[4*MDIM]*URAD; 232 float disp2 = (matrp2->GetMatrixArray())[4*MDIM]*URAD; 233 x0 = (disp2*x1-disp1*x2)/(disp2*alpha1-disp1*alpha2); 234 return x0; 235 } 236 237 float H_RecRPObject::computeY0() { 238 if(energy==NOT_YET_COMPUTED) { 239 cout<<"Please first compute energy using your favourite method"<<endl; 240 return NOT_YET_COMPUTED; 241 } 242 float gamma1 = (matrp1->GetMatrixArray())[2*MDIM+2]; 243 float gamma2 = (matrp2->GetMatrixArray())[2*MDIM+2]; 244 float delta1 = (matrp1->GetMatrixArray())[3*MDIM+2]; 245 float delta2 = (matrp2->GetMatrixArray())[3*MDIM+2]; 246 y0 = (delta2*y1-delta1*y2)/(delta2*gamma1-delta1*gamma2); 247 return y0; 248 } 249 250 float H_RecRPObject::computeTX() { 251 if(energy==NOT_YET_COMPUTED) { 252 cout<<"Please first compute energy using your favourite method"<<endl; 253 return NOT_YET_COMPUTED; 254 } 255 float beta1 = (matrp1->GetMatrixArray())[1*MDIM]; 256 float beta2 = (matrp2->GetMatrixArray())[1*MDIM]; 257 float disp1 = (matrp1->GetMatrixArray())[4*MDIM]*URAD; 258 float disp2 = (matrp2->GetMatrixArray())[4*MDIM]*URAD; 259 // computes thx (murad) 260 thx = (x1*disp2-x2*disp1)/(beta1*disp2-beta2*disp1); 261 return thx; 262 } 263 264 float H_RecRPObject::computeTY() { 265 if(energy==NOT_YET_COMPUTED) { 266 cout<<"Please first compute energy using your favourite method"<<endl; 267 return NOT_YET_COMPUTED; 268 } 269 float gamma1 = (matrp1->GetMatrixArray())[2*MDIM+2]; 270 float gamma2 = (matrp2->GetMatrixArray())[2*MDIM+2]; 271 float delta1 = (matrp1->GetMatrixArray())[3*MDIM+2]; 272 float delta2 = (matrp2->GetMatrixArray())[3*MDIM+2]; 273 // computes thy (murad) 274 thy = (y1*gamma2-y2*gamma1)/(delta1*gamma2-delta2*gamma1); 275 return thy; 276 } 277 278 float H_RecRPObject::computeE_TM() { 279 // computes the emitted particle energy, from the trivial method 280 float disp = ((thebeam->getPartialMatrix("rp1",0,MP,QP))->GetMatrixArray())[4*MDIM]*URAD; 281 energy = -x1/disp; 282 // corrects for nonlinear effects 283 energy = (1+corr1_TM)*energy + corr2_TM*energy*energy; 284 // sets the rp matrices at obtained energy 285 delete matrp1; 286 delete matrp2; 287 matrp1 = new TMatrix(*(thebeam->getPartialMatrix("rp1",energy,MP,QP))); 288 matrp2 = new TMatrix(*(thebeam->getPartialMatrix("rp2",energy,MP,QP))); 289 // returns ... 290 return energy; 291 } 292 293 float H_RecRPObject::computeE_AM() { 294 // computes the emitted particle energy, from the angle compensation method, iterative way 295 const int N = 10; 296 delete matrp1; 297 delete matrp2; 298 matrp1 = new TMatrix(*(thebeam->getPartialMatrix("rp1",0,MP,QP))); 299 float disp = (matrp1->GetMatrixArray())[4*MDIM]*URAD; 300 delete matrp1; 301 energy = -x1/disp; 302 matrp1 = new TMatrix(*(thebeam->getPartialMatrix("rp1",energy,MP,QP))); 303 matrp2 = new TMatrix(*(thebeam->getPartialMatrix("rp2",energy,MP,QP))); 304 float beta1 = (matrp1->GetMatrixArray())[1*MDIM]; 305 float beta2 = (matrp2->GetMatrixArray())[1*MDIM]; 306 float disp1 = (matrp1->GetMatrixArray())[4*MDIM]*URAD; 307 float disp2 = (matrp2->GetMatrixArray())[4*MDIM]*URAD; 148 308 for(int i = 0; i < N; i++) { 149 e_i[i] = emin + i * (emax - emin)/((double)N-1); 150 // 151 // the bug seems to be linked to the delete operator of the TMatrixT class in root. 152 // valgrind shows memory problems at that point. 153 // for unknwown reasons, copying the matrix gets around this bug. 154 // valgrind (related) ouptut : 155 // 156 // ==13029== Invalid read of size 4 157 // ==13029== at 0x5E75DB6: H_RecRPObject::initialize() (in /home/jdf/GGamma/Hector/lib/libHector.so) 158 // ==13029== by 0x804A2D8: intelligentreco_rpo(double, double, double, double, std::string, int) (H_IntelligentReco.cpp:314) 159 // ==13029== by 0x804A937: main (H_IntelligentReco.cpp:542) 160 // ==13029== Address 0x64AC740 is 80 bytes inside a block of size 144 free'd 161 // ==13029== at 0x4021D18: operator delete[](void*) (vg_replace_malloc.c:256) 162 // ==13029== by 0x5651F57: TMatrixT<float>::Delete_m(int, float*&) (in /home/jdf/root/5.12/lib/libMatrix.so) 163 // ==13029== by 0x565CC5D: TMatrixT<float>::~TMatrixT() (in /home/jdf/root/5.12/lib/libMatrix.so) 164 // ==13029== by 0x5E75D25: H_RecRPObject::initialize() (in /home/jdf/GGamma/Hector/lib/libHector.so) 165 // 166 // 167 TMatrix el_mattt(b1->getPartialMatrix("rp1",e_i[i],MP,QP)); 168 const float *el_mat1 = el_mattt.GetMatrixArray(); 169 // 170 // conclusion : first line of el_mat1 is completely messed-up if it is taken directly from 171 // the return of getpartialmatrix like this : 172 // const float *el_mat1 = (b1->getPartialMatrix("rp1",e_i[i],MP,QP)).GetMatrixArray() 173 // 174 // long-term solution (apart noticing root staff) : replacing el_mat1[i] by the equivalent el_mattt(j,k) 175 // which is anyway more transparent for the reader. 176 // 177 f_1i[i] = el_mat1[0]; 178 g_1i[i] = el_mat1[1*MDIM]; 179 d_1i[i] = MEGA*el_mat1[4*MDIM]; 180 k_1i[i] = el_mat1[2*MDIM+2]; 181 l_1i[i] = el_mat1[3*MDIM+2]; 182 const float* el_mat2 = (b1->getPartialMatrix("rp2",e_i[i],MP,QP).GetMatrixArray()); 183 f_2i[i] = el_mat2[0]; 184 g_2i[i] = el_mat2[1*MDIM]; 185 d_2i[i] = MEGA*el_mat2[4*MDIM]; 186 k_2i[i] = el_mat2[2*MDIM+2]; 187 l_2i[i] = el_mat2[3*MDIM+2]; 188 } 189 TGraph gf_1(N,e_i,f_1i); 190 TGraph gg_1(N,e_i,g_1i); 191 TGraph gd_1(N,e_i,d_1i); 192 TGraph gf_2(N,e_i,f_2i); 193 TGraph gg_2(N,e_i,g_2i); 194 TGraph gd_2(N,e_i,d_2i); 195 TGraph gk_1(N,e_i,k_1i); 196 TGraph gl_1(N,e_i,l_1i); 197 TGraph gk_2(N,e_i,k_2i); 198 TGraph gl_2(N,e_i,l_2i); 199 200 // functions get their final shape 201 gf_1.Fit("f_1","Q"); 202 gg_1.Fit("g_1","Q"); 203 gd_1.Fit("d_1","Q"); 204 gf_2.Fit("f_2","Q"); 205 gg_2.Fit("g_2","Q"); 206 gd_2.Fit("d_2","Q"); 207 gk_1.Fit("k_1","Q"); 208 gl_1.Fit("l_1","Q"); 209 gk_2.Fit("k_2","Q"); 210 gl_2.Fit("l_2","Q"); 211 212 // cleaning the rest 213 delete b1; 214 215 // the end 216 return; 217 } 218 219 void H_RecRPObject::setDetPos(const float ss1, const float ss2) { 309 energy = -(beta2*x1-beta1*x2)/(beta2*disp1-beta1*disp2); 310 delete matrp1; 311 delete matrp2; 312 matrp1 = new TMatrix(*(thebeam->getPartialMatrix("rp1",energy,MP,QP))); 313 matrp2 = new TMatrix(*(thebeam->getPartialMatrix("rp2",energy,MP,QP))); 314 beta1 = (matrp1->GetMatrixArray())[1*MDIM]; 315 beta2 = (matrp2->GetMatrixArray())[1*MDIM]; 316 disp1 = (matrp1->GetMatrixArray())[4*MDIM]*URAD; 317 disp2 = (matrp2->GetMatrixArray())[4*MDIM]*URAD; 318 } 319 // returns ... 320 return energy; 321 } 322 323 float H_RecRPObject::computeE_PM() { 324 cout<<"Not yet implemented, nothing done"<<endl; 220 325 energy = NOT_YET_COMPUTED; 221 s1 = ss1; 222 s2 = ss2; 223 if(ss1==ss2) cout<<"<H_RecRPObject> WARNING : detectors are on same position"<<endl; 224 return; 225 } 226 227 void H_RecRPObject::setPositions(const float xx1, const float xx2, const float yy1, const float yy2) { 228 energy = NOT_YET_COMPUTED; 229 x1 = xx1; 230 x2 = xx2; 231 y1 = yy1; 232 y2 = yy2; 233 return; 234 } 235 236 void H_RecRPObject::setPosition_det1(const float xx1, const float yy1) { 237 energy = NOT_YET_COMPUTED; 238 x1 = xx1; 239 y1 = yy1; 240 } 241 242 void H_RecRPObject::setPosition_det2(const float xx2, const float yy2) { 243 energy = NOT_YET_COMPUTED; 244 x2 = xx2; 245 y2 = yy2; 246 } 247 248 void H_RecRPObject::setERange(const float eemin, const float eemax) { 249 energy = NOT_YET_COMPUTED; 250 emin = eemin; 251 emax = eemax; 252 f_1->SetRange(emin,emax); 253 f_2->SetRange(emin,emax); 254 g_1->SetRange(emin,emax); 255 g_2->SetRange(emin,emax); 256 d_1->SetRange(emin,emax); 257 d_2->SetRange(emin,emax); 258 k_1->SetRange(emin,emax); 259 k_2->SetRange(emin,emax); 260 l_1->SetRange(emin,emax); 261 l_2->SetRange(emin,emax); 262 return; 263 } 264 265 266 void H_RecRPObject::computeERange() { 267 // optional method to determine the energy range of the FIRST detector 268 // in order to refine the fits and get maximum precision. 269 energy = NOT_YET_COMPUTED; 270 H_AbstractBeamLine * b1 = thebeam->clone(); 271 H_RomanPot * rp1 = new H_RomanPot("rp1",s1,0); 272 b1->add(rp1); 273 float max = 1; 274 // number of energies to check 275 const int N = 1000; 276 for(int i=0; i<N; i++) { 277 H_BeamParticle p; 278 p.setE(BE - (emin + i*(BE-emin)/((float)N))); 279 p.computePath(b1); 280 if(p.stopped(b1)) { 281 if(p.getStoppingElement()->getName()=="rp1") { 282 max = emin + i*(BE-emin)/((float)N); 283 } 284 } 285 } 286 cout<<"<H_RecRPObject> Valid energy losses run from 0 (default) to "<<max+20.<<" GeV"<<endl; 287 setERange(0,max+20.); 288 delete b1; 289 return; 290 } 291 292 void H_RecRPObject::computeAll() { 293 // The big game : 294 // computing E, tx, ty, Q2 and Pt and filling the variables. 295 // 296 // The root TF1 class features nice bugs, which explains the 297 // seemingly-dumb structures happening sometimes here as 298 // workarounds for these bugs. The overall thing works very 299 // well but will be cleaned later anyway. 300 301 if(energy!=NOT_YET_COMPUTED) { 302 cout<<"<H_RecRPObject> already computed variables, skipping ..."<<endl; 303 return; 304 } 305 306 TF1 par0("par0","[0]",emin,emax); 307 par0.SetParameter(0,-x1); 308 TF1 par2("par2","[0]",emin,emax); 309 par2.SetParameter(0,-y1); 310 TF1 par1("par1","[0]",emin,emax); 311 par1.SetParameter(0,-x2); 312 TF1 par3("par3","[0]",emin,emax); 313 par3.SetParameter(0,-y2); 314 315 // angle compensating method : 316 TF1 xx_E("xx_E","(g_2*(par0-d_1*x)-g_1*(par1-d_2*x))/(f_2*g_1-f_1*g_2)",emin,emax); 317 TF1 yy_E("yy_E","(par2*l_2 - par3*l_1) / (k_2*l_1 - k_1*l_2)",emin,emax); 318 TF1 xp_E("xp_E","(f_2*(par0-d_1*x)-f_1*(par1-d_2*x))/(g_2*f_1-g_1*f_2)",emin,emax); 319 TF1 yp_E("yp_E","(par2*k_2-par3*k_1)/(l_2*k_1-l_1*k_2)",emin,emax); 320 // it is possible to refine study using y info, but effect was not tested. 321 // TF1 p_xy_E("p_xy_E","(-xx_E*xx_E-yy_E*yy_E)",emin,emax); 322 TF1 p_xy_E("p_xy_E","(-xx_E*xx_E)",emin,emax); 323 324 energy = p_xy_E.GetMaximumX(emin,emax); 325 txip = xp_E.Eval(energy); 326 tyip = yp_E.Eval(energy); 327 pt = sqrt(BE*(BE-energy)*(txip*txip+tyip*tyip)/(MEGA*MEGA)); 328 329 return; 330 } 331 332 float H_RecRPObject::getE(int ) { 333 // put for backward compatibility 334 if(energy==NOT_YET_COMPUTED) { computeAll(); }; 335 return energy; 336 } // to be removed !!!!! 337 338 float H_RecRPObject::getE() { 339 if(energy==NOT_YET_COMPUTED) { computeAll(); }; 340 return energy; 341 } 342 343 float H_RecRPObject::getTX() { 344 if(energy==NOT_YET_COMPUTED) { computeAll(); }; 345 return txip; 346 } 347 348 float H_RecRPObject::getTY() { 349 if(energy==NOT_YET_COMPUTED) { computeAll(); }; 350 return tyip; 351 } 352 353 float H_RecRPObject::getQ2() { 354 if(energy==NOT_YET_COMPUTED) { computeAll(); }; 355 cout<<"<H_RecRPObject::getQ2> Not implemented yet"<<endl; 356 return 0; 357 } 358 359 float H_RecRPObject::getPt() { 360 if(energy==NOT_YET_COMPUTED) { computeAll(); }; 361 return pt; 362 } 363 364 std::ostream& operator<< (std::ostream& os, const H_RecRPObject& rp) { 365 os << "e_min=" << rp.emin << "\t e_max= " << rp.emax << endl; 366 os << "x1=" << rp.x1 << "\t x2= " << rp.x2 << "\t y1=" << rp.y1 << "\t y2=" << rp.y2 367 << "\t s1=" << rp.s1 << "\t s2=" << rp.s2 << endl; 368 os << "txip=" << rp.txip << "\t tyip=" << rp.tyip << "\t energy=" << rp.energy << "\t q2=" << rp.q2 << "\t pt=" << rp.pt << endl; 369 return os; 370 } 371 326 return energy; 327 } 328 329 float H_RecRPObject::computeQ2() { 330 // computes the emitted particle virtuality 331 // energy should be teconstructed first 332 if(energy==NOT_YET_COMPUTED) { 333 cout<<"Please first compute energy using your favourite method"<<endl; 334 return NOT_YET_COMPUTED; 335 } 336 // getting parameters for reconstructed particle energy 337 float beta1 = (matrp1->GetMatrixArray())[1*MDIM]; 338 float beta2 = (matrp2->GetMatrixArray())[1*MDIM]; 339 float gamma1 = (matrp1->GetMatrixArray())[2*MDIM+2]; 340 float gamma2 = (matrp2->GetMatrixArray())[2*MDIM+2]; 341 float delta1 = (matrp1->GetMatrixArray())[3*MDIM+2]; 342 float delta2 = (matrp2->GetMatrixArray())[3*MDIM+2]; 343 float disp1 = (matrp1->GetMatrixArray())[4*MDIM]*URAD; 344 float disp2 = (matrp2->GetMatrixArray())[4*MDIM]*URAD; 345 // angles reconstruction 346 float rec_thx = (x1*disp2-x2*disp1)/(beta1*disp2-beta2*disp1)/URAD; 347 float rec_thy = (y1*gamma2-y2*gamma1)/(delta1*gamma2-delta2*gamma1)/URAD; 348 // Q² reconstruction 349 virtuality = BE*(BE-energy)*(rec_thx*rec_thx+rec_thy*rec_thy); 350 // returns ... 351 return virtuality; 352 } 353 -
trunk/external/Hector/H_RecRPObject.h
r1360 r1365 2 2 #define _H_RecRPObject_ 3 3 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 5 * * 6 * --<--<-- A fast simulator --<--<-- * 7 * / --<--<-- of particle --<--<-- * 8 * ----HECTOR----< * 9 * \ -->-->-- transport through -->-->-- * 10 * -->-->-- generic beamlines -->-->-- * 11 * * 12 * JINST 2:P09005 (2007) * 13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 14 * http://www.fynu.ucl.ac.be/hector.html * 15 * * 16 * Center for Cosmology, Particle Physics and Phenomenology * 17 * Universite catholique de Louvain * 18 * Louvain-la-Neuve, Belgium * 19 * * 20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 4 /* 5 ---- Hector the simulator ---- 6 A fast simulator of particles through generic beamlines. 7 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 8 9 http://www.fynu.ucl.ac.be/hector.html 10 11 Centre de Physique des Particules et de Phénoménologie (CP3) 12 Université Catholique de Louvain (UCL) 13 */ 21 14 22 15 /// \file H_RecRPObject.h 23 /// \brief Reconstructed information from objects detected by two roman pots 16 /// \brief Reconstructed information from objects detected by two roman pots 24 17 25 #include "TF1.h" 26 #include <cmath> 27 28 #include "H_BeamLine.h" 18 // local #includes 19 #include "H_Parameters.h" 20 #include "H_AbstractBeamLine.h" 21 #include <math.h> 22 using namespace std; 29 23 30 24 #define NOT_YET_COMPUTED -666 25 // trivial (TM), angle compensation (AM) and position compensation (PM) methods 26 #define TM 1 27 #define AM 2 28 #define PM 3 31 29 32 /// Reconstruct ion of parameters at IP using measurements withroman pots30 /// Reconstructed information from objects detected by two roman pots 33 31 class H_RecRPObject { 32 /// Uses (s1,x1,y1) and (s2,x2,y2) to compute the energy, the virtuality of the event, as well as the positions (x,y) and angles (thx, thy) at IP. 34 33 public: 34 /// Constructors, destructor and operators 35 //@{ 35 36 H_RecRPObject(); 36 H_RecRPObject(const float, const float, const H_AbstractBeamLine* ); 37 H_RecRPObject(const float, const float, const H_AbstractBeamLine& ); // old-fashioned -- to be deleted 37 H_RecRPObject(const float, const float, const H_AbstractBeamLine& ); 38 38 H_RecRPObject(const H_RecRPObject&); 39 39 H_RecRPObject& operator=(const H_RecRPObject&); 40 ~H_RecRPObject() {delete f_1; delete f_2; delete g_1; delete g_2; 41 delete d_1; delete d_2; delete k_1; delete k_2; 42 delete l_1; delete l_2; 43 delete thebeam;}; 40 ~H_RecRPObject() {delete matrp1; delete matrp2; delete thebeam; return;}; 41 //@} 44 42 45 inline float getX1() const {return x1;}; 46 inline float getX2() const {return x2;}; 47 inline float getY1() const {return y1;}; 48 inline float getY2() const {return y2;}; 49 inline float getS1() const {return s1;}; 50 inline float getS2() const {return s2;}; 43 /// Getters 44 //@{ 45 inline float getX1() const {return x1; /*horizontal position at first roman pot, in \mu m */} 46 inline float getY1() const {return y1; /*vertical position at first roman pot, in \mu m */} 47 inline float getS1() const {return s1; /*longitudinal position of the first roman pot, in m, from IP*/} 48 inline float getX2() const {return x2; /*horizontal position at second RP in \mu m*/} 49 inline float getY2() const {return y2; /*vertical position at second RP in \mu m*/} 50 inline float getS2() const {return s2; /*longitudinal position of the second RP, in m, from IP*/} 51 inline float getTXRP() const {return URAD*atan((x2-x1)/((s2-s1)*URAD)); /* horizontal angle at first RP, in \mu rad*/ } 52 inline float getTYRP() const {return URAD*atan((y2-y1)/((s2-s1)*URAD)); /* vertical angle at first RP, in \mu rad*/ } 53 float getX0() {return (x0==NOT_YET_COMPUTED) ? computeX0():x0; /*reconstructed horizontal position at IP*/} 54 float getY0() {return (y0==NOT_YET_COMPUTED) ? computeY0():y0; /*reconstructed vertical position at IP*/} 55 float getTXIP() {return (thx==NOT_YET_COMPUTED) ? computeTX():thx; /*reconstructed horizontal angle at IP*/} 56 float getTYIP() {return (thy==NOT_YET_COMPUTED) ? computeTY():thy; /*reconstructed vertical angle at IP*/} 57 float getE(const unsigned int); /*returns the reconstructed energy*/ 58 float getE(); /*returns the reconstructed energy if already computed*/ 59 float getQ2() {return (virtuality==NOT_YET_COMPUTED) ? computeQ2():virtuality; /*returns the reconstructed virtuality*/} 60 //@} 51 61 52 float getTX(); 53 float getTY(); 54 float getE(); 55 float getE(int ); 56 float getQ2(); 57 float getPt(); 62 // Sets the proton hit positions 63 void setPositions(const float, const float, const float, const float); 64 // Shows the variable content. 65 void printProperties() const; 66 protected: 67 /// Measured particle coordinates at RP (X - horizontal and Y - vertical in [\f$ \mu \f$m], S -longitudinal in [m]) 68 //@{ 69 float x1, x2, y1, y2, s1, s2; 70 //@} 58 71 59 void setPositions(const float, const float, const float, const float); 60 void setPosition_det1(const float, const float); 61 void setPosition_det2(const float, const float); 62 void setDetPos(const float, const float); 63 void setERange(const float, const float); 64 void computeERange(); 65 void initialize(); 66 void computeAll(); 67 68 protected: 69 70 float emin, emax; 71 float x1, x2, y1, y2, s1, s2; 72 float txip, tyip, energy, q2, pt; 73 H_AbstractBeamLine* thebeam; 74 TF1* f_1, *f_2; 75 TF1* g_1, *g_2; 76 TF1* d_1, *d_2; 77 TF1* k_1, *k_2; 78 TF1* l_1, *l_2; 79 friend std::ostream& operator<< (std::ostream& os, const H_RecRPObject& rp); 72 /// Reconstructed positions and angles at IP in [\f$ \mu m\f$] and [\f$ \mu rad\f$] 73 //@{ 74 float x0, y0, thx, thy; 75 //@} 76 77 /// Reconstructed energy and virtuality at IP in GeV and GeV\f$^2\f$ 78 //@{ 79 float energy, virtuality; 80 //@} 81 82 float computeX0(); 83 float computeY0(); 84 float computeTX(); 85 float computeTY(); 86 /// Energy reconstruction : trivial method 87 float computeE_TM(); 88 /// Energy reconstruction : angle compensation method 89 float computeE_AM(); 90 /// Energy reconstruction : position compensation method 91 float computeE_PM(); 92 /// Virtuality reconstruction. Energy should be reconstructed before. 93 float computeQ2(); 94 /// Calibrates the energy reconstruction with respect to the chromaticity of the transfer matrices 95 float getECorrectionFactor(const unsigned int, const unsigned int ); 96 /// The beamline : 97 H_AbstractBeamLine * thebeam; 98 /// The matrices 99 //@{ 100 TMatrix * matrp1; 101 TMatrix * matrp2; 102 //@} 103 /// The correction factors 104 //@{ 105 float corr1_TM, corr2_TM, corr1_AM, corr2_AM; 106 //@} 80 107 }; 81 108 82 109 #endif 83 -
trunk/external/Hector/H_RectEllipticAperture.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 18 11 19 12 /// \file H_RectEllipticAperture.cc … … 27 20 28 21 // ROOT #includes 29 #include "TPolyLine.h" 22 //#include "TEllipse.h" 23 //#include "TPave.h" 30 24 31 25 // local #includes … … 33 27 using namespace std; 34 28 35 H_RectEllipticAperture::H_RectEllipticAperture(const float l, const float h, const float L, const float H, const float posx, const float posy) :H_Aperture(RECTELLIPSE,((l==0)?L:l),((h==0)?H:h),L,H,posx,posy) {36 /// @param l, h, L, Hare the geometrical parameters of the rect-ellipse shape29 H_RectEllipticAperture::H_RectEllipticAperture(const float x1, const float x2, const float x3, const float x4, const float posx, const float posy) :H_Aperture(RECTELLIPSE,((x1==0)?x3:x1),((x2==0)?x4:x2),x3,x4,posx,posy) { 30 /// @param x1, x2, x3, x4 are the geometrical parameters of the rect-ellipse shape 37 31 /// @param posx, posy defines the (x,y) of the center of the shape 32 type= RECTELLIPSE; 38 33 } 39 34 40 H_RectEllipticAperture* H_RectEllipticAperture::clone() const { 41 return new H_RectEllipticAperture(x1,x2,x3,x4,fx,fy); 42 } 43 44 45 TPolyLine * rectellipse(const float a_e = 2, const float b_e = 1, const float a_r = 1, const float b_r = 2, const float center_x = 0, const float center_y =0) { 46 const int n = 20; // number of points per segment 47 const int N = 4*n; // there are 4 segments 48 float x[N+1], y[N+1]; 49 50 if(a_e>a_r) { 51 // a rectellipse has 4 segments 52 // 1) upper one 53 for (int i=0; i<n; i++) { 54 x[i] = -a_r + i*(2*a_r)/(float)n; 55 y[i] = b_e * sqrt(1-pow(x[i]/a_e,2)); 56 } 57 58 // 2) right vertical segment 59 // upper right corner 60 const float y2 = b_e * sqrt(1-pow(a_r/a_e,2)); 61 // lower right corner 62 const float y3 = -b_e * sqrt(1-pow(a_r/a_e,2)); 63 for (int i=n; i<2*n; i++) { 64 x[i] = a_r; 65 y[i] = y2 - (i-n)*(2*y2)/(float)n; 66 } 67 68 // 3) lower side 69 for (int i=2*n; i<3*n; i++) { 70 x[i] = a_r - (i-2*n)*(2*a_r)/(float)n; 71 y[i] = -b_e * sqrt(1-pow(x[i]/a_e,2)); 72 } 73 74 // 4) left vertical segment 75 // lower left corner 76 const float y4 = y3; 77 for (int i=3*n; i<4*n; i++) { 78 x[i] = -a_r; 79 y[i] = y4 + (i-3*n)*(2*y2)/(float)n; 80 } 81 } else { 82 // 1) upper one : flat 83 const float x1 = -a_e * sqrt(1-pow(b_r/b_e,2)); 84 const float x2 = -x1; 85 for (int i=0; i<n; i++) { 86 y[i] = b_r; 87 x[i] = x1 + i * (x2-x1)/(float)n; 88 } 89 90 // 2) right curved border 91 for (int i=n; i<2*n; i++) { 92 y[i] = b_r - (i-n) * (2*b_r)/(float)n; 93 x[i] = a_e * sqrt(1-pow(y[i]/b_e,2)); 94 } 95 96 // 3) lower side : flat 97 for (int i=2*n; i<3*n; i++) { 98 y[i] = -b_r; 99 x[i] = x2 - (i-2*n) * (2*x2)/(float)n; 100 } 101 102 // 4) left curved border 103 for (int i=3*n; i<4*n; i++) { 104 y[i] = -b_r + (i-3*n) * (2*b_r)/(float)n; 105 x[i] = -a_e * sqrt(1-pow(y[i]/b_e,2)); 106 } 107 } 108 109 // closing the polyline 110 x[N] = x[0]; 111 y[N] = y[0]; 112 113 // shifting the center 114 for (int i=0; i<N+1; i++) { 115 x[i] += center_x; 116 y[i] += center_y; 117 } 118 119 return new TPolyLine(N+1,x,y); 120 } 121 122 123 void H_RectEllipticAperture::draw(const float scale) const { 124 TPolyLine * re = rectellipse(x3*scale, x4*scale, x1*scale, x2*scale, fx*scale, fy*scale); 125 re->SetLineColor(39); 126 re->SetLineWidth(2); 127 re->Draw("l"); 35 void H_RectEllipticAperture::draw() const { 36 /* TEllipse* te = new TEllipse(fx,fy,x3,x4); 37 TPave* tp = new TPave(fx-x1,fy-x2,fx+x1,fy+x2,1); 38 te->SetLineColor(2); 39 te->SetFillColor(2); 40 te->SetFillStyle(3004); 41 te->Draw(); 42 tp->SetLineColor(2); 43 tp->SetFillColor(2); 44 tp->SetFillStyle(3005); 45 tp->Draw(); 128 46 return; 47 */ 129 48 } 130 49 … … 133 52 } 134 53 135 std::ostream& operator<< (std::ostream& os, const H_RectEllipticAperture& ap){136 os << "Aperture shape:" << ap.aptypestring << ", parameters " << ap.x1 <<", "<< ap.x2 <<", "<< ap.x3 <<", "<< ap.x4<< endl;137 os << " \t Center : " << ap.fx <<", "<< ap.fy<<endl;138 return os;54 void H_RectEllipticAperture::printProperties() const { 55 cout << "Aperture shape:" << getTypeString() << ", parameters " <<x1<<", "<<x2<<", "<<x3<<", "<<x4<< endl; 56 cout << " \t Center : "<<fx<<", "<<fy<<endl; 57 return; 139 58 } -
trunk/external/Hector/H_RectEllipticAperture.h
r1360 r1365 2 2 #define _H_RectEllipticAperture_ 3 3 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * *5 * *6 * --<--<-- A fast simulator --<--<-- *7 * / --<--<-- of particle --<--<-- *8 * ----HECTOR----< *9 * \ -->-->-- transport through -->-->-- *10 * -->-->-- generic beamlines -->-->-- *11 * *12 * JINST 2:P09005 (2007) *13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) *14 * http://www.fynu.ucl.ac.be/hector.html *15 * *16 * Center for Cosmology, Particle Physics and Phenomenology *17 * Universite catholique de Louvain *18 * Louvain-la-Neuve, Belgium *19 * *20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */21 22 4 /// \file H_RectEllipticAperture.h 23 5 /// \brief Defines the Rect-Elliptic aperture of beamline elements. 6 7 /* 8 ---- Hector the simulator ---- 9 A fast simulator of particles through generic beamlines. 10 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 11 12 http://www.fynu.ucl.ac.be/hector.html 13 14 Centre de Physique des Particules et de Phénoménologie (CP3) 15 Université Catholique de Louvain (UCL) 16 */ 24 17 25 18 // local #includes … … 34 27 H_RectEllipticAperture():H_Aperture(RECTELLIPSE,0,0,0,0,0,0) {} 35 28 H_RectEllipticAperture(const float,const float,const float,const float, const float, const float); 36 ~H_RectEllipticAperture() {}; 37 H_RectEllipticAperture* clone() const; 29 ~H_RectEllipticAperture() {return;}; 38 30 //@} 31 virtual void printProperties() const; 39 32 /// Checks whether the point is inside the aperture or not 40 33 virtual bool isInside(const float, const float) const; 41 34 /// Draws the aperture shape. 42 virtual void draw(const float scale=1) const; 43 friend std::ostream& operator<< (std::ostream& os, const H_RectEllipticAperture& ap); 35 virtual void draw() const; 44 36 }; 45 37 -
trunk/external/Hector/H_RectangularAperture.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 18 11 19 12 /// \file H_RectangularAperture.cc … … 27 20 28 21 // ROOT #includes 29 #include "TPave.h"22 //#include "TPave.h" 30 23 31 24 // local #includes … … 33 26 using namespace std; 34 27 35 H_RectangularAperture::H_RectangularAperture(const float l, const float h, const float posx, const float posy) :H_Aperture(RECTANGULAR,l,h,0,0,posx,posy) {36 /// @param l, hare the length and height of the rectangular shape28 H_RectangularAperture::H_RectangularAperture(const float x1, const float x2, const float posx, const float posy) :H_Aperture(RECTANGULAR,x1,x2,0,0,posx,posy) { 29 /// @param x1, x2 are the length and height of the rectangular shape 37 30 /// @param posx, posy are the (x,y) coordinates of the center of the rectangular shape 38 }39 40 H_RectangularAperture* H_RectangularAperture::clone() const {41 return new H_RectangularAperture(x1,x2,fx,fy);42 31 } 43 32 … … 47 36 } 48 37 49 void H_RectangularAperture::draw( const float scale) const {50 TPave* tp = new TPave((fx-x1)*scale,(fy-x2)*scale,(fx+x1)*scale,(fy+x2)*scale,1);38 void H_RectangularAperture::draw() const { 39 /* TPave* tp = new TPave(fx-x1,fy-x2,fx+x1,fy+x2,1); 51 40 tp->SetFillStyle(3003); 52 41 tp->SetLineColor(2); … … 54 43 tp->Draw(); 55 44 return; 45 */ 56 46 } 57 58 std::ostream& operator<< (std::ostream& os, const H_RectangularAperture& ap) { 59 os << "Aperture shape:" << ap.aptypestring << ", rectangle sides : " << ap. x1 <<", " << ap.x2 <<endl; 60 os << " \t Center : " << ap.fx << "," << ap.fy << endl; 61 return os; 47 void H_RectangularAperture::printProperties() const { 48 cout << "Aperture shape:" << getTypeString() << ", rectangle Sides : "<<x1<<", "<<x2<<endl; 49 cout << " \t Center : " << fx << "," << fy << endl; 50 return; 62 51 } 63 -
trunk/external/Hector/H_RectangularAperture.h
r1360 r1365 2 2 #define _H_RectangularAperture_ 3 3 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 5 * * 6 * --<--<-- A fast simulator --<--<-- * 7 * / --<--<-- of particle --<--<-- * 8 * ----HECTOR----< * 9 * \ -->-->-- transport through -->-->-- * 10 * -->-->-- generic beamlines -->-->-- * 11 * * 12 * JINST 2:P09005 (2007) * 13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 14 * http://www.fynu.ucl.ac.be/hector.html * 15 * * 16 * Center for Cosmology, Particle Physics and Phenomenology * 17 * Universite catholique de Louvain * 18 * Louvain-la-Neuve, Belgium * 19 * * 20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 4 /* 5 ---- Hector the simulator ---- 6 A fast simulator of particles through generic beamlines. 7 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 8 9 http://www.fynu.ucl.ac.be/hector.html 10 11 Centre de Physique des Particules et de Phénoménologie (CP3) 12 Université Catholique de Louvain (UCL) 13 */ 21 14 22 15 /// \file H_RectangularAperture.h … … 34 27 H_RectangularAperture():H_Aperture(RECTANGULAR,0,0,0,0,0,0) {} 35 28 H_RectangularAperture(const float,const float,const float,const float); 36 ~H_RectangularAperture() {}; 37 H_RectangularAperture* clone() const; 29 ~H_RectangularAperture() {return;}; 38 30 //@} 39 31 /// Checks whether the point is inside the aperture or not 40 32 virtual bool isInside(const float, const float) const; 41 33 /// Draws the aperture shape. 42 virtual void draw( const float scale=1) const;43 friend std::ostream& operator<< (std::ostream& os, const H_RectangularAperture& ap);34 virtual void draw() const; 35 virtual void printProperties() const; 44 36 }; 45 37 -
trunk/external/Hector/H_RectangularCollimator.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 18 11 19 12 /// \file H_RectangularCollimator.cc … … 31 24 } 32 25 33 H_RectangularCollimator::H_RectangularCollimator(const double s, const double l) :H_Drift(s,l){ 34 type = RCOLLIMATOR; 35 init(); 26 H_RectangularCollimator::H_RectangularCollimator(const double s, const double l) :H_OpticalElement(RCOLLIMATOR,s,0.,l){ 27 init(); 36 28 return; 37 29 } 38 30 39 H_RectangularCollimator::H_RectangularCollimator(const string& nameE, const double s, const double l) :H_Drift(nameE,s,l){ 40 type = RCOLLIMATOR; 41 init(); 31 H_RectangularCollimator::H_RectangularCollimator(const string nameE, const double s, const double l) :H_OpticalElement(nameE,RCOLLIMATOR,s,0.,l){ 32 init(); 42 33 return; 43 34 } 44 35 45 std::ostream& operator<< (std::ostream& os, const H_RectangularCollimator& el) { 46 os << el.typestring << el.name << "\t\t at s = " << el.fs << "\t length = " << el.element_length <<endl; 47 if(el.element_aperture->getType()!=NONE) { 48 os << *(el.element_aperture) << endl; 36 void H_RectangularCollimator::printProperties() const { 37 cout << typestring << name; 38 cout << "\t\t at s = " << fs; 39 cout << "\t length = " << element_length; 40 cout<<endl; 41 if(element_aperture->getType()!=NONE) { 42 cout <<"\t aperture type = " << element_aperture->getTypeString(); 43 element_aperture->printProperties(); 49 44 } 50 45 51 if(el.element_length<0 && VERBOSE) os <<"<H_RectangularCollimator> ERROR : Interpenetration of elements !"<<endl; 52 else if(el.element_length==0 && VERBOSE) os <<"<H_RectangularCollimator> WARNING : 0-length "<< el.name << " !" << endl; 53 return os; 46 if(element_length<0) { if(VERBOSE) cout<<"\t ERROR : Interpenetration of elements !"<<endl; } 47 if(element_length==0) { if(VERBOSE) cout<<"\t WARNING : 0-length "<< name << " !" << endl; } 54 48 } 55 49 56 //void H_RectangularCollimator::setMatrix(const float eloss, const float p_mass, const float p_charge){57 void H_RectangularCollimator::setMatrix(const float , const float , const float ) { 58 element_mat = driftmat(element_length);50 void H_RectangularCollimator::setMatrix(const float eloss, const float p_mass, const float p_charge) const { 51 // cout<<"\t WARNING : H_RectangularCollimator matrices not yet implemented ! " << endl; 52 *element_mat = driftmat(element_length); 59 53 return ; 60 54 } 61 62 H_RectangularCollimator* H_RectangularCollimator::clone() const {63 H_RectangularCollimator* temp_coll = new H_RectangularCollimator(name,fs,element_length);64 temp_coll->setAperture(element_aperture);65 temp_coll->setX(xpos);66 temp_coll->setY(ypos);67 temp_coll->setTX(txpos);68 temp_coll->setTY(typos);69 temp_coll->setBetaX(betax);70 temp_coll->setBetaY(betay);71 return temp_coll;72 }73 -
trunk/external/Hector/H_RectangularCollimator.h
r1360 r1365 2 2 #define _H_RectangularCollimator_ 3 3 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 5 * * 6 * --<--<-- A fast simulator --<--<-- * 7 * / --<--<-- of particle --<--<-- * 8 * ----HECTOR----< * 9 * \ -->-->-- transport through -->-->-- * 10 * -->-->-- generic beamlines -->-->-- * 11 * * 12 * JINST 2:P09005 (2007) * 13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 14 * http://www.fynu.ucl.ac.be/hector.html * 15 * * 16 * Center for Cosmology, Particle Physics and Phenomenology * 17 * Universite catholique de Louvain * 18 * Louvain-la-Neuve, Belgium * 19 * * 20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 4 /* 5 ---- Hector the simulator ---- 6 A fast simulator of particles through generic beamlines. 7 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 8 9 http://www.fynu.ucl.ac.be/hector.html 10 11 Centre de Physique des Particules et de Phénoménologie (CP3) 12 Université Catholique de Louvain (UCL) 13 */ 21 14 22 15 /// \file H_RectangularCollimator.h … … 24 17 25 18 // local #includes 26 #include "H_ Drift.h"19 #include "H_OpticalElement.h" 27 20 28 21 /// R-Collimators for the LHC beamline. 29 class H_RectangularCollimator : public H_ Drift {22 class H_RectangularCollimator : public H_OpticalElement { 30 23 31 24 public: 32 25 /// Constructors and destructor 33 26 //@{ 34 H_RectangularCollimator():H_ Drift() {type = RCOLLIMATOR;init();}27 H_RectangularCollimator():H_OpticalElement(RCOLLIMATOR,0.,0.,0.) {init();} 35 28 H_RectangularCollimator(const double, const double ); 36 H_RectangularCollimator(const string &, const double, const double );37 ~H_RectangularCollimator() { };29 H_RectangularCollimator(const string, const double, const double ); 30 ~H_RectangularCollimator() { return; }; 38 31 //@} 39 H_RectangularCollimator* clone() const;32 virtual void printProperties() const; 40 33 void init(); 34 41 35 private: 42 36 virtual void setTypeString() {typestring=RCOLLIMATORNAME;}; 43 virtual void setMatrix(const float, const float, const float) ; 44 45 friend std::ostream& operator<< (std::ostream& os, const H_RectangularCollimator& el); 37 virtual void setMatrix(const float, const float, const float) const ; 38 46 39 }; 47 40 -
trunk/external/Hector/H_RectangularDipole.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 18 11 19 12 /// \file H_RectangularDipole.cc … … 23 16 #include "H_TransportMatrices.h" 24 17 25 void H_RectangularDipole::setMatrix(const float eloss, const float p_mass, const float p_charge) {26 if (fk !=0 ) element_mat = rdipmat(element_length,fk,eloss,p_mass,p_charge);18 void H_RectangularDipole::setMatrix(const float eloss, const float p_mass, const float p_charge) const { 19 if (fk !=0 ) *element_mat = rdipmat(element_length,fk,eloss,p_mass,p_charge); 27 20 else { 28 element_mat = driftmat(element_length);29 if(VERBOSE) cout<<" <H_RectangularDipole>WARNING : k0= 0, drift-like dipole (" << name << ") !" << endl;21 *element_mat = driftmat(element_length); 22 if(VERBOSE) cout<<"\t WARNING : k0= 0, drift-like dipole (" << name << ") !" << endl; 30 23 } 31 24 return ; 32 25 } 33 34 H_RectangularDipole* H_RectangularDipole::clone() const {35 H_RectangularDipole* temp_dip = new H_RectangularDipole(name,fs,fk,element_length);36 temp_dip->setAperture(element_aperture);37 temp_dip->setX(xpos);38 temp_dip->setY(ypos);39 temp_dip->setTX(txpos);40 temp_dip->setTY(typos);41 temp_dip->setBetaX(betax);42 temp_dip->setBetaY(betay);43 return temp_dip;44 }45 -
trunk/external/Hector/H_RectangularDipole.h
r1360 r1365 2 2 #define _H_RectangularDipole_ 3 3 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * *5 * *6 * --<--<-- A fast simulator --<--<-- *7 * / --<--<-- of particle --<--<-- *8 * ----HECTOR----< *9 * \ -->-->-- transport through -->-->-- *10 * -->-->-- generic beamlines -->-->-- *11 * *12 * JINST 2:P09005 (2007) *13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) *14 * http://www.fynu.ucl.ac.be/hector.html *15 * *16 * Center for Cosmology, Particle Physics and Phenomenology *17 * Universite catholique de Louvain *18 * Louvain-la-Neuve, Belgium *19 * *20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */21 22 4 /// \file H_RectangularDipole.h 23 5 /// \brief Classes aiming at simulating LHC beam rectangular dipoles. 6 7 /* 8 ---- Hector the simulator ---- 9 A fast simulator of particles through generic beamlines. 10 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 11 12 http://www.fynu.ucl.ac.be/hector.html 13 14 Centre de Physique des Particules et de Phénoménologie (CP3) 15 Université Catholique de Louvain (UCL) 16 */ 24 17 25 18 #include "H_Dipole.h" … … 32 25 //@{ 33 26 H_RectangularDipole():H_Dipole(RDIPOLE,0.,0.,0.) {init();} 34 35 H_RectangularDipole(const string&nameE, const double s, const double k, const double l) :H_Dipole(nameE,RDIPOLE,s,k,l){init();}36 ~H_RectangularDipole() {};27 H_RectangularDipole(const double s, const double k, const double l) :H_Dipole(RDIPOLE,s,k,l){init();} 28 H_RectangularDipole(const string nameE, const double s, const double k, const double l) :H_Dipole(nameE,RDIPOLE,s,k,l){init();} 29 ~H_RectangularDipole() {return;}; 37 30 //@} 38 H_RectangularDipole* clone() const ;39 31 private: 40 41 virtual void setMatrix(const float, const float, const float);32 virtual void setTypeString() { typestring = RDIPOLENAME;}; 33 virtual void setMatrix(const float, const float, const float) const ; 42 34 }; 43 35 -
trunk/external/Hector/H_RomanPot.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 18 11 19 12 /// \file H_RomanPot.cc … … 32 25 } 33 26 34 H_RomanPot::H_RomanPot(const double s, const double app) :H_Drift(s,RP_LENGTH){ 35 type = RP; 27 H_RomanPot::H_RomanPot(const double s, const double app) :H_OpticalElement(RP,s,0.,RP_LENGTH){ 36 28 init(); 37 if(element_aperture) delete element_aperture;38 element_aperture = new H_RectangularAperture(app,RP_HEIGHT,0,0);29 H_RectangularAperture* rapp = new H_RectangularAperture(app,RP_HEIGHT,0,0); 30 setAperture(rapp); 39 31 } 40 32 41 H_RomanPot::H_RomanPot(const string& nameE, const double s, const double app) :H_Drift(nameE,s,RP_LENGTH){ 42 type = RP; 33 H_RomanPot::H_RomanPot(const string nameE, const double s, const double app) :H_OpticalElement(nameE,RP,s,0.,RP_LENGTH){ 43 34 init(); 44 if(element_aperture) delete element_aperture;45 element_aperture = new H_RectangularAperture(app,RP_HEIGHT,0,0);35 H_RectangularAperture* rapp = new H_RectangularAperture(app,RP_HEIGHT,0,0); 36 setAperture(rapp); 46 37 } 47 38 48 std::ostream& operator<< (std::ostream& os, const H_RomanPot& el) { 49 os << el.typestring << el.name << "\t\t at s = " << el.fs; 50 if(el.element_aperture->getType()!=NONE) { 51 os << *(el.element_aperture) << endl; 52 } 53 return os; 39 void H_RomanPot::printProperties() const { 40 cout << typestring << name; 41 cout << "\t\t at s = " << fs; 42 if(element_aperture->getType()!=NONE) { 43 cout <<"\t aperture type = " << element_aperture->getTypeString(); 44 element_aperture->printProperties(); 45 } 46 47 cout<<endl; 54 48 } 55 49 56 void H_RomanPot::setMatrix(const float eloss, const float p_mass, const float p_charge) {57 element_mat = driftmat(0);50 void H_RomanPot::setMatrix(const float eloss, const float p_mass, const float p_charge) const { 51 *element_mat = driftmat(0); 58 52 return ; 59 53 } 60 61 H_RomanPot* H_RomanPot::clone() const {62 H_RomanPot* temp_rp = new H_RomanPot(name,fs,element_length);63 temp_rp->setAperture(element_aperture);64 temp_rp->setX(xpos);65 temp_rp->setY(ypos);66 temp_rp->setTX(txpos);67 temp_rp->setTY(typos);68 temp_rp->setBetaX(betax);69 temp_rp->setBetaY(betay);70 return temp_rp;71 }72 -
trunk/external/Hector/H_RomanPot.h
r1360 r1365 1 1 #ifndef _H_RomanPot_ 2 2 #define _H_RomanPot_ 3 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * *5 * *6 * --<--<-- A fast simulator --<--<-- *7 * / --<--<-- of particle --<--<-- *8 * ----HECTOR----< *9 * \ -->-->-- transport through -->-->-- *10 * -->-->-- generic beamlines -->-->-- *11 * *12 * JINST 2:P09005 (2007) *13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) *14 * http://www.fynu.ucl.ac.be/hector.html *15 * *16 * Center for Cosmology, Particle Physics and Phenomenology *17 * Universite catholique de Louvain *18 * Louvain-la-Neuve, Belgium *19 * *20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */21 3 22 4 /// \file H_RomanPot.h 23 5 /// \brief Roman pot class 24 6 7 /* 8 ---- Hector the simulator ---- 9 A fast simulator of particles through generic beamlines. 10 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 11 12 http://www.fynu.ucl.ac.be/hector.html 13 14 Centre de Physique des Particules et de Phénoménologie (CP3) 15 Université Catholique de Louvain (UCL) 16 */ 17 25 18 // local #includes 26 #include "H_ Drift.h"19 #include "H_OpticalElement.h" 27 20 28 21 #define RP_LENGTH 0.0001 29 22 #define RP_HEIGHT 10000 30 23 /// Roman pot as an optics element. 31 class H_RomanPot : public H_ Drift {24 class H_RomanPot : public H_OpticalElement { 32 25 33 26 public: 34 27 /// Constructors and destructor 35 28 //@{ 36 H_RomanPot():H_ Drift() {type = RP;init();}37 H_RomanPot(const string &, const double, const double);29 H_RomanPot():H_OpticalElement(RP,0.,0.,RP_LENGTH) {init();} 30 H_RomanPot(const string, const double, const double); 38 31 H_RomanPot(const double, const double); 39 ~H_RomanPot() { };32 ~H_RomanPot() { return; }; 40 33 //@} 41 H_RomanPot* clone() const;34 virtual void printProperties() const; 42 35 void init(); 36 43 37 private: 44 38 virtual void setTypeString() {typestring=RPNAME;}; 45 virtual void setMatrix(const float, const float, const float) ;46 friend std::ostream& operator<< (std::ostream& os, const H_RomanPot& el); 39 virtual void setMatrix(const float, const float, const float) const; 40 47 41 }; 48 42 -
trunk/external/Hector/H_SectorDipole.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 18 11 19 12 /// \file H_SectorDipole.cc … … 23 16 #include "H_TransportMatrices.h" 24 17 25 void H_SectorDipole::setMatrix(const float eloss, const float p_mass, const float p_charge) {26 if (fk !=0 ) element_mat = sdipmat(element_length,fk,eloss,p_mass,p_charge);18 void H_SectorDipole::setMatrix(const float eloss, const float p_mass, const float p_charge) const { 19 if (fk !=0 ) *element_mat = sdipmat(element_length,fk,eloss,p_mass,p_charge); 27 20 else { 28 element_mat = driftmat(element_length);29 if(VERBOSE) cout<<" <H_SectorDipole>WARNING : k0= 0, drift-like dipole (" << name << ") !" << endl;21 *element_mat = driftmat(element_length); 22 if(VERBOSE) cout<<"\t WARNING : k0= 0, drift-like dipole (" << name << ") !" << endl; 30 23 } 31 24 return ; 32 25 } 33 34 H_SectorDipole* H_SectorDipole::clone() const {35 H_SectorDipole* temp_dip = new H_SectorDipole(name,fs,fk,element_length);36 temp_dip->setAperture(element_aperture);37 temp_dip->setX(xpos);38 temp_dip->setY(ypos);39 temp_dip->setTX(txpos);40 temp_dip->setTY(typos);41 temp_dip->setBetaX(betax);42 temp_dip->setBetaY(betay);43 return temp_dip;44 }45 -
trunk/external/Hector/H_SectorDipole.h
r1360 r1365 1 /// \file H_SectorDipole.h 2 /// \brief Classes aiming at simulating sector dipoles. 3 1 4 #ifndef _H_SectorDipole_ 2 5 #define _H_SectorDipole_ 3 6 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 5 * * 6 * --<--<-- A fast simulator --<--<-- * 7 * / --<--<-- of particle --<--<-- * 8 * ----HECTOR----< * 9 * \ -->-->-- transport through -->-->-- * 10 * -->-->-- generic beamlines -->-->-- * 11 * * 12 * JINST 2:P09005 (2007) * 13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 14 * http://www.fynu.ucl.ac.be/hector.html * 15 * * 16 * Center for Cosmology, Particle Physics and Phenomenology * 17 * Universite catholique de Louvain * 18 * Louvain-la-Neuve, Belgium * 19 * * 20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 7 /* 8 ---- Hector the simulator ---- 9 A fast simulator of particles through generic beamlines. 10 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 21 11 22 /// \file H_SectorDipole.h 23 /// \brief Classes aiming at simulating sector dipoles. 12 http://www.fynu.ucl.ac.be/hector.html 13 14 Centre de Physique des Particules et de Phénoménologie (CP3) 15 Université Catholique de Louvain (UCL) 16 */ 24 17 25 18 #include "H_Dipole.h" … … 33 26 H_SectorDipole():H_Dipole(SDIPOLE,0.,0.,0.) {init();} 34 27 H_SectorDipole(const double s, const double k, const double l) :H_Dipole(SDIPOLE,s,k,l){init();} 35 H_SectorDipole(const string &nameE, const double s, const double k, const double l) :H_Dipole(nameE,SDIPOLE,s,k,l){init();}36 ~H_SectorDipole() { };28 H_SectorDipole(const string nameE, const double s, const double k, const double l) :H_Dipole(nameE,SDIPOLE,s,k,l){init();} 29 ~H_SectorDipole() {return;}; 37 30 //@} 38 H_SectorDipole* clone() const ;39 31 private: 40 32 virtual void setTypeString() { typestring = SDIPOLENAME;}; 41 virtual void setMatrix(const float ,const float, const float) ;33 virtual void setMatrix(const float ,const float, const float) const ; 42 34 }; 43 35 -
trunk/external/Hector/H_TransportMatrices.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 18 11 19 12 /// \file H_TransportMatrices.cc … … 36 29 37 30 extern double omega(const double k, const double l) { 38 // [l] = [m] and [k] = [1/m ^2] for quadrupoles31 // [l] = [m] and [k] = [1/mï¿œ] for quadrupoles 39 32 // [omega] = [1] 40 33 return sqrt(fabs(k))*l; … … 42 35 43 36 extern double radius(const double k) { 44 // [k] = [1/m ^2] for quadrupoles37 // [k] = [1/mï¿œ] for quadrupoles 45 38 // [k] = [1/m] for dipoles 46 39 // [radius(k)] = [m] 47 if(k==0 && VERBOSE) cout<<" <H_TransportMatrices>ERROR : Dipole has no effect : results will be corrupted"<<endl;40 if(k==0 && VERBOSE) cout<<"ERROR : Dipole has no effect : results will be corrupted"<<endl; 48 41 // this is protected by the "if(k==0) -> driftmat" in every matrix below (ex vquatmat) 49 42 return (k==0) ? 1 : 1/k; 50 43 } 51 44 52 extern void printMatrix(TMatrix TMat) {45 extern void printMatrix(TMatrix * TMat) { 53 46 char temp[20]; 54 47 float * el = new float[MDIM*MDIM]; 55 el = (TMat .GetMatrixArray());48 el = (TMat->GetMatrixArray()); 56 49 57 50 cout << endl << "\t"; … … 156 149 double simp = r*2*sin(l/(2*r))*sin(l/(2*r))/BE; 157 150 double psy = ke*l/2.; 158 float tefmat[MDIM*MDIM] = {1., tan(psy)*ke, 0., 0., 0., 0.,151 float tefmat[MDIM*MDIM] = {1., (float)(tan(psy)*ke), 0., 0., 0., 0., 159 152 0., 1., 0., 0., 0., 0., 160 0., 0., 1., -tan(psy)*ke, 0., 0.,153 0., 0., 1., (float)(-tan(psy)*ke), 0., 0., 161 154 0., 0., 0., 1., 0., 0., 162 155 0., 0., 0., 0., 1., 0., … … 167 160 0.,0.,1.,0., 0., 0., 168 161 0.,0.,l,1., 0., 0., 169 simp, sin(l/r)/BE, 0., 0., 1., 0.,162 (float)simp, (float)(sin(l/r)/BE), 0., 0., 1., 0., 170 163 0., 0., 0., 0., 0., 1. }; 171 164 for(int i=0;i<MDIM*MDIM;i++) { … … 216 209 0.,0.,1.,0., 0., 0., 217 210 0.,0.,l,1., 0., 0., 218 simp, sin(l/r)/BE, 0., 0., 1., 0.,211 simp, (float)(sin(l/r)/BE), 0., 0., 1., 0., 219 212 0., 0., 0., 0., 0., 1. 220 213 }; … … 265 258 0.,0.,l ,1.,0.,0., 266 259 0.,0.,0.,0.,1.,0., 267 l*tan(ke)/2.,ke, 0., 0., 0., 1.260 (float)(l*tan(ke)/2.),ke, 0., 0., 0., 1. 268 261 }; 269 262 … … 292 285 0.,0.,l ,1.,0.,0., 293 286 0.,0.,0.,0.,1.,0., 294 0.,0., l*tan(ke)/2.,ke, 0., 1.287 0.,0.,(float)(l*tan(ke)/2.),ke, 0., 1. 295 288 }; 296 289 -
trunk/external/Hector/H_TransportMatrices.h
r1360 r1365 2 2 #define _H_TransportMatrices_ 3 3 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 5 * * 6 * --<--<-- A fast simulator --<--<-- * 7 * / --<--<-- of particle --<--<-- * 8 * ----HECTOR----< * 9 * \ -->-->-- transport through -->-->-- * 10 * -->-->-- generic beamlines -->-->-- * 11 * * 12 * JINST 2:P09005 (2007) * 13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 14 * http://www.fynu.ucl.ac.be/hector.html * 15 * * 16 * Center for Cosmology, Particle Physics and Phenomenology * 17 * Universite catholique de Louvain * 18 * Louvain-la-Neuve, Belgium * 19 * * 20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 4 /* 5 ---- Hector the simulator ---- 6 A fast simulator of particles through generic beamlines. 7 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 8 9 http://www.fynu.ucl.ac.be/hector.html 10 11 Centre de Physique des Particules et de Phénoménologie (CP3) 12 Université Catholique de Louvain (UCL) 13 */ 21 14 22 15 /** \file H_TransportMatrices.h … … 56 49 57 50 /// Prints the matrix 58 extern void printMatrix(TMatrix );51 extern void printMatrix(TMatrix * ); 59 52 60 53 /// \brief Returns the matrix for a vertically focussing quadrupole (H_VerticalQuadrupole) -
trunk/external/Hector/H_VerticalKicker.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 18 11 19 12 /// \file H_VerticalKicker.cc … … 26 19 #include "H_TransportMatrices.h" 27 20 28 void H_VerticalKicker::setMatrix(const float eloss, const float p_mass, const float p_charge) {21 void H_VerticalKicker::setMatrix(const float eloss, const float p_mass, const float p_charge) const { 29 22 extern int kickers_on; 30 23 if(kickers_on) { 31 element_mat = vkickmat(element_length,fk,eloss,p_mass,p_charge);24 *element_mat = vkickmat(element_length,fk,eloss,p_mass,p_charge); 32 25 } else{ 33 element_mat = driftmat(element_length);26 *element_mat = driftmat(element_length); 34 27 } 35 28 return ; 36 29 } 37 38 H_VerticalKicker* H_VerticalKicker::clone() const {39 H_VerticalKicker* temp_kick = new H_VerticalKicker(name,fs,fk,element_length);40 temp_kick->setAperture(element_aperture);41 temp_kick->setX(xpos);42 temp_kick->setY(ypos);43 temp_kick->setTX(txpos);44 temp_kick->setTY(typos);45 temp_kick->setBetaX(betax);46 temp_kick->setBetaY(betay);47 return temp_kick;48 }49 50 -
trunk/external/Hector/H_VerticalKicker.h
r1360 r1365 2 2 #define _H_VerticalKicker_ 3 3 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 5 * * 6 * --<--<-- A fast simulator --<--<-- * 7 * / --<--<-- of particle --<--<-- * 8 * ----HECTOR----< * 9 * \ -->-->-- transport through -->-->-- * 10 * -->-->-- generic beamlines -->-->-- * 11 * * 12 * JINST 2:P09005 (2007) * 13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 14 * http://www.fynu.ucl.ac.be/hector.html * 15 * * 16 * Center for Cosmology, Particle Physics and Phenomenology * 17 * Universite catholique de Louvain * 18 * Louvain-la-Neuve, Belgium * 19 * * 20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 4 /* 5 ---- Hector the simulator ---- 6 A fast simulator of particles through generic beamlines. 7 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 8 9 http://www.fynu.ucl.ac.be/hector.html 10 11 Centre de Physique des Particules et de Phénoménologie (CP3) 12 Université Catholique de Louvain (UCL) 13 */ 21 14 22 15 /// \file H_VerticalKicker.h … … 37 30 H_VerticalKicker():H_Kicker(VKICKER,0.,0.,0.) {init();} 38 31 H_VerticalKicker(const double s, const double k, const double l) :H_Kicker(VKICKER,s,k,l){init();} 39 H_VerticalKicker(const string &nameE, const double s, const double k, const double l) :H_Kicker(nameE,VKICKER,s,k,l){init();}40 ~H_VerticalKicker() { };32 H_VerticalKicker(const string nameE, const double s, const double k, const double l) :H_Kicker(nameE,VKICKER,s,k,l){init();} 33 ~H_VerticalKicker() {return;}; 41 34 //@} 42 H_VerticalKicker* clone() const ;43 35 private: 44 36 virtual void setTypeString() {typestring=VKICKERNAME;}; 45 virtual void setMatrix(const float, const float, const float) ;37 virtual void setMatrix(const float, const float, const float) const ; 46 38 }; 47 39 -
trunk/external/Hector/H_VerticalQuadrupole.cc
r1360 r1365 1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * * 3 * --<--<-- A fast simulator --<--<-- * 4 * / --<--<-- of particle --<--<-- * 5 * ----HECTOR----< * 6 * \ -->-->-- transport through -->-->-- * 7 * -->-->-- generic beamlines -->-->-- * 8 * * 9 * JINST 2:P09005 (2007) * 10 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 11 * http://www.fynu.ucl.ac.be/hector.html * 12 * * 13 * Center for Cosmology, Particle Physics and Phenomenology * 14 * Universite catholique de Louvain * 15 * Louvain-la-Neuve, Belgium * 16 * * 17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 1 /* 2 ---- Hector the simulator ---- 3 A fast simulator of particles through generic beamlines. 4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 5 6 http://www.fynu.ucl.ac.be/hector.html 7 8 Centre de Physique des Particules et de Phénoménologie (CP3) 9 Université Catholique de Louvain (UCL) 10 */ 18 11 19 12 /// \file H_VerticalQuadrupole.cc … … 23 16 #include "H_VerticalQuadrupole.h" 24 17 25 void H_VerticalQuadrupole::setMatrix(const float eloss, const float p_mass, const float p_charge) {26 if (fk<0) { if(VERBOSE) cout<<" <H_VerticalQuadrupole>ERROR : k1 should be > 0 for H_VerticalQuadrupole (" << name << ")!"<<endl; }27 if (fk !=0 ) element_mat = vquadmat(element_length,fk,eloss, p_mass,p_charge);18 void H_VerticalQuadrupole::setMatrix(const float eloss, const float p_mass, const float p_charge) const { 19 if (fk<0) { if(VERBOSE) cout<<"\t ERROR : k1 should be > 0 for H_VerticalQuadrupole (" << name << ")!"<<endl; } 20 if (fk !=0 ) *element_mat = vquadmat(element_length,fk,eloss, p_mass,p_charge); 28 21 else { 29 element_mat = driftmat(element_length);30 if(VERBOSE) cout<<" <H_VerticalQuadrupole>WARNING : k1= 0, drift-like quadrupole (" << name << ") !" << endl;22 *element_mat = driftmat(element_length); 23 if(VERBOSE) cout<<"\t WARNING : k1= 0, drift-like quadrupole (" << name << ") !" << endl; 31 24 } 32 25 return ; 33 26 } 34 35 H_VerticalQuadrupole* H_VerticalQuadrupole::clone() const {36 H_VerticalQuadrupole* temp_quad = new H_VerticalQuadrupole(name,fs,fk,element_length);37 temp_quad->setAperture(element_aperture);38 temp_quad->setX(xpos);39 temp_quad->setY(ypos);40 temp_quad->setTX(txpos);41 temp_quad->setTY(typos);42 temp_quad->setBetaX(betax);43 temp_quad->setBetaY(betay);44 return temp_quad;45 } -
trunk/external/Hector/H_VerticalQuadrupole.h
r1360 r1365 2 2 #define _H_VerticalQuadrupole_ 3 3 4 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * 5 * * 6 * --<--<-- A fast simulator --<--<-- * 7 * / --<--<-- of particle --<--<-- * 8 * ----HECTOR----< * 9 * \ -->-->-- transport through -->-->-- * 10 * -->-->-- generic beamlines -->-->-- * 11 * * 12 * JINST 2:P09005 (2007) * 13 * X Rouby, J de Favereau, K Piotrzkowski (CP3) * 14 * http://www.fynu.ucl.ac.be/hector.html * 15 * * 16 * Center for Cosmology, Particle Physics and Phenomenology * 17 * Universite catholique de Louvain * 18 * Louvain-la-Neuve, Belgium * 19 * * 20 * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 4 /* 5 ---- Hector the simulator ---- 6 A fast simulator of particles through generic beamlines. 7 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be 21 8 9 http://www.fynu.ucl.ac.be/hector.html 10 11 Centre de Physique des Particules et de Phénoménologie (CP3) 12 Université Catholique de Louvain (UCL) 13 */ 22 14 23 15 /// \file H_VerticalQuadrupole.h … … 35 27 H_VerticalQuadrupole():H_Quadrupole(VQUADRUPOLE,0.,0.,0.) {init();} 36 28 H_VerticalQuadrupole(const double s, const double k, const double l):H_Quadrupole(VQUADRUPOLE,s,k,l) {init();} 37 H_VerticalQuadrupole(const string& nameE, const double s, const double k, const double l):H_Quadrupole(nameE,VQUADRUPOLE,s,k,l) {init();} 38 ~H_VerticalQuadrupole() {}; 39 //@} 40 H_VerticalQuadrupole* clone() const ; 41 private: 29 H_VerticalQuadrupole(string nameE, const double s, const double k, const double l):H_Quadrupole(nameE,VQUADRUPOLE,s,k,l) {init();} 30 ~H_VerticalQuadrupole() {return;}; 31 //@} 32 private: 42 33 virtual void setTypeString() {typestring = VQUADRUPOLENAME;} ; 43 virtual void setMatrix(const float, const float, const float) ;34 virtual void setMatrix(const float, const float, const float) const; 44 35 }; 45 36
Note:
See TracChangeset
for help on using the changeset viewer.