Fork me on GitHub

Changeset 405 in svn for trunk/src


Ignore:
Timestamp:
May 20, 2009, 10:39:01 AM (16 years ago)
Author:
Xavier Rouby
Message:

bug removed in Hector implementation for RP220/FP420

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/VeryForward.cc

    r389 r405  
    3131
    3232#include "VeryForward.h"
     33#include "PdgParticle.h"
    3334#include "H_RomanPot.h"
    34 #include "PdgParticle.h"
     35
    3536#include <iostream>
     37#include <fstream>
    3638#include<cmath>
    3739
    3840using namespace std;
    3941
     42/* Notes on the correct initialisation for Hector
     43 * -- these notes apply to the LHC beamlines
     44 *
     45 *  beam1 : forward direction is for increasing 's' values
     46 *    beamline1 = new H_BeamLine(1,...);
     47 *    beamline1->fill(DET->RP_beam1Card,1,DET->RP_IP_name);
     48 *  beam2 : forward direction is for decreasing 's' values
     49 *    beamline2 = new H_BeamLine(-1,...);
     50 *    beamline2->fill(DET->RP_beam2Card,-1,DET->RP_IP_name);
     51 *
     52 *  relative_energy should be false -- kickers_on should be 1
     53 *
     54 */
    4055
    4156//------------------------------------------------------------------------------
    4257VeryForward::VeryForward() :
    4358   DET(new RESOLution()), d_max(1.+std::max(DET->RP_420_s,DET->RP_220_s)),
    44    beamline1(new H_BeamLine(1,d_max)), beamline2(new H_BeamLine(1,d_max)),
    45    relative_energy(true), // should always be true
    46    kickers_on(1)         // should always be 1
     59   beamline1(new H_BeamLine(1,d_max)), beamline2(new H_BeamLine(-1,d_max)),
     60   rel_energy(true), // should always be true
     61   kickers(1)         // should always be 1
    4762   {
    4863   init(); //Initialisation of Hector
    … …  
    5570   const float d_max = 1.+std::max(DET->RP_420_s,DET->RP_220_s);
    5671   beamline1 = new H_BeamLine(1,d_max);
    57    beamline2 = new H_BeamLine(1,d_max);
     72   beamline2 = new H_BeamLine(-1,d_max);
    5873   init(); //Initialisation of Hector
    59    relative_energy = true; // should always be true
    60    kickers_on = 1;         // should always be 1
     74   rel_energy = true; // should always be true
     75   kickers = 1;         // should always be 1
    6176}
    6277
    6378VeryForward::VeryForward(const RESOLution * DetDatacard) :
    6479        DET(new RESOLution(*DetDatacard)), d_max(1.+std::max(DET->RP_420_s,DET->RP_220_s)),
    65         beamline1(new H_BeamLine(1,d_max)), beamline2(new H_BeamLine(1,d_max)),
    66         relative_energy(true), // should always be true
    67         kickers_on(1)          // should always be 1
     80        beamline1(new H_BeamLine(1,d_max)), beamline2(new H_BeamLine(-1,d_max)),
     81        rel_energy(true), // should always be true
     82        kickers(1)          // should always be 1
    6883   {
    6984   init();  //Initialisation of Hector
    … …  
    7388   DET(new RESOLution(*(vf.DET))),   d_max(vf.d_max),
    7489   beamline1(new H_BeamLine(*(vf.beamline1))), beamline2(new H_BeamLine(*(vf.beamline2))),
    75    relative_energy(vf.relative_energy),
    76    kickers_on(vf.kickers_on) {
     90   rel_energy(vf.rel_energy),
     91   kickers(vf.kickers) {
    7792}
    7893
    … …  
    8398   beamline1 = new H_BeamLine(*(vf.beamline1));
    8499   beamline2 = new H_BeamLine(*(vf.beamline2));
    85    relative_energy =vf.relative_energy;
    86    kickers_on = vf.kickers_on;
     100   rel_energy =vf.rel_energy;
     101   kickers = vf.kickers;
    87102   return *this;
    88103}
    … …  
    93108  static unsigned int counter;
    94109  counter =0;
    95   relative_energy = true; // should always be true
    96   kickers_on = 1;         // should always be 1
     110  extern bool relative_energy;
     111  extern int kickers_on;
     112  relative_energy = rel_energy; // should always be true
     113  kickers_on = kickers;         // should always be 1
    97114  beamline1->fill(DET->RP_beam1Card,1,DET->RP_IP_name);                               
    98   beamline1->offsetElements(DET->RP_offsetEl_s,-DET->RP_offsetEl_x);
    99   H_RomanPot * rp220_1 = new H_RomanPot("rp220_1",DET->RP_220_s,DET->RP_220_x*1E6); // RP 220m, 2mm, beam 1 
    100   H_RomanPot * rp420_1 = new H_RomanPot("rp420_1",DET->RP_420_s,DET->RP_420_x*1E6); // RP 420m, 4mm, beam 1 
     115  beamline1->offsetElements(DET->RP_offsetEl_s,-DET->RP_offsetEl_x); // relative energy: does not change anything
     116  H_RomanPot * rp220_1 = new H_RomanPot("rp220_1",DET->RP_220_s,DET->RP_220_x*1E6);
     117  // RP 220m, 2mm, beam 1 
     118  H_RomanPot * rp420_1 = new H_RomanPot("rp420_1",DET->RP_420_s,DET->RP_420_x*1E6);
     119  // RP 420m, 4mm, beam 1 
     120  //rp220_1->printProperties();
     121  //rp420_1->printProperties();
    101122  beamline1->add(rp220_1);
    102123  beamline1->add(rp420_1);
    103124 
    104125  beamline2->fill(DET->RP_beam2Card,-1,DET->RP_IP_name);                             
    105   beamline2->offsetElements(DET->RP_offsetEl_s,+DET->RP_offsetEl_x);
    106   H_RomanPot * rp220_2 = new H_RomanPot("rp220_2",DET->RP_220_s,DET->RP_220_x*1E6);// RP 220m, 2mm, beam 2     
    107   H_RomanPot * rp420_2 = new H_RomanPot("rp420_2",DET->RP_420_s,DET->RP_420_x*1E6);// RP 420m, 4mm, beam 2     
     126  beamline2->offsetElements(DET->RP_offsetEl_s,-DET->RP_offsetEl_x); // relative energy: does not change anything
     127  H_RomanPot * rp220_2 = new H_RomanPot("rp220_2",DET->RP_220_s,DET->RP_220_x*1E6);
     128  // RP 220m, 2mm, beam 2     
     129  H_RomanPot * rp420_2 = new H_RomanPot("rp420_2",DET->RP_420_s,DET->RP_420_x*1E6);
     130  // RP 420m, 4mm, beam 2     
     131  //rp220_2->printProperties();
     132  //rp420_2->printProperties();
    108133  beamline2->add(rp220_2);
    109134  beamline2->add(rp420_2);
    110135  // rp220_1, rp220_2, rp420_1 and rp420_2 will be deallocated in ~H_AbstractBeamLine
    111136  // do not put explicit delete
     137}
     138
     139
     140float VeryForward::time_of_flight(TRootGenParticle *particle, const float detector_s, const float detector_etamin, const float detector_t_resolution) {
     141    // time of flight t is t = T + d/[ cos(theta) v ]
     142    float cos_theta = 1;     //very good approximation, if detector_etamin >3
     143    if (detector_etamin<3) { // if smaller eta -> make the complete calculation
     144       double tx = atan(particle->Px/particle->Pz);
     145       double ty = atan(particle->Py/particle->Pz);
     146       double theta = sqrt( pow(tx,2) + pow(ty,2) );
     147       // cout << "tx = " << tx << " ty = " << ty << " theta  = " << theta << " cos(theta) = " << cos(theta) << endl;
     148       // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1
     149       // eta = 2.6 <-> cos(theta) = 0.99
     150       // eta = 3.0 <-> cos(theta) = 0.995
     151       cos_theta = cos(theta);
     152    }
     153    // units from StdHEP : Z [mm] T[mm/c]
     154    // units from Delphes : VFD_s_zdc [m] speed_of_light [m/s]
     155    double flight_distance = (detector_s - particle->Z*(1E-3))/cos_theta ;
     156    // assumes highly relativistic particles
     157    double flight_time = (flight_distance + 1E-3 * particle->T )/speed_of_light;
     158    double timeS = gRandom->Gaus(flight_time,detector_t_resolution);
     159    return timeS;
    112160}
    113161
    … …  
    125173    elementZdc = (TRootZdcHits*) branchZDC->NewEntry();
    126174
    127 
    128     elementZdc->pid = particle->PID;
    129 
    130      
    131     // for compatibility with 'old' version
    132175    TLorentzVector genMomentum;
    133176    genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
    134     elementZdc->Set(genMomentum);
    135     // ******************
    136    
    137  
    138     //particle->print();
     177    elementZdc->Set(genMomentum);  // initialises the gen-level data
     178    elementZdc->pid = particle->PID;
    139179 
    140180    // 1) energy smearing
    … …  
    149189                                                    pow(DET->HAD_Szdc*sqrt(particle->E),2) ));
    150190    elementZdc->E = energyS;
    151  
    152 
    153     // 2) time of flight t is t = T + d/[ cos(theta) v ]
    154     float cos_theta = 1;      //very good approximation, if eta_zdc >3
    155     if (DET->VFD_min_zdc<3) { // if smaller eta -> make the complete calculation
    156        double tx = atan(particle->Px/particle->Pz);
    157        double ty = atan(particle->Py/particle->Pz);
    158        double theta = sqrt( pow(tx,2) + pow(ty,2) );
    159        //cout << "tx = " << tx << " ty = " << ty << " theta  = " << theta << " cos(theta) = " << cos(theta) << endl;
    160        // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1
    161        // eta = 2.6 <-> cos(theta) = 0.99
    162        // eta = 3.0 <-> cos(theta) = 0.995
    163        cos_theta = cos(theta);
    164     }
    165     // units from StdHEP : Z [mm] T[mm/c]
    166     // units from Delphes : VFD_s_zdc [m] speed_of_light [m/s]
    167     double flight_distance = (DET->VFD_s_zdc - particle->Z*(1E-3))/cos_theta ;
    168     double flight_time = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s]
    169     double timeS = gRandom->Gaus(flight_time,DET->ZDC_T_resolution);
    170     elementZdc->T = timeS;
     191
     192    // 2) time of flight t is t = T + d/[ cos(theta) v ] + detector smearing on time
     193    elementZdc->T = time_of_flight(particle, DET->VFD_s_zdc, DET->VFD_min_zdc, DET->ZDC_T_resolution);
    171194
    172195    // 3) side: which ZDC has been hit?
    … …  
    174197
    175198    // 4) object nature : e.m. (photon) or had (neutron) ?
    176     //elementZdc->hadronic_hit = (bool) (particle->PID==pN);
    177   }
     199    elementZdc->hadronic_hit = (bool) (particle->PID!=pGAMMA);
     200  } // if neutrons or photons over E_threshold
     201
    178202}
    179203
    … …  
    181205void VeryForward::RomanPots(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchRP220,ExRootTreeBranch *branchFP420,TRootGenParticle *particle)
    182206{
     207  if(particle->Status != 1) return; // reject particles that are not final ones
     208  extern bool relative_energy;
     209  relative_energy = rel_energy;
     210  extern int kickers_on;
     211  kickers_on = kickers;
     212
    183213  float charge = particle->Charge, mass = particle->M;
     214  //float charge, mass, ctau;
     215  //charge = mass = ctau = UNDEFINED;
    184216  if (mass<-999) { // unitialised!
    185217          PdgParticle pdg_part = DET->PDGtable[particle->PID];
    186           charge  = pdg_part.charge();
    187           mass    = pdg_part.mass();
     218          charge  = pdg_part.charge();  // e+
     219          mass    = pdg_part.mass();    // GeV
     220  //        ctau    = pdg_part.ctau();    // m
     221  // cout << "ctau = " << ctau << endl;
    188222  }
    189   //if(particle->Charge!=1) return; // only particles with Q=+1 can hope to reach RP200/FP420
    190  
     223
     224
     225  if(particle->Charge==0) return; // only particles with Q=+1 can hope to reach RP200/FP420
     226  //cout << "particle ("<< particle->PID << "): m = " << mass << " \t Q= " << charge << endl;
    191227  TRootRomanPotHits* elementRP220;
    192228  //TRootForwardTaggerHits* elementFP420;
    … …  
    197233 
    198234  // to go faster, why not rejecting particles already going into the ZDC?
    199   if( particle->PID == pP)
    200   if( (particle->Status == 1) &&  (fabs(genMomentum.Eta()) > DET->CEN_max_calo_fwd) )
     235 
     236  // K_L^0 has a ctau of ~16m ; only pi+ and p+ can beat it ;
     237  // so if RP/FP too far away, the particle must be a proton or a mu+
     238  if( std::min(DET->RP_420_s,DET->RP_220_s) > 17 && (particle->PID != pP && particle->PID != 13)) return;
     239
     240  if( fabs(genMomentum.Eta()) > DET->CEN_max_calo_fwd )
    201241    {
    202       //cout << "VeryForward :: M = " << mass << "\t Q = " << charge << "\t\t " << particle->PID << endl;
    203242      H_BeamParticle p1(mass,charge);
    204       p1.smearAng();   p1.smearPos(); // vertex smearing
    205         cout << "x = " << p1.getX() + DET->RP_cross_x
    206              << " y= " << p1.getY() + DET->RP_cross_y
    207              << " tx= " << p1.getTX() - kickers_on*DET->RP_cross_ang
    208              << " ty=" << p1.getTY() << endl;
    209       p1.setPosition(p1.getX()+DET->RP_cross_x,p1.getY()+DET->RP_cross_y,p1.getTX()-kickers_on*DET->RP_cross_ang,p1.getTY(),0);
    210       //p1.set4Momentum(particle->Px,particle->Py,particle->Pz,particle->E);
     243      double tx = 1E6*atan(particle->Px/particle->Pz); // in microrad
     244      double ty = 1E6*atan(particle->Py/particle->Pz); // in microrad
     245      //p1.smearAng();   p1.smearPos(); // vertex smearing do not put it here !!!
     246
     247      /* cout << "x = "  << particle->X << " + " << p1.getX()   << " + " <<  DET->RP_cross_x
     248             << " y= "  << particle->Y << " + " << p1.getY()   << " + " << DET->RP_cross_y
     249             << " tx= " << tx << " + " << p1.getTX() << " - " << kickers_on*DET->RP_cross_ang_x
     250             << " ty= " << ty << " + " << p1.getTY()  << " - " << kickers_on*DET->RP_cross_ang_y 
     251             << " z= "  << particle->Z << endl;*/
     252
     253      // here below, p1.getX(), p1.getY(), p1.getTX() and p1.getTY() =0 unless some smearing is done
     254      // all in micrometers or microradians
     255      p1.setPosition((1E3)*particle->X + p1.getX() + DET->RP_cross_x,
     256                     (1E3)*particle->Y + p1.getY() + DET->RP_cross_y,
     257                     tx + p1.getTX()- kickers_on*DET->RP_cross_ang_x,
     258                     ty + p1.getTY()- kickers_on*DET->RP_cross_ang_y,
     259                     (1E3)*particle->Z);
    211260      p1.setE(particle->E);
    212261       
    … …  
    218267     
    219268      if(p1.stopped(beamline)) {
     269        // roman pots at 220 m
    220270        if (p1.getStoppingElement()->getName()=="rp220_1" || p1.getStoppingElement()->getName()=="rp220_2") {
    221           static unsigned int counter;
    222           counter++;
    223           if (counter==1) {
    224                 //p1.getPath(0,"p1path.txt");
    225                 cout << "RP : " << particle->PID << "\t" << charge << "=" << particle->Charge
    226                      << "\t" << mass << "=" << particle->M
    227                      << "\t E=" << particle->E
    228                      << endl;           
    229           }
     271
     272        /*static unsigned int counter;
     273        counter++;
     274        if (counter==1) {
     275           p1.getPath(0,"p1path.txt");
     276           cout << "RP : " << particle->PID << "\t" << charge << "=" << particle->Charge
     277                << "\t" << mass << "=" << particle->M << "\t E=" << particle->E << endl;
     278          }*/
     279
    230280          p1.propagate(DET->RP_220_s);
    231281          elementRP220 = (TRootRomanPotHits*) branchRP220->NewEntry();
     282         
     283          // detector measurements
    232284          elementRP220->X  = (1E-6)*p1.getX();  // [m]
    233285          elementRP220->Y  = (1E-6)*p1.getY();  // [m]
    … …  
    235287          elementRP220->Ty = (1E-6)*p1.getTY(); // [rad]
    236288          elementRP220->S = p1.getS();          // [m]
    237 
    238             /* time of flight t is t = T + d/[ cos(theta) v ]
    239             // nb: here we assume a straight path to the detector, which is not the case!
    240             // this time estimate is always underestimated (while exact for the ZDC case)
    241             float cos_theta = 1;      //very good approximation, if CEN_max_calo_fwd >3
    242             if (DET->CEN_max_calo_fwd<3) { // if smaller eta -> make the complete calculation
    243                double tx = atan(particle->Px/particle->Pz);
    244                double ty = atan(particle->Py/particle->Pz);
    245                double theta = sqrt( pow(tx,2) + pow(ty,2) );
    246                //cout << "tx = " << tx << " ty = " << ty << " theta  = " << theta << " cos(theta) = " << cos(theta) << endl;
    247                // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1
    248                // eta = 2.6 <-> cos(theta) = 0.99
    249                // eta = 3.0 <-> cos(theta) = 0.995
    250                cos_theta = cos(theta);
    251             }
    252             // units from StdHEP : Z [mm] T[mm/c]
    253             // units from Delphes : p1.getS [m]   speed_of_light [m/s]
    254             //double flight_distance = (p1.getS() - particle->Z*(1E-3))/cos_theta ;
    255             //elementRP220->T = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s]
    256             */
    257             elementRP220->E = p1.getE();        // not yet implemented
    258             elementRP220->q2 = -1;                // not yet implemented
    259             elementRP220->side = sign(particle->Eta);
    260 
    261 
    262                 elementRP220->pid = particle->PID;
     289          elementRP220->T = time_of_flight(particle, DET->RP_220_s, DET->CEN_max_calo_fwd, DET->RP220_T_resolution);
     290          elementRP220->side = sign(particle->Eta);
    263291         
    264         } else if (p1.getStoppingElement()->getName()=="rp420_1" || p1.getStoppingElement()->getName()=="rp420_2") {
     292          // reconstructed data 
     293          float sE = p1.getE(); // apply the smearing here!!!
     294          elementRP220->E = sE;           // not yet implemented
     295          elementRP220->q2 = UNDEFINED;   // not yet implemented
     296
     297          // generator level data
     298          elementRP220->pid = particle->PID;
     299          elementRP220->Set(genMomentum);
     300        } // if RP220
     301
     302        // proton taggers at 420 m
     303        else if (p1.getStoppingElement()->getName()=="rp420_1" || p1.getStoppingElement()->getName()=="rp420_2") {
     304
    265305          p1.propagate(DET->RP_420_s);
     306          elementFP420 = (TRootRomanPotHits*) branchFP420->NewEntry();
    266307          //elementFP420 = (TRootForwardTaggerHits*) branchFP420->NewEntry();
    267           elementFP420 = (TRootRomanPotHits*) branchFP420->NewEntry();
     308
     309          // detector measurements
    268310          elementFP420->X  = (1E-6)*p1.getX();  // [m]
    269311          elementFP420->Y  = (1E-6)*p1.getY();  // [m]
    … …  
    271313          elementFP420->Ty = (1E-6)*p1.getTY(); // [rad]
    272314          elementFP420->S = p1.getS();          // [m]
    273 
    274             // time of flight t is t = T + d/[ cos(theta) v ]
    275             // nb: here we assume a straight path to the detector, which is not the case!
    276             // this time estimate is always underestimated (while exact for the ZDC case)
    277             float cos_theta = 1;      //very good approximation, if CEN_max_calo_fwd >3
    278             if (DET->CEN_max_calo_fwd<3) { // if smaller eta -> make the complete calculation
    279                double tx = atan(particle->Px/particle->Pz);
    280                double ty = atan(particle->Py/particle->Pz);
    281                double theta = sqrt( pow(tx,2) + pow(ty,2) );
    282                //cout << "tx = " << tx << " ty = " << ty << " theta  = " << theta << " cos(theta) = " << cos(theta) << endl;
    283                // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1
    284                // eta = 2.6 <-> cos(theta) = 0.99
    285                // eta = 3.0 <-> cos(theta) = 0.995
    286                cos_theta = cos(theta);
    287             }
    288             // units from StdHEP : Z [mm] T[mm/c]
    289             // units from Delphes : p1.getS [m]   speed_of_light [m/s]
    290             double flight_distance = (p1.getS() - particle->Z*(1E-3))/cos_theta ;
    291             elementFP420->T = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s]
     315          elementFP420->T = time_of_flight(particle, DET->RP_420_s, DET->CEN_max_calo_fwd, DET->RP420_T_resolution);
     316          elementFP420->side = sign(particle->Eta);
     317
     318          // reconstructed data
    292319          elementFP420->E = p1.getE();          // not yet implemented
    293           elementFP420->q2 = -1;                // not yet implemented
    294           elementFP420->side = sign(particle->Eta);
    295 
    296 
    297                 elementFP420->pid = particle->PID;
    298         }
    299 
    300       }
    301       //  if(p1.stopped(beamline) && (p1.getStoppingElement()->getS() > 100))
    302       //     cout << "Eloss ="  << 7000.-p1.getE() << " ; " << p1.getStoppingElement()->getName() << endl;
     320          elementFP420->q2 = UNDEFINED;         // not yet implemented
     321
     322          // generator level data
     323          elementFP420->pid = particle->PID;
     324          elementFP420->Set(genMomentum);
     325        }  // if FP420
     326
     327      } // if stopped
    303328    } // if forward proton
    304 }
    305 
    306     // Forward particles in CASTOR ?
    307     //    if (particle->Status == 1 && (fabs(particle->Eta) > DET->MIN_CALO_VFWD)
    308     //                              && (fabs(particle->Eta) < DET->MAX_CALO_VFWD)) {
    309     //
    310     //
    311     //    } // CASTOR
    312     // */
    313 
     329
     330}
Note: See TracChangeset for help on using the changeset viewer.