Fork me on GitHub

Changeset 374 in svn for trunk/src/VeryForward.cc


Ignore:
Timestamp:
May 10, 2009, 8:23:32 PM (15 years ago)
Author:
Xavier Rouby
Message:

added Resolution terms for energy and timing for ZDC

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/VeryForward.cc

    r361 r374  
    44**         |  Delphes, a framework for the fast simulation  |         **
    55**         |       of a  generic collider experiment        |         **
    6 **          \----------------------------------------------/          **
     6**          \-------------  arXiv:0903.2225v1  ------------/          **
    77**                                                                    **
    88**                                                                    **
    … …  
    3939
    4040//------------------------------------------------------------------------------
    41 
    42 VeryForward::VeryForward() {
    43    DET = new RESOLution();
    44    beamline1 = new H_BeamLine(1,500.);
    45    beamline2 = new H_BeamLine(1,500.);
    46    init();
    47   //Initialisation of Hector
    48   relative_energy = true; // should always be true
    49   kickers_on = 1;         // should always be 1
    50 
    51 }
    52 
    53 VeryForward::VeryForward(const string& DetDatacard) {
    54    DET = new RESOLution();
     41VeryForward::VeryForward() :
     42   DET(new RESOLution()), d_max(1.+std::max(DET->RP_420_s,DET->RP_220_s)),
     43   beamline1(new H_BeamLine(1,d_max)), beamline2(new H_BeamLine(1,d_max)),
     44   relative_energy(true), // should always be true
     45   kickers_on(1)         // should always be 1
     46   {
     47   init(); //Initialisation of Hector
     48}
     49
     50VeryForward::VeryForward(const string& DetDatacard) :
     51           DET(new RESOLution())
     52   {
    5553   DET->ReadDataCard(DetDatacard);
    56    beamline1 = new H_BeamLine(1,500.);
    57    beamline2 = new H_BeamLine(1,500.);
    58    init();
    59   //Initialisation of Hector
    60   relative_energy = true; // should always be true
    61   kickers_on = 1;         // should always be 1
    62 
    63 }
    64 
    65 VeryForward::VeryForward(const RESOLution * DetDatacard) {
    66    DET = new RESOLution(*DetDatacard);
    67    beamline2 = new H_BeamLine(1,500.);
    68    beamline1 = new H_BeamLine(1,500.);
    69 
    70    init();
    71   //Initialisation of Hector
    72   relative_energy = true; // should always be true
    73   kickers_on = 1;         // should always be 1
    74 
    75 }
    76 
    77 VeryForward::VeryForward(const VeryForward& vf) {
    78    DET = new RESOLution(*(vf.DET));
    79    beamline1 = new H_BeamLine(*(vf.beamline1));
    80    beamline2 = new H_BeamLine(*(vf.beamline2));
     54   const float d_max = 1.+std::max(DET->RP_420_s,DET->RP_220_s);
     55   beamline1 = new H_BeamLine(1,d_max);
     56   beamline2 = new H_BeamLine(1,d_max);
     57   init(); //Initialisation of Hector
     58   relative_energy = true; // should always be true
     59   kickers_on = 1;         // should always be 1
     60}
     61
     62VeryForward::VeryForward(const RESOLution * DetDatacard) :
     63        DET(new RESOLution(*DetDatacard)), d_max(1.+std::max(DET->RP_420_s,DET->RP_220_s)),
     64        beamline1(new H_BeamLine(1,d_max)), beamline2(new H_BeamLine(1,d_max)),
     65        relative_energy(true), // should always be true
     66        kickers_on(1)          // should always be 1
     67   {
     68   init();  //Initialisation of Hector
     69}
     70
     71VeryForward::VeryForward(const VeryForward& vf) :
     72   DET(new RESOLution(*(vf.DET))),   d_max(vf.d_max),
     73   beamline1(new H_BeamLine(*(vf.beamline1))), beamline2(new H_BeamLine(*(vf.beamline2))),
     74   relative_energy(vf.relative_energy),
     75   kickers_on(vf.kickers_on) {
    8176}
    8277
    … …  
    8479   if (this==&vf) return *this;
    8580   DET = new RESOLution(*(vf.DET));
     81   d_max = vf.d_max;
    8682   beamline1 = new H_BeamLine(*(vf.beamline1));
    8783   beamline2 = new H_BeamLine(*(vf.beamline2));
     84   relative_energy =vf.relative_energy;
     85   kickers_on = vf.kickers_on;
    8886   return *this;
    8987}
    … …  
    9492  relative_energy = true; // should always be true
    9593  kickers_on = 1;         // should always be 1
    96   // user should provide : (1) optics file for each beamline, and IPname,
    97   // and offset data (s,x) for optical elements
    9894  beamline1->fill(DET->RP_beam1Card,1,DET->RP_IP_name);                               
    9995  beamline1->offsetElements(DET->RP_offsetEl_s,-DET->RP_offsetEl_x);
    … …  
    114110
    115111
    116 void VeryForward::ZDC(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchZDC,TRootGenParticle *particle)
     112void VeryForward::ZDC(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchZDC, TRootGenParticle *particle)
    117113{
    118   int pid=abs(particle->PID);
    119114  TRootZdcHits *elementZdc;
    120   TLorentzVector genMomentum;
     115  float energy = particle->E;
     116  //TLorentzVector genMomentum;
     117
    121118  // Zero degree calorimeter, for forward neutrons and photons
    122   if (particle->Status ==1 && (pid == pN || pid == pGAMMA ) && fabs(particle->Eta) > DET->VFD_min_zdc ) {
    123     genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
    124     // !!!!!!!!! vérifier que particle->Z est bien en micromÚtres!!!
    125     // !!!!!!!!! vérifier que particle->T est bien en secondes!!!
    126     // !!!!!!!!! pas de smearing ! on garde trop d'info !
     119  if (particle->Status ==1 && ( (particle->PID==pN     && energy>DET->ZDC_n_E) ||
     120                                (particle->PID==pGAMMA && energy>DET->ZDC_gamma_E) )
     121                           && fabs(particle->Eta) > DET->VFD_min_zdc ) {
    127122    elementZdc = (TRootZdcHits*) branchZDC->NewEntry();
    128     elementZdc->Set(genMomentum);
    129    
    130     // time of flight t is t = T + d/[ cos(theta) v ]
    131     //double tx = acos(particle->Px/particle->Pz);
    132     //double ty = acos(particle->Py/particle->Pz);
    133     //double theta = (1E-6)*sqrt( pow(tx,2) + pow(ty,2) );
    134     //double flight_distance = (DET->ZDC_S - particle->Z*(1E-6))/cos(theta) ; // assumes that Z is in micrometers
    135     double flight_distance = DET->VFD_s_zdc - particle->Z*(1E-6);
    136     // assumes also that the emission angle is so small that 1/(cos theta) = 1
    137     elementZdc->T = particle->T + flight_distance/speed_of_light; // assumes highly relativistic particles
    138 //cout << "ZDC: T = " << particle->T << " ; " << flight_distance/speed_of_light << endl;
     123 
     124    // 1) energy smearing
     125    float energyS = -1.;
     126    if (particle->PID == pGAMMA)
     127        energyS = gRandom->Gaus(particle->E, sqrt( pow(DET->ELG_Nzdc,2) +
     128                                                   pow(DET->ELG_Czdc*particle->E,2) +
     129                                                   pow(DET->ELG_Szdc*sqrt(particle->E),2) ));
     130    else // smearing with hadronic resolution
     131        energyS =  gRandom->Gaus(particle->E, sqrt( pow(DET->HAD_Nzdc,2) +
     132                                                    pow(DET->HAD_Czdc*particle->E,2) +
     133                                                    pow(DET->HAD_Szdc*sqrt(particle->E),2) ));
     134    elementZdc->E = energyS;
     135 
     136
     137    // 2) time of flight t is t = T + d/[ cos(theta) v ]
     138    float cos_theta = 1;      //very good approximation, if eta_zdc >3
     139    if (DET->VFD_min_zdc<3) { // if smaller eta -> make the complete calculation
     140       double tx = atan(particle->Px/particle->Pz);
     141       double ty = atan(particle->Py/particle->Pz);
     142       double theta = sqrt( pow(tx,2) + pow(ty,2) );
     143       //cout << "tx = " << tx << " ty = " << ty << " theta  = " << theta << " cos(theta) = " << cos(theta) << endl;
     144       // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1
     145       // eta = 2.6 <-> cos(theta) = 0.99
     146       // eta = 3.0 <-> cos(theta) = 0.995
     147       cos_theta = cos(theta);
     148    }
     149    // units from StdHEP : Z [mm] T[mm/c]
     150    // units from Delphes : VFD_s_zdc [m] speed_of_light [m/s]
     151    double flight_distance = (DET->VFD_s_zdc - particle->Z*(1E-3))/cos_theta ;
     152    double flight_time = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s]
     153    double timeS = gRandom->Gaus(flight_time,DET->ZDC_T_resolution);
     154    elementZdc->T = timeS;
     155
     156    // 3) side: which ZDC has been hit?
    139157    elementZdc->side = sign(particle->Eta);
    140    
    141    
     158
     159    // 4) object nature : e.m. (photon) or had (neutron) ?
     160    elementZdc->hadronic_hit = (bool) (particle->PID==pN);
    142161  }
    143162 
    … …  
    148167 
    149168  TRootRomanPotHits* elementRP220;
    150   TRootRomanPotHits* elementFP420;
     169  TRootForwardTaggerHits* elementFP420;
    151170 
    152171  TLorentzVector genMomentum;
    … …  
    177196          elementRP220->Ty = (1E-6)*p1.getTY(); // [rad]
    178197          elementRP220->S = p1.getS();          // [m]
    179           // in first approximation only ! this number is always lower than the real distance-of-flight
    180           double flight_distance = p1.getS() - particle->Z*(1E-6);
    181           elementRP220->T = particle->T + flight_distance/speed_of_light; // assumes highly relativistic particles
    182 //cout << "T = " << particle->T << " ; " << flight_distance/speed_of_light << endl;
    183           elementRP220->E = p1.getE();          // not yet implemented
    184           elementRP220->q2 = -1;                // not yet implemented
    185           elementRP220->side = sign(particle->Eta);
     198
     199            /* time of flight t is t = T + d/[ cos(theta) v ]
     200            // nb: here we assume a straight path to the detector, which is not the case!
     201            // this time estimate is always underestimated (while exact for the ZDC case)
     202            float cos_theta = 1;      //very good approximation, if CEN_max_calo_fwd >3
     203            if (DET->CEN_max_calo_fwd<3) { // if smaller eta -> make the complete calculation
     204               double tx = atan(particle->Px/particle->Pz);
     205               double ty = atan(particle->Py/particle->Pz);
     206               double theta = sqrt( pow(tx,2) + pow(ty,2) );
     207               //cout << "tx = " << tx << " ty = " << ty << " theta  = " << theta << " cos(theta) = " << cos(theta) << endl;
     208               // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1
     209               // eta = 2.6 <-> cos(theta) = 0.99
     210               // eta = 3.0 <-> cos(theta) = 0.995
     211               cos_theta = cos(theta);
     212            }
     213            // units from StdHEP : Z [mm] T[mm/c]
     214            // units from Delphes : p1.getS [m]   speed_of_light [m/s]
     215            //double flight_distance = (p1.getS() - particle->Z*(1E-3))/cos_theta ;
     216            //elementRP220->T = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s]
     217            */
     218            elementRP220->E = p1.getE();        // not yet implemented
     219            elementRP220->q2 = -1;                // not yet implemented
     220            elementRP220->side = sign(particle->Eta);
    186221         
    187222        } else if (p1.getStoppingElement()->getName()=="rp420_1" || p1.getStoppingElement()->getName()=="rp420_2") {
    188223          p1.propagate(DET->RP_420_s);
    189           elementFP420 = (TRootRomanPotHits*) branchFP420->NewEntry();
     224          elementFP420 = (TRootForwardTaggerHits*) branchFP420->NewEntry();
    190225          elementFP420->X  = (1E-6)*p1.getX();  // [m]
    191226          elementFP420->Y  = (1E-6)*p1.getY();  // [m]
    … …  
    193228          elementFP420->Ty = (1E-6)*p1.getTY(); // [rad]
    194229          elementFP420->S = p1.getS();          // [m]
    195           // in first approximation only ! this number is always lower than the real distance-of-flight
    196           double flight_distance = p1.getS() - particle->Z*(1E-6);
    197 //cout << "T = " << particle->T << " ; " << flight_distance/speed_of_light << endl;
    198           elementFP420->T = particle->T + flight_distance/speed_of_light;
     230
     231            // time of flight t is t = T + d/[ cos(theta) v ]
     232            // nb: here we assume a straight path to the detector, which is not the case!
     233            // this time estimate is always underestimated (while exact for the ZDC case)
     234            float cos_theta = 1;      //very good approximation, if CEN_max_calo_fwd >3
     235            if (DET->CEN_max_calo_fwd<3) { // if smaller eta -> make the complete calculation
     236               double tx = atan(particle->Px/particle->Pz);
     237               double ty = atan(particle->Py/particle->Pz);
     238               double theta = sqrt( pow(tx,2) + pow(ty,2) );
     239               //cout << "tx = " << tx << " ty = " << ty << " theta  = " << theta << " cos(theta) = " << cos(theta) << endl;
     240               // NB: in practice, eta= 8 <-> theta 0.038° <-> 7x10^-4 rad <-> cos(theta) ~1
     241               // eta = 2.6 <-> cos(theta) = 0.99
     242               // eta = 3.0 <-> cos(theta) = 0.995
     243               cos_theta = cos(theta);
     244            }
     245            // units from StdHEP : Z [mm] T[mm/c]
     246            // units from Delphes : p1.getS [m]   speed_of_light [m/s]
     247            double flight_distance = (p1.getS() - particle->Z*(1E-3))/cos_theta ;
     248            elementFP420->T = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s]
    199249          elementFP420->E = p1.getE();          // not yet implemented
    200250          elementFP420->q2 = -1;                // not yet implemented
Note: See TracChangeset for help on using the changeset viewer.