VetoRecon.cxx Source File

00001 // $Header: /nfs/slac/g/glast/ground/cvs/reconstruction/src/VetoRecon.cxx,v 1.13 2001/03/21 14:38:21 heather Exp $
00002 
00003 
00004 #include "reconstruction/VetoRecon.h"
00005 
00006 #include "geometry/Ray.h"
00007 #include "geometry/Box.h"
00008 //#include "instrument/Scintillator.h"
00009 
00010 #if 0 //THB
00011 #include "glastmedia/ACDTopMed.h"
00012 #include "glastmedia/ACDSideMed.h"
00013 #endif
00014 
00015 #include "DetectorData.h"
00016 #include "idents/AcdId.h"
00017 
00018 #include "reconstruction/CalRecon.h"
00019 #include "reconstruction/ReconVisitor.h"
00020 #include "reconstruction/TriggerRecon.h"
00021 #include "reconstruction/TrackerRecon.h"
00022 
00023 #include "gui/DisplayRep.h"
00024 #include <algorithm>
00025 #include <cstdio>
00026 
00027 using namespace std;
00028 
00029 double VetoRecon::threshold_energy;
00030 double VetoRecon::tile_width_x;
00031 double VetoRecon::tile_width_y;
00032 const unsigned int VetoRecon::firstModId = 0;
00033 unsigned int VetoRecon::lastModId;
00034 
00035 
00036 static double scintillator_threshold_energy = 0.00069999997504056; 
00037 static int xNum = 4;
00038 static int yNum = 4;
00039 
00040 static double tower_mod_width = 38.970000000000;
00041 static double tower_wall_thickness  = 0.15000000000000;
00042 static double glast_wall_gap = 0.10000000000000;
00043 
00044 VetoRecon::VetoRecon (CalRecon* cal, TrackerRecon* trk)
00045 :Recon()
00046 ,cal(cal), track(trk), vetoData(0), m_DOCA(200)
00047 {
00048     getParameters();
00049     setTitle("Glast VetoRecon: ver 2.0");
00050 
00051     // Initiate Labelled Data Array...
00052     i_vto_Sys       = data.addElement("Veto_System");
00053     i_vto_TKR_E     = data.addElement("Veto_TKR_E");
00054     i_vto_Energy    = data.addElement("Veto_Energy");
00055     i_vto_Count     = data.addElement("No_Vetos_Hit");
00056     i_vto_vkeep     = data.addElement("Trig_vtoKeep");
00057     i_vto_TkrKeep  = data.addElement("Trig_tkrKeep");
00058     i_vto_DOCA      = data.addElement("Veto_DOCA");
00059     i_vto_DOCA_TOP  = data.addElement("Veto_DOCA_Top");
00060     i_vto_DOCA_S1   = data.addElement("Veto_DOCA_S1");
00061     i_vto_DOCA_S2   = data.addElement("Veto_DOCA_S2");
00062     i_vto_DOCA_S3   = data.addElement("Veto_DOCA_S3");
00063 }
00064 
00065 
00066 void VetoRecon::getParameters ()
00067 {
00068     threshold_energy = scintillator_threshold_energy;
00069     lastModId = (DetectorData::Glast::xNum * yNum) - 1;
00070 
00071     // Compute the tile widths - in the future modify ACDTopMed and ACDSideMed to
00072     // allow access to variables storing these values.
00073     unsigned int totY, totX;
00074     if( false /*THB ACDTopMed::align == 1 */) {
00075            totY = yNum;
00076            totX = xNum;
00077     } else {
00078            totX = xNum + 1;
00079            totY = yNum + 1;
00080     }
00081     double modWid = tower_mod_width + 2.*tower_wall_thickness + glast_wall_gap;
00082     tile_width_x = 31.496000000000; //((Glast::xNum*modWid) / totX ) - ACDTopMed::veto_gap;
00083     tile_width_y = 31.496000000000; //((Glast::yNum*modWid) / totY ) - ACDTopMed::veto_gap;
00084 }
00085 
00086 void VetoRecon::accept (ReconVisitor& rv)
00087 {
00088     rv.visit(this);
00089 }
00090 
00091 
00092 void VetoRecon::reconstruct (const IVetoData* v)
00093 {
00094     if(state() == DONE) return;
00095 
00096     vetoData = v;
00097 
00098     std::vector<double> TKR_E; // store the energy in tile(s) pierced by TKR recon
00099     std::vector<idents::AcdId> TKR_ID;  // store tile id too
00100 
00101     int count = 0;      // number of hit tiles
00102     double totEnergy = 0.0; 
00103 
00104     // Flags denoting whether Recon pierced lit tiles
00105     int calKeep = 1, calKeep2 = 1;  // init to "keep" this event
00106     int trackKeep = 1, trackKeep2 = 1; 
00107 
00108     // Nearest Neighbor - if Nearest Neighbor tiles next to a Tower with 3 trays in a row are lit
00109     //  then toss this event out.  Init to SAVE the event
00110     int NNflag = 1;  // Flag denoting whether or not NN criteria fails
00111 
00112     Box *tile = 0;
00113 
00114     const double maxDist = 200.0;  // between a tile and the ray representing a trajectory
00115 
00116     // Create Rays for the TKR and CAL recons
00117     bool haveCalRecon = true, haveTkrRecon = true;
00118     Vector dir = cal->direction();
00119     if ( (dir.x() == 0.) && (dir.y() == 0.) && (dir.z() == 0.) ) haveCalRecon = false;
00120     Ray calRay1(cal->center(), dir);
00121     Vector calOppDir(-(dir));
00122     Ray calRay2(cal->center(), calOppDir); 
00123 
00124     // Check to see if there is a valid TKR reconstruction, otherwise return
00125     // This is equivalent to No_Tracks in ntuple, must have at least 1 x and 1 y track
00126     if ( (track->xtrackList.size() <= 0) || (track->ytrackList.size() <= 0) ) haveTkrRecon = false;
00127     Ray tkrRay1(track->center(), track->direction());
00128     Vector tkrOppDir(-(track->direction()));
00129     Ray tkrRay2(track->center(), tkrOppDir);
00130 
00131     //THB: duplicate here from instrument.xml -- should fix this!
00132     double side_tile_heights[]= {25.0,20.0,15.0,15.0};
00133 
00134     // iterate through all hit ACD tiles
00135     for( IVetoData::const_iterator it= v->begin(); it != v->end(); ++it) {
00136         
00137         if ((*it).energy() < threshold_energy) continue; // toss out hits below threshold
00138 
00139         count++;
00140         double energy = (*it).energy();
00141         totEnergy += energy;
00142         idents::AcdId type = (*it).type();
00143 
00144         // Continue only if this is a top tile or upper side tile
00145         if ( (type.top() == false) && (type.row() != 0) ) continue;
00146 
00147         // Create a box to represent the current ACD tile
00148         if (type.top()) {   // top
00149             tile = new Box(tile_width_x, tile_width_y, 1.0 /*ACDTopMed::veto_thickness*/);
00150             tile->move(Vector((*it).position().x(), (*it).position().y(), (*it).position().z()));
00151         } else {            // side
00152             int index = 4 /*ACDSideMed::side_tile_heights.size()*/ - type.row() - 1;
00153             if ( (type.face() == 2) || (type.face() == 4) ) { // Face 2 or 4
00154                 tile = new Box(1.0 /*ACDTopMed::veto_thickness*/, tile_width_x, side_tile_heights[index]);
00155                 tile->rotateZ(M_PI * 0.5);
00156             } else {        // Face 1 or 3
00157                 tile = new Box(1.0 /*ACDTopMed::veto_thickness*/, tile_width_y, side_tile_heights[index]);
00158             }
00159             tile->move(Vector((*it).position().x(),(*it).position().y(),(*it).position().z()));
00160         }
00161 
00162 
00163         // Check to see if the CAL recon pierces this lit tile
00164         double dist;
00165         if (haveCalRecon) {
00166             dist = tile->distanceToEnter(calRay1, maxDist);
00167             if (dist < FLT_MAX) {
00168                 if (type.layer() == 0) calKeep = 0;
00169                 else calKeep2 = 0;
00170             }
00171             dist = tile->distanceToEnter(calRay2, maxDist);
00172             if (dist < FLT_MAX) {
00173                 if (type.layer() == 0) calKeep = 0;
00174                 else calKeep2 = 0;
00175             }
00176         }
00177 
00178         // Check to see if the TKR recon pierces this lit tile
00179         if (haveTkrRecon) {
00180             dist = tile->distanceToEnter(tkrRay1, maxDist);
00181             if (dist < FLT_MAX) {
00182                 if (type.layer() == 0) trackKeep = 0;
00183                 else trackKeep2 = 0;
00184                 TKR_E.push_back(energy);
00185                 TKR_ID.push_back(type);
00186             }
00187             dist = tile->distanceToEnter(tkrRay2, maxDist);
00188             if (dist < FLT_MAX) {
00189                 if (type.layer() == 0) trackKeep = 0;
00190                 else trackKeep2 = 0;
00191                 TKR_E.push_back(energy);
00192                 TKR_ID.push_back(type);
00193             }
00194         }
00195 
00196         delete tile;
00197 
00198         // LEVEL 2 TRIGGER (temporary calculation until readout 
00199         //    project is further along)
00200         // Determine if Nearest Neighbor Tiles are lit up, next to Tower w/ 3 trays in a row
00201         if(NNflag == 1) TriggeredTowers(type, &NNflag); 
00202     }
00203 
00204     // LEVEL 2 TRIGGER (temporary location until readout project
00205     //   is further along)
00206     // A value of 1 means that the reconstruction DIDN'T hit a lit tile
00207     // A value of 0 means that recon DID HIT a lit tile
00208     // 1s place means cal recon hit a lit tile on the 1st layer
00209     // 10s place means track recon hit a lit time on the 1st layer
00210     // 100s place means cal recon hit a lit tile on 2nd layer (if it exists)
00211     // 1000s place means track recon hit a lit tile on 2nd layer (if it exists)
00212     int vkeep = /*THB (ACDTopMed::two_layer_top || ACDSideMed::two_layer_side)*trackKeep2*1000 +
00213         (ACDTopMed::two_layer_top || ACDSideMed::two_layer_side)*calKeep2*100 + */
00214         trackKeep*10 + calKeep;
00215 
00216     double maxE = 0.0;
00217     std::vector<double>::iterator loc_TKR_E = max_element(TKR_E.begin(), TKR_E.end());
00218     if (TKR_E.size() > 0) maxE = *loc_TKR_E;
00219 
00220     int TKR_type = -1;
00221     // Find the type of tile "hit" by TKR recon
00222     for(unsigned int i = 0; i < TKR_E.size(); i++) {
00223         if (TKR_E[i] == *loc_TKR_E) TKR_type = TKR_ID[i].id();
00224     }
00225     
00226 
00227     // Report which Veto Tile "hit" by TKR recon and energy deposited in tile 
00228     //  (if 2 tiles are pierced, then choose the one with greater energy
00229     // Recall that TKR_type is just the type of the tile pierced by TKR recon, it may or
00230     //  may not actually have lit.  TKR_type == -1 when there is NO tile pierced by TKR recon
00231     //  TKR_E is the energy deposited in this tile, and will
00232     //  denote whether or not this tile was actually lit (above threshold)
00233     data.addData(i_vto_Sys, TKR_type);
00234     data.addData(i_vto_TKR_E, maxE);
00235     data.addData(i_vto_Energy, totEnergy);
00236     data.addData(i_vto_Count, count);
00237     data.addData(i_vto_vkeep, vkeep);
00238     data.addData(i_vto_TkrKeep, NNflag);
00239 
00240  // Section to determine distance to nearest hit veto counter
00241  // Loop over all x/y combo's of tracks..
00242     int numx = track->xtrackList.size();
00243     int numy = track->ytrackList.size();
00244     m_DOCA = 200.; // default value for no hit found case
00245 
00246     double m_DOCA_TOP = 200., m_DOCA_S1 = 200., 
00247         m_DOCA_S2 = 200., m_DOCA_S3 = 200.;
00248 
00249     for (int ix=0; ix<numx; ix++) {
00250         KalTrack xtrk = track->xtrackList[ix];
00251         for(int iy=0; iy<numy; iy++) {
00252             KalTrack ytrk = track->ytrackList[iy];
00253             Vector dir = Vector(xtrk.slope(), ytrk.slope(), 1.).unit();
00254             Point x0(xtrk.positionAtZ(0), ytrk.positionAtZ(0), 0.);
00255             float testDOCA = vetoDOCAFor(Ray(x0, dir), m_DOCA_TOP, m_DOCA_S1, 
00256             m_DOCA_S2, m_DOCA_S3);
00257             if(testDOCA < m_DOCA) m_DOCA = testDOCA;
00258         }
00259     }
00260     data.addData(i_vto_DOCA, m_DOCA);
00261 
00262     data[i_vto_DOCA_TOP] = m_DOCA_TOP;
00263     data[i_vto_DOCA_S1] = m_DOCA_S1;
00264     data[i_vto_DOCA_S2] = m_DOCA_S2;
00265     data[i_vto_DOCA_S3] = m_DOCA_S3;
00266 
00267     setState(DONE);
00268 
00269     return;
00270 }
00271 
00272 void VetoRecon::clear ()
00273 {
00274     Recon::clear();
00275     data.clear();
00276 }
00277 
00278 float VetoRecon::vetoDOCAFor (const Ray& t, double &top, double &s1, double &s2, double  &s3)
00279 {
00280     // This section looks for close-by hits in the veto systems
00281 
00282     float vetoDOCA = 200.;
00283     float dist;
00284 
00285     Point x0 = t.position(0.);
00286     Vector t0 = t.direction();
00287 
00288     for( IVetoData::const_iterator it= vetoData->begin(); it != vetoData->end(); ++it) {
00289 
00290         // Protect against low-energy backsplash
00291                 float veto_energy = (*it).energy();
00292                 if(veto_energy < threshold_energy) continue;
00293 
00294         Vector dX = (*it).position() - x0;
00295 
00296         float prod = dX * t0;
00297         dist = sqrt(dX.mag2() - prod*prod);
00298         if(dist < vetoDOCA){
00299             vetoDOCA = dist;
00300             m_hit_tile = *it;
00301         }
00302         
00303         idents::AcdId type = (*it).type();
00304 
00305         // Pick up the min. dist. from each type of tile
00306         // i.e. top, and each type of side row tile
00307         if (type.top() && dist < top) top = dist;
00308         if (type.side()) {
00309             if (type.row() == 0 && dist < s1) s1 = dist;
00310             if (type.row() == 1 && dist < s2) s2 = dist;
00311             if (type.row() == 2 && dist < s3) s3 = dist;
00312         }
00313     }
00314 
00315     return vetoDOCA;
00316 }
00317 
00318 void VetoRecon::draw (gui::DisplayRep& v)
00319 {
00320     if( m_DOCA>100. ) return;
00321     v.set_color("red");
00322     v.markerAt(m_hit_tile.position());
00323 #if 1 // disable to not plot the DOCA value
00324     v.move_to(m_hit_tile.position()); 
00325     char text[10]; sprintf(text, "%5g", m_DOCA);
00326     v.drawText(text);
00327 #endif
00328     v.set_color("black");
00329 }
00330 
00331 void VetoRecon::TriggeredTowers(IVetoData::Id tileId, int *tkrKeep) {
00332 #if 0
00333    // This routine compares a hit tile (and module ID) 
00334   // to the Towers that had 3 trays in a row for this event.
00335   // If the Tile and triggered Tower are neigbors, set trkKeep = 0 
00336   // (meaning this event would be rejected due to the lit veto tile).
00337     for(unsigned int TrigId = firstModId; TrigId <= lastModId; TrigId++) { // Loop through all possible Towers
00338       // Did this tower trigger, and is it a neighbor to the tower with lit veto?
00339         if( (TriggerRecon::ThreeInARow(TrigId)) && (tileId.association (TrigId) != 0) ) { 
00340             *tkrKeep = 0;
00341             return;
00342        }
00343     }
00344 
00345 #endif
00346 }