Micromegas LAPP reconstruction

00001 /*  @version $Revision: 1789 $ * @modifiedby $Author: lcdet $ * @lastmodified $Date: 2012-06-20 16:45:37 +0200 (Wed, 20 Jun 2012) $ */
00002 
00003 
00004 //#include <c++/4.0.0/bits/stl_map.h>
00005 
00006 
00007 
00008 #include "event/ChannelHit.hh"
00009 #include "event/Event.hh"
00010 #include "event/Run.hh"
00011 #include "geometry/Chip.hh"
00012 #include "geometry/Dif.hh"
00013 #include "geometry/Board.hh"
00014 #include "geometry/Chamber.hh"
00015 #include "geometry/Detector.hh"
00016 #include "geometry/Channel.hh"
00017 
00018 #include "tools/MicroException.hh"
00019 #include "tools/Log.hh"
00020 #include "slowControl/SlowControlEntry.hh"
00021 #include "slowControl/ChamberSlowControl.hh"
00022 
00023 
00024 #include <string>
00025 /*
00026 #include "UTIL/LCTime.h"
00027 #include "IMPL/LCFlagImpl.h"
00028 #include "IMPL/LCCollectionVec.h"
00029 #include "UTIL/CellIDEncoder.h"
00030 #include "IMPL/RawCalorimeterHitImpl.h"
00031 #include "IMPL/CalorimeterHitImpl.h"
00032 #include "IMPL/LCTOOLS.h"
00033 
00034 #include "IO/LCWriter.h"
00035 #include "EVENT/LCIO.h"
00036 #include "DATA/LCFloatVec.h"
00037 #include "DATA/LCIntVec.h"
00038 
00039 #include "IMPL/LCEventImpl.h"
00040 #include "IMPL/LCRunHeaderImpl.h"
00041 #include "IMPL/LCCollectionVec.h"
00042 #include "IMPL/RawCalorimeterHitImpl.h"
00043 #include "IMPL/SimTrackerHitImpl.h"
00044 #include "IMPL/MCParticleImpl.h"
00045 #include "IMPL/LCFlagImpl.h"
00046 #include "IMPL/LCTOOLS.h"
00047 
00048 #include "IMPL/TrackerRawDataImpl.h"
00049 #include "IMPL/TrackerDataImpl.h"
00050 #include "IMPL/TrackerPulseImpl.h"
00051 
00052 #include "UTIL/LCRelationNavigator.h"
00053 #include "UTIL/LCTime.h"
00054 #include "UTIL/CellIDEncoder.h"
00055 #include "UTIL/LCTypedVector.h"
00056 
00057 #include "UTIL/LCSplitWriter.h"
00058 
00059 static const int NEVENT = 10 ; // events
00060 static const int NMCPART = 10 ;  // mc particles per event
00061 static const int NHITS = 5 ;  // calorimeter hits per event
00062 */
00063 
00064 
00065 
00066 
00067 
00068 
00069 
00070 
00071 using namespace std;
00072 
00073 
00074   std::ostream& operator <<(std::ostream &out, const Event &x) {
00075     return(x.operator <<(out));
00076   }
00077 
00078 
00079 Event::Event(Run& aRun):run(aRun), slowControlEntry(NULL)
00080 {
00081   //FILE_LOG(logDEBUG) << "-Event contructor"<< endl;
00082 //  initChannelHitVector();
00083 timestamp = 0;
00084 difSynchro = 0;
00085 id = 0;
00086 slowControlEntry = NULL;
00087 crcIsCorrect = true;
00088 valid = VALID;
00089 }
00090 
00091  Event::Event(Run& aRun, i64 aTimestamp, i32 aId):run(aRun), timestamp(aTimestamp), id(aId), slowControlEntry(NULL)
00092  {
00093   timestamp = 0;
00094   difSynchro = 0;
00095   id = 0;
00096   slowControlEntry = NULL;
00097   crcIsCorrect = true;
00098   valid = VALID;
00099      //FILE_LOG(logDEBUG) << "-Event contructor with id: "<< aId << endl;
00100  }
00101 
00102 Event::~Event()
00103 {
00104   //FILE_LOG(logDEBUG) << "-Event destructor"<< endl;
00105 
00106    map<ui16, multiset<ChannelHit*,sorter>* >::iterator iter;
00107 
00108   for (iter=channelHits.begin(); iter!=channelHits.end(); iter++)
00109   {
00110       multiset<ChannelHit*,sorter>* list = (*iter).second;
00111       for (multiset<ChannelHit*,sorter>::iterator iterList = list->begin(); iterList!=list->end(); iterList++ )
00112       {
00113           delete (*iterList);
00114       }
00115       delete list;
00116   }
00117 
00118 }
00119 
00120 /*---------------------------------------------------*/
00121 ostream& Event::operator <<(ostream &out) const {
00122   return(out << "-Event id["<< id << "]" << endl
00123              << "-Event timestamp[" << timestamp << "]"  << endl);
00124 }
00125 
00126 /*---------------------------------------------------*/
00127 void Event::initChannelHitVector()
00128     {
00129       ChamberMap_t chambers = this->run.getDetector().getChambers();
00130 
00131       for (ChamberMap_t::iterator iter = chambers.begin() ; iter != chambers.end() ; iter++)
00132       {
00133         i32 chamberId  = (*iter).first;
00134         multiset<ChannelHit*, sorter> *list = new multiset<ChannelHit*, sorter>();
00135              FILE_LOG(logINFO) << "INIT map size BEFORE INSERT:[" << channelHits.size() << "]" << endl;
00136         channelHits[chamberId] = list;  // a ChannilHits* vector per Chamber
00137              FILE_LOG(logINFO) << "INIT map size AFTER INSERT:[" << channelHits.size() << "]" << endl;
00138       }
00139     }
00140 /*---------------------------------------------------*/
00141     const multiset<ChannelHit*,sorter>& Event::getChannelHitVector(const Chamber& aChamber) const
00142     {
00143 
00144          map<ui16, multiset<ChannelHit*,sorter>* >::const_iterator iter = channelHits.find(aChamber.getId());
00145          if( iter != channelHits.end() ) {
00146              return (*iter->second);
00147          }
00148 
00149          throw ( MicroException("No channelHit vector for this chamber"));
00150     //    return channelHits.find(aChamber.getId())->second;
00151     }
00152 /*------------------------------------------------------------------------------------*/
00153     void Event::insertHit(Chamber& aChamber, ChannelHit *aChannelHit)
00154     {
00155          ui16 chamberId  = aChamber.getId();
00156 
00157          multiset<ChannelHit*,sorter>* hits = NULL;
00158 
00159          map<ui16, multiset<ChannelHit*,sorter> *>::iterator iter = channelHits.find(chamberId);
00160 
00161          if( iter == channelHits.end() ) {
00162              hits = new  multiset<ChannelHit*, sorter>();
00163              channelHits.insert(pair<i32, multiset<ChannelHit*,sorter>* > (chamberId,hits));
00164          }
00165 
00166          else
00167          {
00168             hits = iter->second;
00169          }
00170 
00171          hits->insert(aChannelHit);
00172 
00173 
00174     }
00175 /*------------------------------------------------------------------------------------*/
00176 
00177  const ChannelHit& Event::getChannelHitByOrder( ui64 order ,const Chamber& aChamber, i16 aThreshold)
00178  {
00179      if ( getChannelHitByOrder(order+1,aChamber).getAnalogValue() <= aThreshold )
00180      {
00181          return getChannelHitByOrder(order,aChamber);
00182      }
00183      else
00184      {
00185          throw MicroException("No Channel respond for this event, this order, this threshold and this chamber ");
00186      }
00187  }
00188 
00189 /*------------------------------------------------------------------------------------*/
00190  const ChannelHit& Event::getChannelHitByOrder(ui64 order ,const Chamber& aChamber)
00191  {
00192         multiset<ChannelHit*,sorter>* chamberChannelHit = this->channelHits.find(aChamber.getId())->second;//channelHits[aChamber.getId()];
00193         ui64 items=0;
00194                 // Find first chamberChannelHit for select chamber
00195         for (multiset<ChannelHit*,sorter>::const_iterator iter=chamberChannelHit->begin(); iter!=chamberChannelHit->end(); iter++)
00196         {
00197                 if ( (*iter)->getChannel().getChamber().getId() == aChamber.getId() )
00198                 {
00199                     if ( order == items)
00200                     {
00201                         return *(*iter);
00202                     }
00203 
00204                     items++;
00205                 }
00206 
00207         }
00208   throw MicroException("No Channel respond or this event, this order and this chamber ");
00209 
00210  }
00211  /*
00212  //------------------------------------------------------------------------------------
00213  ui64 Event::getChannelHitOrder(const ChannelHit& aChannelHit ,const Chamber& aChamber, bool sort)
00214  {
00215      return this->getChannelHitOrder(aChannelHit.getChannel(),aChamber, sort);
00216  }
00217  //------------------------------------------------------------------------------------
00218  ui64 Event::getChannelHitOrder(const Channel& achannel  ,const Chamber& aChamber)
00219  {
00220      multiset< ChannelHit*,sorter>* list = this->channelHits.find(aChamber.getId())->second;//channelHits[aChamber.getId()];
00221      ui64 order=0;
00222   // Find corresponding channelHit
00223      for (multiset<ChannelHit*,sorter>::const_iterator iter=list->begin(); iter!=list->end(); iter++)
00224      {
00225              if ( (*iter)->getChannel().getSoftId() == achannel.getSoftId())
00226              {
00227                  return order;
00228              }
00229              order++;  // incremente le numero d'ordre pour les channela de la chambre
00230 
00231     }
00232   throw MicroException("No Channel respond for this event, this order and this chamber ");
00233 
00234  }
00235         */
00236 /*------------------------------------------------------------------------------------*/
00237 void Event::setSlowControlParam(const SlowControlEntry* params)
00238 {
00239   slowControlEntry = params;
00240   FILE_LOG(logDEBUG1) << "Event temperature["<< slowControlEntry->getTemperature() << "]" << endl;
00241   FILE_LOG(logDEBUG1) << "Event overPressure["<< slowControlEntry->getOverPressure() << "]" << endl;
00242   FILE_LOG(logDEBUG1) << "Event Pressure["<< slowControlEntry->getPressure() << "]" << endl;
00243 }
00244 
00245 /*------------------------------------------------------------------------------------*/
00246 
00247 void Event::print(void) const
00248 {
00249   FILE_LOG(logINFO) << "-Event id["<< id << "]" << endl;
00250   FILE_LOG(logINFO) << "-Event timestamp[" << timestamp << "]"  << endl ;
00251 }
00252 
00253 /*------------------------------------------------------------------------------------*/
00254 const float Event::getPressure() const
00255 {
00256     if ( slowControlEntry != NULL)
00257     {
00258         return slowControlEntry->getPressure();
00259     }
00260     else
00261     {
00262         return 0;
00263     }
00264 
00265 }
00266 /*------------------------------------------------------------------------------------*/
00267 const float Event::getTemperature() const
00268 {
00269     if ( slowControlEntry != NULL)
00270     {
00271         return slowControlEntry->getTemperature();
00272     }
00273     else
00274     {
00275         return 0;
00276     }
00277 
00278 }
00279 /*------------------------------------------------------------------------------------*/
00280 const BoardSlowControl& Event::getBoardSlowControl(const ui32 chamberId,const ui32 difId, const ui32 boardId) const
00281 {
00282     if ( slowControlEntry != NULL)
00283     {
00284        try
00285        {
00286            return slowControlEntry->getBoard(chamberId,difId,boardId);
00287        }
00288        catch ( MicroException e) {}
00289     }
00290 
00291     throw MicroException("Event::getBoardSlowControl: board not found");
00292 }
00293 /*------------------------------------------------------------------------------------*/
00294 const float Event::getBoardDriftVolt(const ui32 chamberId,const ui32 difId, const ui32 boardId) const
00295 {
00296     if ( slowControlEntry != NULL)
00297     {
00298        try
00299        {
00300            return slowControlEntry->getBoard(chamberId,difId,boardId).getDrift().voltage;
00301        }
00302        catch ( MicroException e) {}
00303     }
00304     throw MicroException("Event::getBoardSlowControl: board not found");
00305 
00306 }
00307 /*------------------------------------------------------------------------------------*/
00308 const float Event::getBoardMeshVolt(const ui32 chamberId,const ui32 difId, const ui32 boardId) const
00309 {
00310     if ( slowControlEntry != NULL)
00311     {
00312        try
00313        {
00314            return slowControlEntry->getBoard(chamberId,difId,boardId).getMesh().voltage;
00315        }
00316        catch ( MicroException e) {}
00317     }
00318     //throw ( MicroException("No mesh voltage for this chamber"));
00319     return -1;
00320 
00321 }
00322 //------------
00323 ui64 Event::getChannelHitSize() const
00324 { 
00325     ui64 size = 0;
00326     for ( std::map<ui16, std::multiset<ChannelHit*,sorter>* >::const_iterator it=channelHits.begin(); it != channelHits.end() ; it++)
00327     {
00328             size +=it->second->size();
00329     }
00330     return size;
00331 }
00332 
00333 //------------
00334 void Event::setTemperature(ui16 difId, float asu1, float asu2, float dif)
00335 {
00336   // Search or create difs container
00337     if ( temperatureAsu1.find(difId) != temperatureAsu1.end())
00338     { // dif already exit/ replace value
00339        temperatureAsu1[difId] = asu1;
00340     }
00341     else
00342     { // does not exit create it
00343       temperatureAsu1.insert(make_pair(difId,asu1));
00344     }
00345   // Search or create difs container
00346     if ( temperatureAsu2.find(difId) != temperatureAsu2.end())
00347     { // dif already exit/ replace value
00348        temperatureAsu2[difId] = asu2;
00349     }
00350     else
00351     { // does not exit create it
00352       temperatureAsu2.insert(make_pair(difId,asu2));
00353     }
00354   // Search or create difs container
00355     if ( temperatureDif.find(difId) != temperatureDif.end())
00356     { // dif already exit/ replace value
00357        temperatureDif[difId] = dif;
00358     }
00359     else
00360     { // does not exit create it
00361       temperatureDif.insert(make_pair(difId,dif));
00362     }
00363 
00364 }
00365 
00366 
00367 
00368 /*------------------------------------------------------------------------------------*/
00369 /*
00370 LCEventImpl* Event::getLcioImpl() 
00371 {
00372 try {
00373         LCEventImpl* evt = new LCEventImpl();
00374 
00375         evt->setRunNumber(  1   ) ;
00376     evt->setDetectorName( "m2" ) ;
00377     evt->setWeight( 1 ) ;
00378 
00379 
00380     evt->setEventNumber( this->id ) ;
00381     // Event timestamp express in millisecond
00382         // But lcio timestamp express in second
00383     evt->setTimeStamp( this->timestamp * 1000000  ) ;
00384 
00385 
00386     evt->parameters().setValue("Dif_syncro",this->difSynchro ) ;
00387     evt->parameters().setValue("Temperature",this->getTemperature() ) ;
00388     evt->parameters().setValue("Pressure",this->getPressure()) ;
00389     evt->parameters().setValue("Global trigger counter",(int)this->globalTriggerCounter ) ;
00390     evt->parameters().setValue("Dif trigger counter",(int)this->difTriggerCounter ) ;
00391     //evt->parameters().setValue("SetDriftVoltage",this->difTriggerCounter ) ;
00392     //evt->parameters().setValue("SetMeshVoltage",this->difTriggerCounter ) ;
00393 
00394 
00395 
00396     // now add some raw calorimeter hits
00397     LCCollectionVec* calVec = new LCCollectionVec( LCIO::RAWCALORIMETERHIT )  ;
00398 
00399     LCFlagImpl chFlag(0) ;
00400     //chFlag.setBit( LCIO::CHBIT_ID1 ) ;
00401 
00402     chFlag.setBit( LCIO::RCHBIT_LONG )   ;  // long(1) - short(0) , incl./excl. position
00403     chFlag.setBit( LCIO::RCHBIT_BARREL ) ;  // barrel(1) - endcap(0)
00404     //chFlag.setBit( LCIO::RCHBIT_TIME   ) ; // 1: time information stored
00405     //chFlag.setBit( LCIO::RCHBIT_ID1   )  ;  // cellid1 stored
00406     //chFlag.setBit( LCIO::RCHBIT_NO_PTR ) ; // 1: pointer tag not added
00407     //chFlag.setBit( LCIO::RCHBIT_ENERGY_ERROR )  ;   // 1: store energy error
00408     calVec->setFlag( chFlag.getFlag()) ;
00409 
00410 
00411 //
00412 //OFFSET_CHAMBER *NB_CHANNEL* GetChamberId()
00413 //+ OFFSET_DIF *NB_CHANNEL*  GetDifId()
00414 //+ OFFSET_CHIP*NB_CHANNEL * GetChipId ()
00415 //+  channelId 
00416 //
00417         // Chambre 1024 -> 2 exp 10
00418         // Dif     128  -> 2 exp 7
00419         // Chip    1024 -> 2 exp 10
00420         // Channel 1024 -> 2 exp 10
00421     std::string cellIDEncoding( "CHAMBER:10,DIF:7,CHIP:10,CHANNEL:10") ;
00422     CellIDEncoder<RawCalorimeterHitImpl> cellEncoder( cellIDEncoding , calVec ) ;
00423 
00424 
00425     ChamberMap_t chamberList = this->getRun().getDetector().getChambers();
00426     for (ChamberMap_t::iterator iter=chamberList.begin(); iter!=chamberList.end(); iter++)
00427     {
00428       Chamber& chamber = *(iter->second);
00429       try
00430       {
00431           const  multiset<ChannelHit*,sorter>& list = this->getChannelHitVector(chamber);
00432           int order = 0 ;
00433           for (set<ChannelHit*,sorter>::iterator iter=list.begin(); iter!=list.end(); iter++)
00434           {
00435             ChannelHit& channelHit = *(*iter);
00436             const Channel& channel = channelHit.getChannel();
00437                 RawCalorimeterHitImpl* rawHit = NULL ;
00438             try
00439             {
00440                           // Create RawCalorimeterHit
00441                   rawHit = new RawCalorimeterHitImpl() ;
00442 
00443               if ( chamber.getType().compare(GASSIPLEXCHAMBER1)==0 ||  chamber.getType().compare(GASSIPLEXCHAMBER4)==0)
00444               {
00445                     channelHit.setTimestamp(this->getTimeStamp());
00446                           }
00447                         
00448                         
00449                   cellEncoder["CHAMBER"] = channel.getChamber().getId(); 
00450                   cellEncoder["DIF"] = channel.getChip().getBoard().getDif()->getId(); 
00451                   cellEncoder["CHIP"] = channel.getChip().getId(); 
00452                   cellEncoder["CHANNEL"] = channel.getHardId(); 
00453                           cellEncoder.setCellID( rawHit ) ;
00454                         
00455                   rawHit->setAmplitude( channelHit.getDigitalValue()   )  ;
00456                   //rawHit->setTimeStamp( channelHit.getTimestamp() / 1000) ; 
00457 
00458               order++;  // because aEvent have already been  sorted
00459                   calVec->push_back( rawHit ) ;
00460             } // end try block
00461             catch (MicroException e)
00462             {
00463             }
00464 
00465           } // end for channelHits for this chamber
00466       }
00467       catch (MicroException e)
00468       {
00469                 // get ChamberHitVector for next chamber
00470       }
00471 
00472     }
00473     //LCTOOLS::printRawCalorimeterHits(calVec);
00474     // add all collections to the event
00475 
00476     evt->addCollection( calVec , "m2" ) ;
00477         return evt;
00478 
00479 
00480   }
00481   catch( Exception& ex){
00482 
00483     cout << " an excpetion occured: " << endl ;
00484     cout << "  " << ex.what() << endl ;
00485     return NULL;
00486   }
00487 
00488     return NULL;
00489 
00490 
00491 }*/
00492 bool Event::isValid(void)
00493 {
00494   if (( valid & VALID ) == VALID) return true;
00495   return false;
00496 }
00497 
00498 //_____________________________________________________________________________
00499 void Event::setValidFlag(ui16 value,bool isValid) 
00500 {
00501     valid = valid | value; 
00502     if ( ! isValid ) { valid = valid & 0xFFFE ; }
00503 }
00504 //_____________________________________________________________________________
00505 string Event::getValidFlagString(void ) const
00506 {
00507 string result;
00508   if ( (valid & VALID )               == VALID )               { result = result + " " + "VALID" ;}
00509   if ( (valid & UNKNOW_NOT_VALID )    == UNKNOW_NOT_VALID )    { result = result + " " + "UNKNOW_NOT_VALID" ;}
00510   if ( (valid & BAD_CHIP_ID )         == BAD_CHIP_ID )         { result = result + " " + "BAD_CHIP_ID" ;}
00511   if ( (valid & BAD_BOARD_ID )        == BAD_BOARD_ID )        { result = result + " " + "BAD_BOARD_ID" ;}
00512   if ( (valid & BAD_DIF_ID )          == BAD_DIF_ID )          { result = result + " " + "BAD_DIF_ID" ;}
00513   if ( (valid & CRC_ERROR )           == CRC_ERROR )           { result = result + " " + "CRC_ERROR" ;}
00514   if ( (valid & ONLY_FFFF)            == ONLY_FFFF )           { result = result + " " + "ONLY_FFFF" ;}
00515   if ( (valid & ONLY_AAAA)            == ONLY_AAAA )           { result = result + " " + "ONLY_AAAA" ;}
00516   if ( (valid & ONLY_5555)            == ONLY_5555 )           { result = result + " " + "ONLY_5555" ;}
00517   if ( (valid & GLOBAL_HEADER_ERROR)  == GLOBAL_HEADER_ERROR ) { result = result + " " + "GLOBAL_HEADER_ERROR" ;}
00518   if ( (valid & FRAME_HEADER_ERROR)   == FRAME_HEADER_ERROR )  { result = result + " " + "FRAME_HEADER_ERROR" ;}
00519   if ( (valid & TRAILER_C3_FOUND)     == TRAILER_C3_FOUND )    { result = result + " " + "TRAILER_C3_FOUND" ;}
00520   if ( (valid & END_OF_FILE_ERROR)    == END_OF_FILE_ERROR )   { result = result + " " + "END_OF_FILE_ERROR" ;}
00521 
00522 
00523   return result;
00524 }
00525
/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/event/Event.cpp
