Fork me on GitHub

Ignore:
Timestamp:
Apr 16, 2014, 3:56:14 PM (11 years ago)
Author:
Pavel Demin
Message:

switch to a more stable Hector version

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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*/
    1811
    1912/// \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"
    2323
    2424// local #includes
    … …  
    2828using namespace std;
    2929
    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  {}
     30H_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
     52H_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
     83H_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}
    94104
    95105H_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
     129float 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
     178void 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;
    122189       
    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
     193void 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
     206float 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
     214float 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
     224float 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
     237float 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
     250float 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
     264float 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
     278float 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
     293float 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;
    148308        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
     323float H_RecRPObject::computeE_PM() {
     324        cout<<"Not yet implemented, nothing done"<<endl;
    220325        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
     329float 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
Note: See TracChangeset for help on using the changeset viewer.