ACDthrottle.cxx Source File

00001 // $Header: /nfs/slac/g/glast/ground/cvs/trigger/src/ACDthrottle.cxx,v 1.3 2001/10/02 13:12:38 heather Exp $
00002 
00003 // Include files
00004 // Gaudi system includes
00005 #include "GaudiKernel/MsgStream.h"
00006 #include "GaudiKernel/AlgFactory.h"
00007 #include "GaudiKernel/IDataProviderSvc.h"
00008 #include "GaudiKernel/SmartDataPtr.h"
00009 #include "GaudiKernel/Algorithm.h"
00010 
00011 // ntuple interface
00012 #include "ntupleWriterSvc/INTupleWriterSvc.h"
00013 
00014 
00015 // TDS class declarations: input data, and McParticle tree
00016 #include "GlastEvent/data/TdGlastData.h"
00017 #include "GlastEvent/MonteCarlo/McVertex.h"
00018 #include "GlastEvent/MonteCarlo/McParticle.h"
00019 #include "data/IVetoData.h"
00020 // for access to instrument.ini
00021 #include "GlastSvc/GlastDetSvc/IGlastDetSvc.h"
00022 #include "xml/IFile.h"
00023 
00024 // access to the id classes
00025 #include "idents/AcdId.h"
00026 
00027 
00028 #include <cassert>
00029 
00030 // Define the class here instead of in a header file: not needed anywhere but here!
00031 //------------------------------------------------------------------------------
00039 class ACDthrottle : public Algorithm {
00040 public:
00041     ACDthrottle(const std::string& name, ISvcLocator* pSvcLocator);
00042     StatusCode initialize();
00043     StatusCode execute();
00044     StatusCode finalize();
00045     
00046 private: 
00048     void getParameters();
00050     StatusCode fillNtuple();
00052     float calcTotEnergy(const IVetoData* acdDigiData);
00054     bool threeInARow(const SiData* tkrDigiData, unsigned int towerId);
00056     unsigned long TKRtowerbitpattern(const SiData* tkrDigiData, unsigned int towerId);
00058     IGlastDetSvc*    m_detSvc; 
00060     xml::IFile * m_ini; 
00062     int m_numTowers, m_xNumTowers, m_yNumTowers, m_xNumTopTiles, m_yNumTopTiles, m_nplanes;
00063     double m_veto_threshold;
00064     double m_mod_width;
00065     
00066     // elements to store in the ntuple
00067     float m_nhitface0,m_nhitface1, m_nhitface2, m_nhitface3, m_nhitface4;
00068     float m_nhitsiderow0, m_nhitsiderow1, m_nhitsiderow2, m_nhitsiderow3;
00069     float m_vetoword; 
00070     
00071     unsigned int assoctile[16][4];
00072     unsigned int assoctileside[16][4];
00073     
00075     INTupleWriterSvc *m_ntupleWriteSvc;
00077     std::string m_tupleName;
00078 };
00079 //------------------------------------------------------------------------
00080 
00081 // necessary to define a Factory for this algorithm
00082 // expect that the xxx_load.cxx file contains a call     
00083 //     DLL_DECL_ALGORITHM( ACDthrottle );
00084 
00085 static const AlgFactory<ACDthrottle>  Factory;
00086 const IAlgFactory& ACDthrottleFactory = Factory;
00087 
00088 //------------------------------------------------------------------------
00090 ACDthrottle::ACDthrottle(const std::string& name, ISvcLocator* pSvcLocator)
00091 :Algorithm(name, pSvcLocator)
00092 , m_detSvc(0), m_ini(0)
00093 {
00094     // declare properties with setProperties calls
00095     declareProperty("tupleName",  m_tupleName="");
00096     
00097 }
00098 //------------------------------------------------------------------------
00100 StatusCode ACDthrottle::initialize(){
00101     StatusCode  sc = StatusCode::SUCCESS;
00102     MsgStream log(msgSvc(), name());
00103     log << MSG::INFO << "initialize" << endreq;
00104     
00105     // Use the Job options service to set the Algorithm's parameters
00106     setProperties();
00107     
00108     if( m_tupleName.empty()) {log << MSG::ERROR << "tupleName property not set!"<<endreq;
00109     return StatusCode::FAILURE;}
00110     // now try to find the GlastDevSvc service
00111     if (service("GlastDetSvc", m_detSvc).isFailure()){
00112         log << MSG::ERROR << "Couldn't find the GlastDetSvc!" << endreq;
00113         return StatusCode::FAILURE;
00114     }
00115     // get the ini file
00116     m_ini = const_cast<xml::IFile*>(m_detSvc->iniFile()); //OOPS!
00117     assert(4==m_ini->getInt("glast", "xNum"));  // simple check
00118     
00119     // Call our routine to extract parameters from the instrument.xml file
00120     getParameters();
00121     
00122     // setup our ntuple...
00123     //    setupNtuple();
00124     
00125     
00126     // set up the array of associated acd front tiles
00127     assoctile[0][0] = 0x0000;
00128     assoctile[0][1] = 0x0001;
00129     assoctile[0][2] = 0x0010;
00130     assoctile[0][3] = 0x0011;
00131     int itow, itil;
00132     for (itow=1;itow<=3;itow++){
00133         for (itil=0;itil<=3;itil++){
00134             assoctile[itow][itil]=assoctile[itow-1][itil]+0x0001;
00135         }
00136     }
00137     for (itow=4;itow<=15;itow++){
00138         for (itil=0;itil<=3;itil++){
00139             assoctile[itow][itil]=assoctile[itow-4][itil]+0x0010;
00140         }
00141     }
00142     //  for (itow=8;itow<=11;itow++){
00143     //          for (itil=0;itil<=3;itil++){
00144     //         assoctile[itow][itil]=assoctile[itow-4][itil]+0x0010;
00145     //          }
00146     //  }
00147     //  for (itow=12;itow<=15;itow++){
00148     //          for (itil=0;itil<=3;itil++){
00149     //        assoctile[itow][itil]=assoctile[itow-4][itil]+0x0010;
00150     //          }
00151     //  }
00152     // set up the array of associated acd side tiles
00153     for (itow=0;itow<=15;itow++){
00154         for (itil=0;itil<=3;itil++){
00155             assoctileside[itow][itil]=0x0999;
00156         }
00157     }
00158     assoctileside[0][0] = 0x0100;
00159     assoctileside[0][1] = 0x0101;
00160     assoctileside[0][2] = 0x0200;
00161     assoctileside[0][3] = 0x0201;
00162     assoctileside[1][0] = 0x0201;
00163     assoctileside[1][1] = 0x0202;
00164     assoctileside[2][0] = 0x0202;
00165     assoctileside[2][1] = 0x0203;
00166     assoctileside[3][0] = 0x0203;
00167     assoctileside[3][1] = 0x0204;
00168     assoctileside[3][2] = 0x0300;
00169     assoctileside[3][3] = 0x0301;
00170     assoctileside[7][0] = 0x0301;
00171     assoctileside[7][1] = 0x0302;
00172     assoctileside[11][0] = 0x0302;
00173     assoctileside[11][1] = 0x0303;
00174     assoctileside[15][0] = 0x0303;
00175     assoctileside[15][1] = 0x0304;
00176     assoctileside[15][2] = 0x0404;
00177     assoctileside[15][3] = 0x0403;
00178     assoctileside[14][0] = 0x0403;
00179     assoctileside[14][1] = 0x0402;
00180     assoctileside[13][0] = 0x0402;
00181     assoctileside[13][1] = 0x0401;
00182     assoctileside[12][0] = 0x0401;
00183     assoctileside[12][1] = 0x0400;
00184     assoctileside[12][2] = 0x0104;
00185     assoctileside[12][3] = 0x0103;
00186     assoctileside[8][0] = 0x0103;
00187     assoctileside[8][1] = 0x0102;
00188     assoctileside[4][0] = 0x0102;
00189     assoctileside[4][1] = 0x0101;
00190     //
00191     //    now let's whack a tile to study increase in rate
00192     //
00193     //    unsigned int iwhack = 0x0000;
00194     //  for (itow=0;itow<=15;itow++){
00195     //          for (itil=0;itil<=3;itil++){
00196     //          if (assoctile[itow][itil]== iwhack) assoctile[itow][itil]=0x0999;
00197     //          }
00198     //  }
00199     //
00200     // get a pointer to our ntupleWriterSvc
00201     if (service("ntupleWriterSvc", m_ntupleWriteSvc).isFailure()) {
00202         log << MSG::ERROR << "writeJunkAlg failed to get the ntupleWriterSvc" << endreq;
00203         return StatusCode::FAILURE;
00204     }
00205     
00206     //setupNtuple();
00207     return sc;
00208 }
00209 
00210 //------------------------------------------------------------------------
00212 StatusCode ACDthrottle::execute()
00213 {
00214     StatusCode  sc = StatusCode::SUCCESS;
00215     MsgStream   log( msgSvc(), name() );
00216     
00226     // Get the pointer to the event
00227     SmartDataPtr<TdGlastData> glastData(eventSvc(),"/Event/TdGlastData");
00228     
00229     // retrieve TKR data pointer for the event
00230     const SiData *tkrDigiData = glastData->getSiData();
00231     if (tkrDigiData == 0) {
00232         log << MSG::INFO << "No TKR Data available" << endreq;
00233     } else {
00234         // find a Tower with three in a row
00235         bool anyTower = false;
00236         unsigned int iTower;
00237         for (iTower = 0; iTower < m_numTowers-1; iTower++) {
00238             if (threeInARow(tkrDigiData, iTower) == true) {
00239                 anyTower = true;
00240                 break;
00241             }
00242         }
00243     }
00244     
00245     
00246     // retrieve CAL data pointer
00247     const CsIData *calDigiData = glastData->getCsIData();
00248     if (calDigiData == 0) {   
00249         log << MSG::INFO << "No CAL Data available" << endreq;
00250     } else {
00251         // check for CAL-Hi and CAL-Lo
00252         bool CALHI = false;
00253         bool CALLO = false;
00254         // skip this for now.  Wait until L1T properly implemented.  Don't need it for L2 studies yet.
00255     }
00256     
00257     // retrieve the ACD data pointer
00258     const IVetoData *acdDigiData = glastData->getVetoData();
00259     // Check to see if there is any ACD data for this event - if not provide a message
00260     if (acdDigiData == 0) log << MSG::INFO << "No ACD Data available" << endreq;
00261     
00262     // if we have ACD data, calculate the number of hit tiles by face and row
00263     // initialize counters number of tiles on each face
00264     int nhitface[5] = {0.,0.,0.,0.,0.};
00265     // initialize counters number of tiles in each row on the sides
00266     int nhitsiderow[5] = {0.,0.,0.,0.,0.};
00267     unsigned long bitpattern;
00268     if (acdDigiData != 0) {
00269         enum {
00270             _layermask = 0x0800,
00271                 _facemask  = 0x0700,
00272                 _rowmask   = 0x00F0,
00273                 _colmask   = 0x000F
00274         };
00275         for( IVetoData::const_iterator it= acdDigiData->begin(); it != acdDigiData->end(); ++it) {
00276             if ((*it).energy() > m_veto_threshold){
00277                 // add up tiles by face and by row on the side
00278                 unsigned int face = ((*it).type() & _facemask)>>8;
00279                 nhitface[face] += 1;
00280                 if (face != 0) {
00281                     unsigned int row = ((*it).type() & _rowmask)>>4;
00282                     nhitsiderow[row] += 1;
00283                 }
00284             }
00285         }
00286         
00287         // now let's examine those L1Ts that easily match ACD tiles
00288         //  get the bit pattern by tower
00289         unsigned int vetoword = 0;
00290         if (tkrDigiData == 0) {
00291             log << MSG::INFO << "No TKR Data available" << endreq;
00292         } else {
00293             unsigned int iTower;
00294             for (iTower = 0; iTower < m_numTowers-1; iTower++) {
00295                 bitpattern = TKRtowerbitpattern(tkrDigiData,iTower);
00296                 if (bitpattern > 0) {
00297                     if ((bitpattern & 7) == 7) {
00298                         // it's a top-layer conversion trigger
00299                         // loop over hit acd tiles and check if one is associated with this tower
00300                         // just the top guys for now                      
00301                         for( IVetoData::const_iterator it= acdDigiData->begin(); it != acdDigiData->end(); ++it) {
00302                             if ((*it).energy() > m_veto_threshold){
00303                                 unsigned int tileno = (*it).type();
00304                                 int itile;
00305                                 for (itile = 0;itile<=3;itile++){
00306                                     if (tileno == assoctile[iTower][itile]) vetoword |= 1;
00307                                     //                            if (assoctile[iTower][itile]==0x0999) {
00308                                     //                              m_guiSvc->guiMgr()->setState(gui::GuiMgr::PAUSED);
00309                                     //                            }
00310                                 }
00311                             }
00312                         }
00313                     } else {
00314                         // not a top conversion.  let's check the side tiles for this tower
00315                         unsigned int iPlane;
00316                         for (iPlane = 1; iPlane < m_nplanes-2; iPlane++) {
00317                             unsigned long mask = (7 << iPlane);
00318                             if ((bitpattern & mask) == mask) {
00319                                 // we found the first XY coincidence for the 3-in-a-row.  
00320                                 // Look at side ACD tiles in the appropriate row.
00321                                 for( IVetoData::const_iterator it= acdDigiData->begin(); it != acdDigiData->end(); ++it) {
00322                                     if ((*it).energy() > m_veto_threshold){
00323                                         unsigned int tileno = (*it).type();
00324                                         int itile;
00325                                         for (itile = 0;itile<=3;itile++){
00326                                             int rowassoc = 0;
00327                                             if (iPlane > 3) rowassoc = 1;
00328                                             if (iPlane > 10) rowassoc = 2;
00329                                             if (iPlane >14) rowassoc = 3;
00330                                             if (tileno == (assoctileside[iTower][itile]+rowassoc*0x0010)) {
00331                                                 vetoword |= 0x0010*(1<<rowassoc);  // record in vetoword which row in ACD vetos
00332                                                 iPlane = m_nplanes+1;  //break out of this big loop.  CHECK THIS.
00333                                                 break; //  break out of this inner loop
00334                                             }
00335                                             // check wider angles if nothing has been found yet
00336                                             rowassoc = 0;
00337                                             if (iPlane > 10) rowassoc = 1;
00338                                             if (iPlane > 14) rowassoc = 2;
00339                                             if (tileno == (assoctileside[iTower][itile]+rowassoc*0x0010)) {
00340                                                 vetoword |= 0x0010*(1<<rowassoc);  // record in vetoword which row in ACD vetos
00341                                                 iPlane = m_nplanes+1;  //break out of this big loop.  CHECK THIS.
00342                                                 break; //  break out of this inner loop
00343                                             }
00344                                         }
00345                                     }
00346                                     if (iPlane > m_nplanes) break; // check to see if we're done with loop
00347                                 }
00348                             } 
00349                         }
00350                     }
00351                 } 
00352             }
00353         }
00354         
00355         // At the end - we fill the tuple with the data from this event
00356         m_vetoword = vetoword;
00357         m_nhitface0 = nhitface[0];
00358         m_nhitface1 = nhitface[1];
00359         m_nhitface2 = nhitface[2];
00360         m_nhitface3 = nhitface[3];
00361         m_nhitface4 = nhitface[4];
00362         m_nhitsiderow0 = nhitsiderow[0];
00363         m_nhitsiderow1 = nhitsiderow[1];
00364         m_nhitsiderow2 = nhitsiderow[2];
00365         m_nhitsiderow3 = nhitsiderow[3];
00366         //      
00367     }
00368     
00369     sc = fillNtuple();
00370     return sc;
00371 }
00372 
00373 //------------------------------------------------------------------------
00375 StatusCode ACDthrottle::finalize(){
00376     StatusCode  sc = StatusCode::SUCCESS;
00377     MsgStream log(msgSvc(), name());
00378     
00379     return sc;
00380 }
00381 
00382 
00384 void ACDthrottle::getParameters ()
00385 {
00386     
00387     MsgStream   log( msgSvc(), name() );
00388     
00389     m_yNumTowers = m_ini->getInt("glast", "yNum");
00390     m_xNumTowers = m_ini->getInt("glast", "xNum");
00391     
00392     m_numTowers = (m_yNumTowers * m_xNumTowers);
00393     
00394     m_nplanes = m_ini->getInt("tracker", "num_trays");
00395     --m_nplanes;  // # of planes is one less than # of trays
00396     
00397     m_xNumTopTiles = m_ini->getInt("veto", "numXtiles");
00398     m_yNumTopTiles = m_ini->getInt("veto", "numYtiles");
00399     m_veto_threshold = m_ini->getDouble("veto", "threshold");
00400     
00401     m_mod_width = m_ini->getDouble("glast", "mod_width");
00402     
00403     log << MSG::INFO << "Retrieved data from the instrument.xml file" << endreq;
00404     log << MSG::INFO << "The number of towers in X = " << m_xNumTowers << endreq;
00405     
00406 }
00407 
00408 StatusCode ACDthrottle::fillNtuple() {
00409     StatusCode  sc = StatusCode::SUCCESS;
00410     if (m_tupleName.empty()) return sc;
00411     
00412            // Here we are adding to our ROOT ntuple
00413     
00414     // setup the entries to our tuple
00415     if ((sc = m_ntupleWriteSvc->addItem(m_tupleName.c_str(), "ACD_Throttle_Bits", m_vetoword)).isFailure()) return sc;
00416     if ((sc = m_ntupleWriteSvc->addItem(m_tupleName.c_str(), "ACD_No_Face0", m_nhitface0)).isFailure()) return sc;
00417     if ((sc = m_ntupleWriteSvc->addItem(m_tupleName.c_str(), "ACD_No_Face1", m_nhitface1)).isFailure()) return sc;
00418     if ((sc = m_ntupleWriteSvc->addItem(m_tupleName.c_str(), "ACD_No_Face2", m_nhitface2)).isFailure()) return sc;
00419     if ((sc = m_ntupleWriteSvc->addItem(m_tupleName.c_str(), "ACD_No_Face3", m_nhitface3)).isFailure()) return sc;
00420     if ((sc = m_ntupleWriteSvc->addItem(m_tupleName.c_str(), "ACD_No_Face4", m_nhitface4)).isFailure()) return sc;
00421     if ((sc = m_ntupleWriteSvc->addItem(m_tupleName.c_str(), "ACD_No_SideRow0", m_nhitsiderow0)).isFailure()) return sc;
00422     if ((sc = m_ntupleWriteSvc->addItem(m_tupleName.c_str(), "ACD_No_SideRow1", m_nhitsiderow1)).isFailure()) return sc;
00423     if ((sc = m_ntupleWriteSvc->addItem(m_tupleName.c_str(), "ACD_No_SideRow2", m_nhitsiderow2)).isFailure()) return sc;
00424     if ((sc = m_ntupleWriteSvc->addItem(m_tupleName.c_str(), "ACD_No_SideRow3", m_nhitsiderow3)).isFailure()) return sc;
00425     
00426     return sc;
00427     
00428 }
00429 
00430 
00431 // just a little routine to demonstrate manipulating the ACD digi data.
00432 float ACDthrottle::calcTotEnergy(const IVetoData *acdDigiData) 
00433 {
00434     double totEnergy = 0.0;
00435     for( IVetoData::const_iterator it= acdDigiData->begin(); it != acdDigiData->end(); ++it) {
00436         totEnergy += (*it).energy();
00437     }
00438     return totEnergy;
00439 }
00440 
00441 // check to see if there is three In a Row in a particular Tower
00442 bool ACDthrottle::threeInARow(const SiData *tkrDigiData, unsigned int towerId) 
00443 {
00444     unsigned int iPlane;
00445     unsigned long bitString = TKRtowerbitpattern(tkrDigiData,towerId);
00446     for (iPlane = 0; iPlane < m_nplanes; iPlane++) {
00447         unsigned long mask = (7 << iPlane);
00448         if ((bitString & mask) == mask) return true;
00449     }
00450     return false;
00451 }
00452 // produce the bit pattern of XY coincidences by plane
00453 unsigned long ACDthrottle::TKRtowerbitpattern(const SiData *tkrDigiData, unsigned int towerId) 
00454 {
00455     unsigned long bitString = 0;
00456     unsigned int iPlane;
00457     for (iPlane = 0; iPlane < m_nplanes; iPlane++) {
00458         int nx = tkrDigiData->nHits(SiData::X, iPlane);
00459         int ny = tkrDigiData->nHits(SiData::Y, iPlane);
00460         
00461         if ( (nx > 0) && (ny > 0) ) {  // we have hits in both x & y
00462             
00463             bool foundXinTower = false;
00464             unsigned iXhit;
00465             for (iXhit = 0; iXhit < nx; iXhit++) {
00466                 if (towerId == tkrDigiData->moduleId(SiData::X, iPlane, iXhit)) {
00467                     foundXinTower = true;
00468                     break;
00469                 }
00470             }
00471             
00472             bool foundYinTower = false;
00473             unsigned iYhit;
00474             for (iYhit = 0; iYhit < ny; iYhit++) {
00475                 if (towerId == tkrDigiData->moduleId(SiData::Y, iPlane, iYhit)) {
00476                     foundYinTower = true;
00477                     break;
00478                 }
00479             }
00480             
00481             if (foundXinTower && foundYinTower) bitString |= (1 << iPlane);
00482             
00483         } else {  // do not have hits in both x & y..continue to the next plane
00484             continue;
00485         }
00486     }
00487     return bitString;
00488 }