Changeset 100 in svn for trunk/src

Timestamp:

Dec 18, 2008, 2:39:26 PM (16 years ago)

Author:

severine ovyn

Message:

Remove datacard bug + CaloTowers OK

Location:

trunk/src

Files:

• BFieldProp.cc (modified) (5 diffs)
• FrogUtil.cc (modified) (5 diffs)
• JetUtils.cc (modified) (8 diffs)
• VeryForward.cc (modified) (6 diffs)

trunk/src/BFieldProp.cc

@@ -25 +25 @@
 //------------------------------------------------------------------------------
 
-TrackPropagation::TrackPropagation() {
+TrackPropagation::TrackPropagation(string DetDatacard) {
 
+ DET = new RESOLution();
+ DET->ReadDataCard(DetDatacard);
  MAXITERATION = 10000;
  MINSEGLENGTH = 70;
@@ -43 +45 @@
   
   //out of trackibg coverage?
-  if(sqrt(Xvertex1*Xvertex1+Yvertex1*Yvertex1) > TRACK_radius){return;}
-  if(fabs(Zvertex1) > TRACK_length){return;}
+  if(sqrt(Xvertex1*Xvertex1+Yvertex1*Yvertex1) > DET->TRACK_radius){return;}
+  if(fabs(Zvertex1) > DET->TRACK_length){return;}
   
   float Px = Part->Px;
@@ -62 +64 @@
      double vz = pz/M;
 
-     double Bx = TRACK_bfield_x;
-     double By = TRACK_bfield_y;
-     double Bz = TRACK_bfield_z;
+     double Bx = DET->TRACK_bfield_x;
+     double By = DET->TRACK_bfield_y;
+     double Bz = DET->TRACK_bfield_z;
 
      double ax =  (q/M)*(Bz*vy - By*vz);
@@ -80 +82 @@
      int k = 0;
      
+     double Radius=DET->TRACK_radius;
+     double Length=DET->TRACK_length;
+
      while(k < MAXITERATION){
         k++;
@@ -99 +104 @@
         z  += vz*dt;
 
-       if( (x*x+y*y) > TRACK_radius*TRACK_radius ){ x /= (x*x+y*y)/(TRACK_radius*TRACK_radius); y /= (x*x+y*y)/(TRACK_radius*TRACK_radius); break;}
-       if( fabs(z)>TRACK_length)break;
+       if( (x*x+y*y) > Radius*Radius ){ x /= (x*x+y*y)/(Radius*Radius); y /= (x*x+y*y)/(Radius*Radius); break;}
+       if( fabs(z)>Length)break;
 
        xold = x;

trunk/src/FrogUtil.cc

@@ -485 +485 @@
   
   double Rayon_Tracker=50;
-  double Rayon_Calo = Rayon_Tracker*2;
-  double Lenght_Tracker=120;
-  double Lenght_Calo=Lenght_Tracker+Lenght_Tracker/2;
-  double Lenght_FWCalo=Lenght_Tracker+Lenght_Tracker;
-  
-  int NumPhi=10;
+  double Rayon_Calo = Rayon_Tracker*1.5;
+  double Rayon_Muon = Rayon_Tracker*2;
+  double Lenght_Tracker=100;
+  double Lenght_Calo=Lenght_Tracker+Lenght_Tracker/2.5;
+  double Lenght_Muon=Lenght_Calo+Lenght_Calo/2.5;
+  double Lenght_FWCalo=Lenght_Tracker+1.5*Lenght_Tracker;
+  
+  int plus=1;
+  int NumPhi=4;
   
   //************************************************Tracker*************************************************
@@ -498 +501 @@
   detector->addDaughter(Tracker);
   unsigned int DetIdCountTracker = 1;
-  for(double ray=0;ray <= Rayon_Tracker;ray +=5){
+  for(double ray=0;ray <= Rayon_Tracker;ray +=plus){
+  //  double ray=Rayon_Tracker;
     double length = ray/tan(EtaToTheta(DET->CEN_max_tracker));
     if(length >= Lenght_Tracker)
       {
-        FROG_Element_Primitive_Cylinder* tracker = new FROG_Element_Primitive_Cylinder(9100000+DetIdCountTracker*10,ray,0,0,0,0,0,Lenght_Tracker,NumPhi,0);
+        FROG_Element_Primitive_Cylinder* tracker = new FROG_Element_Primitive_Cylinder(9100000+DetIdCountTracker*10,ray,0,0,0,0,0,Lenght_Tracker,NumPhi,1);
+//      FROG_Element_Primitive_Cone* tracker = new FROG_Element_Primitive_Cone(9100000+DetIdCountTracker*10,ray,0,0,0,0,0,Lenght_Tracker,NumPhi,1);
         Tracker->addDaughter(tracker);   DetIdCountTracker++;
       }
@@ -518 +523 @@
   unsigned int DetIdCountCalo = 1;
   
-  for(double ray=Rayon_Tracker;ray <= Rayon_Calo ;ray +=5){
+  for(double ray=Rayon_Tracker;ray <= Rayon_Calo ;ray +=plus){
     double length = ray/tan(EtaToTheta(DET->CEN_max_calo_cen));
     double add;
@@ -529 +534 @@
       }
   }
-  for(double ray=0;ray <= Rayon_Tracker;ray +=5){
+  for(double ray=0;ray <= Rayon_Tracker;ray +=plus){
     double length = ray/tan(EtaToTheta(DET->CEN_max_calo_cen));
     double add;
@@ -552 +557 @@
   double eta_Max=DET->CEN_max_calo_fwd,eta_Min=DET->CEN_max_calo_cen;
   
-  for(double ray=0;ray <= 50;ray +=5){
+  for(double ray=0;ray <= 50;ray +=plus){
     double length = ray/tan(EtaToTheta(eta_Max));
     double add;
-    if(length >= Lenght_FWCalo)add=(Lenght_FWCalo-Lenght_Calo)/2;
-    else add=(length-Lenght_Calo)/2;
+    if(length >= Lenght_FWCalo)add=(Lenght_FWCalo-Lenght_Muon)/2;
+    else add=(length-Lenght_Muon)/2;
     
     if(add > 0 && ray < Lenght_FWCalo*tan(EtaToTheta(eta_Min)))
       {
-        FROG_Element_Primitive_Cylinder* caloFW1 = new FROG_Element_Primitive_Cylinder(9300000+DetIdCountFwCalo*10,ray,0,0,(Lenght_Calo+add)/2,0,0,add,NumPhi,0);
+        FROG_Element_Primitive_Cylinder* caloFW1 = new FROG_Element_Primitive_Cylinder(9300000+DetIdCountFwCalo*10,ray,0,0,(Lenght_Muon+add)/2,0,0,add,NumPhi,0);
         FWCALO->addDaughter(caloFW1);   DetIdCountFwCalo++;
-        FROG_Element_Primitive_Cylinder* caloFW2 = new FROG_Element_Primitive_Cylinder(9300000+DetIdCountFwCalo*10,ray,0,0,-(Lenght_Calo+add)/2,0,0,add,NumPhi,0);
+        FROG_Element_Primitive_Cylinder* caloFW2 = new FROG_Element_Primitive_Cylinder(9300000+DetIdCountFwCalo*10,ray,0,0,-(Lenght_Muon+add)/2,0,0,add,NumPhi,0);
         FWCALO->addDaughter(caloFW2);   DetIdCountFwCalo++;
       }
   }
+
+  //***********************************************Muon chambers********************************************
+  //********************************************************************************************************
+  //
+  FROG_Element_Base_With_DetId_And_Name* MUON = new FROG_Element_Base_With_DetId_And_Name(940000000,"MuonDet");
+  detector->addDaughter(MUON);
+  unsigned int DetIdCountMuon = 1;
+
+  for(double ray=Rayon_Calo;ray <= Rayon_Muon ;ray +=plus){
+    double length = ray/tan(EtaToTheta(DET->CEN_max_mu));
+    double add;
+    if(length >= Lenght_Muon)add=Lenght_Muon;
+    else add=ray/tan(EtaToTheta(DET->CEN_max_calo_cen));
+    if(add > 0)
+      {
+        FROG_Element_Primitive_Cylinder* muon = new FROG_Element_Primitive_Cylinder(9400000+DetIdCountMuon*10,ray,0,0,0,0,0,add,NumPhi,0);
+        MUON->addDaughter(muon);   DetIdCountMuon++;
+      }
+  }
+  for(double ray=0;ray <= Rayon_Calo;ray +=plus){
+    double length = ray/tan(EtaToTheta(DET->CEN_max_calo_cen));
+    double add;
+    if(length >= Lenght_Muon)add=(Lenght_Muon-Lenght_Calo)/2;
+    else add=(length-Lenght_Calo)/2;
+    if(add > 0)
+      {
+        FROG_Element_Primitive_Cylinder* muonP = new FROG_Element_Primitive_Cylinder(9400000+DetIdCountMuon*10,ray,0,0,(Lenght_Calo+add)/2,0,0,add,NumPhi,0);
+        MUON->addDaughter(muonP);   DetIdCountMuon++;
+        FROG_Element_Primitive_Cylinder* muonM = new FROG_Element_Primitive_Cylinder(9400000+DetIdCountMuon*10,ray,0,0,-(Lenght_Calo+add)/2,0,0,add,NumPhi,0);
+        MUON->addDaughter(muonM);   DetIdCountMuon++;
+      }
+  }
+
+  
     
   FROG_Geometry* CustomGeom = new FROG_Geometry(prim);

trunk/src/JetUtils.cc

@@ -15 +15 @@
 using namespace std;
 
-JetsUtil::JetsUtil() {
-  
-  switch(JET_jetalgo) {
+JetsUtil::JetsUtil(const string DetDatacard) {
+ 
+ DET = new RESOLution();
+ DET->ReadDataCard(DetDatacard);
+ 
+  switch(DET->JET_jetalgo) {
   default:
   case 1: {
-    
     // set up a CDF midpoint jet definition
 #ifdef ENABLE_PLUGIN_CDFCONES
-    plugins = new fastjet::CDFJetCluPlugin(JET_seed,JET_coneradius,JET_C_adjacencycut,JET_C_maxiterations,JET_C_iratch,JET_overlap);
+    plugins = new fastjet::CDFJetCluPlugin(DET->JET_seed,DET->JET_coneradius,DET->JET_C_adjacencycut,DET->JET_C_maxiterations,DET->JET_C_iratch,DET->JET_overlap);
     jet_def = fastjet::JetDefinition(plugins);
 #else
@@ -34 +36 @@
     // set up a CDF midpoint jet definition
 #ifdef ENABLE_PLUGIN_CDFCONES
-    plugins = new fastjet::CDFMidPointPlugin (JET_seed,JET_coneradius,JET_M_coneareafraction,JET_M_maxpairsize,JET_M_maxiterations,JET_overlap);
+    plugins = new fastjet::CDFMidPointPlugin (DET->JET_seed,DET->JET_coneradius,DET->JET_M_coneareafraction,DET->JET_M_maxpairsize,DET->JET_M_maxiterations,DET->JET_overlap);
     jet_def = fastjet::JetDefinition(plugins);
 #else
@@ -45 +47 @@
     // set up a siscone jet definition
 #ifdef ENABLE_PLUGIN_SISCONE
-    plugins = new fastjet::SISConePlugin (JET_coneradius,JET_overlap,JET_S_npass, JET_S_protojet_ptmin);
+    plugins = new fastjet::SISConePlugin (DET->JET_coneradius,DET->JET_overlap,DET->JET_S_npass, DET->JET_S_protojet_ptmin);
     jet_def = fastjet::JetDefinition(plugins);
 #else
@@ -54 +56 @@
     
   case 4: {
-    jet_def = fastjet::JetDefinition(fastjet::kt_algorithm, JET_coneradius);
+
+    jet_def = fastjet::JetDefinition(fastjet::kt_algorithm, DET->JET_coneradius);
   }
     break;
     
   case 5: {
-    jet_def = fastjet::JetDefinition(fastjet::cambridge_algorithm,JET_coneradius);
+    jet_def = fastjet::JetDefinition(fastjet::cambridge_algorithm,DET->JET_coneradius);
   }
     break;
     
   case 6: {
-    jet_def = fastjet::JetDefinition(fastjet::antikt_algorithm,JET_coneradius);
+    jet_def = fastjet::JetDefinition(fastjet::antikt_algorithm,DET->JET_coneradius);
   }
     break;
@@ -71 +74 @@
 }
 
-vector<fastjet::PseudoJet> JetsUtil::RunJets(const vector<fastjet::PseudoJet> & input_particles)
+vector<fastjet::PseudoJet> JetsUtil::RunJets(const vector<fastjet::PseudoJet>&  input_particles)
 {
   inclusive_jets.clear();
   sorted_jets.clear();
-  
   // run the jet clustering with the above jet definition
   if(input_particles.size()!=0)
@@ -96 +98 @@
   for (unsigned int i = 0; i < sorted_jets.size(); i++) {
     JET.SetPxPyPzE(sorted_jets[i].px(),sorted_jets[i].py(),sorted_jets[i].pz(),sorted_jets[i].E());
-    if(JET.Pt() > PTCUT_jet)
+    if(JET.Pt() > DET->PTCUT_jet)
       {
         elementJet = (TRootJet*) branchJet->NewEntry();
@@ -102 +104 @@
         // b-jets
         bool btag=false;
-        if((fabs(JET.Eta()) < CEN_max_tracker && Btaggedjet(JET, NFCentralQ)))btag=true;
+        if((fabs(JET.Eta()) < DET->CEN_max_tracker && DET->Btaggedjet(JET, NFCentralQ)))btag=true;
         elementJet->Btag = btag;
       }
@@ -116 +118 @@
     JET.SetPxPyPzE(sorted_jets[i].px(),sorted_jets[i].py(),sorted_jets[i].pz(),sorted_jets[i].E());
     // Tau jet identification : 1! track and electromagnetic collimation
-    if(fabs(JET.Eta()) < (CEN_max_tracker - TAU_track_scone)) {
-      double Energie_tau_central = EnergySmallCone(towers,JET.Eta(),JET.Phi());
+    if(fabs(JET.Eta()) < (DET->CEN_max_tracker - DET->TAU_track_scone)) {
+      double Energie_tau_central = DET->EnergySmallCone(towers,JET.Eta(),JET.Phi());
       if(
-         ( Energie_tau_central/JET.E() > TAU_energy_frac ) &&
-         ( NumTracks(TrackCentral,TAU_track_pt,JET.Eta(),JET.Phi()) == 1 ) &&
-         ( JET.Pt() > PTCUT_taujet)
+         ( Energie_tau_central/JET.E() > DET->TAU_energy_frac ) &&
+         ( DET->NumTracks(TrackCentral,DET->TAU_track_pt,JET.Eta(),JET.Phi()) == 1 ) &&
+         ( JET.Pt() > DET->PTCUT_taujet)
          ) {
         elementTauJet = (TRootTauJet*) branchTauJet->NewEntry();

trunk/src/VeryForward.cc

@@ -26 +26 @@
 //------------------------------------------------------------------------------
 
-VeryForward::VeryForward() {
-  
+VeryForward::VeryForward(string DetDatacard) {
+ 
+ DET = new RESOLution();
+ DET->ReadDataCard(DetDatacard);
+ 
   //Initialisation of Hector
   relative_energy = true; // should always be true
@@ -62 +65 @@
   TLorentzVector genMomentum; 
   // Zero degree calorimeter, for forward neutrons and photons
-  if (particle->Status ==1 && (pid == pN || pid == pGAMMA ) && eta > VFD_min_zdc ) {
+  if (particle->Status ==1 && (pid == pN || pid == pGAMMA ) && eta > DET->VFD_min_zdc ) {
     genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
     // !!!!!!!!! vérifier que particle->Z est bien en micromÚtres!!!
@@ -75 +78 @@
     //double theta = (1E-6)*sqrt( pow(tx,2) + pow(ty,2) );
     //double flight_distance = (DET->ZDC_S - particle->Z*(1E-6))/cos(theta) ; // assumes that Z is in micrometers
-    double flight_distance = (VFD_s_zdc - particle->Z*(1E-6));
+    double flight_distance = (DET->VFD_s_zdc - particle->Z*(1E-6));
     // assumes also that the emission angle is so small that 1/(cos theta) = 1
     elementZdc->T = 0*particle->T + flight_distance/speed_of_light; // assumes highly relativistic particles
@@ -95 +98 @@
   genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
   // if forward proton
-  if( (pid == pP)  && (particle->Status == 1) &&  (fabs(genMomentum.Eta()) > CEN_max_calo_fwd) )
+  if( (pid == pP)  && (particle->Status == 1) &&  (fabs(genMomentum.Eta()) > DET->CEN_max_calo_fwd) )
     {
       // !!!!!!!! put here particle->CHARGE and particle->MASS
@@ -112 +115 @@
       if(p1.stopped(beamline)) {
         if (p1.getStoppingElement()->getName()=="rp220_1" || p1.getStoppingElement()->getName()=="rp220_2") {
-          p1.propagate(RP_220_s);
+          p1.propagate(DET->RP_220_s);
           elementRP220 = (TRootRomanPotHits*) branchRP220->NewEntry();
           elementRP220->X  = (1E-6)*p1.getX();  // [m]
@@ -125 +128 @@
           
         } else if (p1.getStoppingElement()->getName()=="rp420_1" || p1.getStoppingElement()->getName()=="rp420_2") {
-          p1.propagate(RP_420_s);
+          p1.propagate(DET->RP_420_s);
           elementFP420 = (TRootRomanPotHits*) branchFP420->NewEntry();
           elementFP420->X  = (1E-6)*p1.getX();  // [m]