/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/lcio/createLcioCalorimeterHit.cc

Go to the documentation of this file.

/* @version $Revision: 1560 $ * @modifiedby $Author: jacquem $ * @lastmodified $Date: 2012-03-06 09:54:17 +0100 (Tue, 06 Mar 2012) $ */

#include "IO/LCWriter.h"
#include "IMPL/LCRunHeaderImpl.h"
#include "IMPL/LCEventImpl.h"
#include "IMPL/LCTOOLS.h"
#include "IMPL/LCEventImpl.h"
#include "IMPL/CalorimeterHitImpl.h"
#include "IMPL/LCCollectionVec.h"
#include "IMPL/LCFlagImpl.h"

#include "lcio.h"

#include "tools/MicroException.hh"
#include "tools/Arguments.hh"

#include <sstream>
#include <iostream>

#include <TRandom.h>
#include <TTree.h>
#include <TFile.h>
#include <TStyle.h>
#include <TKey.h>
#include <TreeClass.hh>

#include <Riostream.h>
#include <sstream>
#include <string>
#include "mysql/Mysql.hh"
#include "MicroException.hh"
#include <iostream>
#include <fstream>

//#include <TROOT.h>
//#include <TRint.h>
#include "TH2I.h"
//#include "TNtuple.h"

using namespace lcio;

void selectAndCopyEvent(const string& inputFileName,const string& outputFileName);
Long64_t buildCaloEvent(const MTEvent& evt, LCEventImpl&, Long64_t bcidAbsPrev, bool &track,TH1F& hResidualX, TH1F& hResidualY);

void usage(void)
{
    cout  << "     : buildCaloHits -f filename outputfileName" << endl;
}

int main(int argc, char** argv) {

  if ( argc < 4 || argc > 5 )
  {
    usage();
    return -1;
  }

  // Get input file name in a input file
  if ( strcmp (argv[1] ,"-f") == 0 )
  {
     if ( argc != 4 ) { usage(); return -1; } 
     string outputFileName = argv[3];

     ifstream myReadFile;
     myReadFile.open(argv[2]);
     string inputRootFile;
     if (myReadFile.is_open()) {
     while (!myReadFile.eof()) {
        getline(myReadFile,inputRootFile);
        if ( inputRootFile.size() > 0 ) 
        {
          selectAndCopyEvent(inputRootFile,outputFileName);
        }
     }
    }
    myReadFile.close();
    return 0;
  }
  else 
  {
    usage();
    return -1;
  }
  return 0;
}


  /*-----------------------------------------------------------------------------------------------*/

void selectAndCopyEvent(const string& inputFileName,const string& outputFileName)
{

 
  TFile readFile(inputFileName.c_str(),"READONLY");

  ui16 runId = 0;

  // lcio File Writer 
  LCWriter* lcWrt = LCFactory::getInstance()->createLCWriter()  ;
  lcWrt = LCFactory::getInstance()->createLCWriter()  ;

  // Lcio Run header
  LCRunHeaderImpl* runHdr = new LCRunHeaderImpl() ;
  runHdr->setDetectorName( "micromegas microroc m2 and Gassiplex telescope" ) ;
  runHdr->setDescription( " testing lcio with micromegas" ) ;
  /*
string ecalName("micromegas_microroc_m2_calorimeter") ;
  runHdr->addActiveSubdetector( ecalName ) ;

  string telName("telecsope") ;
  runHdr->addActiveSubdetector( telName ) ;
*/
  TIter nextkey(readFile.GetListOfKeys());
  TKey *key = NULL;
  
  while (key= (TKey*)nextkey() )
  {
    TTree *readTree = (TTree*)key->ReadObj();

    MTEvent *readEvt =  new MTEvent();

    TBranch *branch= readTree->GetBranch("MTEvent");

    branch->SetAddress(&readEvt);

  // Open lcio new file in write mode
    lcWrt->open( outputFileName , LCIO::WRITE_NEW )  ;

    Double32_t temp=0;
    Double32_t pres=0;

    Long64_t bcIdAbsPrev = -1;
    runHdr->setRunNumber(runId++ ) ;

    // Loop over event in current input TTree
    int nbEntries = readTree->GetEntries();
    ui16 evtNum = 0;
    TH1F hResidualX("hResidualX" , "hResidualX;xtel - xm2 (pad)" ,201,-100,100); 
    TH1F hResidualY("hResidualY" , "hResidualY;ytel - ym2 (pad)" ,201,-100,100); 

    for ( evtNum = 0; evtNum < nbEntries ; evtNum++)
    {
      // lcio event
      LCEventImpl* lcioEvtImpl =  new LCEventImpl();
      // Get input Event throw readEvt
      readTree->GetEntry(evtNum);
      if ( evtNum % 1 == 0)
      {
        cout << "Done for entry " << evtNum  +1 << " / " << nbEntries << "\r" << flush;
      }

      temp = temp + readEvt->GetTemperature();
      pres = pres + readEvt->GetPressure();
      bool track = false;
      bcIdAbsPrev = buildCaloEvent(*readEvt,*lcioEvtImpl,bcIdAbsPrev, track, hResidualX,hResidualY);
      if( track )
      {
        lcWrt->writeEvent( lcioEvtImpl ) ;
      } 
      delete lcioEvtImpl;
    }
    pres = pres / evtNum;
    temp = temp / evtNum;
    //run.SetTemperature(temp);
    //run.SetPressure(pres);
    
    //run.SetOffsetX(-101 + hResidualX.GetMaximumBin()); 
    //run.SetOffsetY(-101 + hResidualY.GetMaximumBin()); 

  
    lcWrt->close() ;
    delete lcWrt;
    delete readEvt;
  }
  readFile.Close();

 cout << endl;
}





//
// Build Tree with position, time and digit of M2 prototype hits
// + cut on time information (last readout less than 3.35 s before)
// + cut on high multiplicity in m2 prototype inside dt window from 0-20 bcid
// + cut on the dt distribution between 5 and 10



  Long64_t buildCaloEvent(const MTEvent& evt, LCEventImpl& caloEvent, Long64_t bcIdAbsPrev ,bool& track,TH1F& hResidualX, TH1F& hResidualY)
  {

  caloEvent.setTimeStamp( evt.GetTimestamp() * 1000000  ) ;
  caloEvent.setDetectorName( "Lapp Microroc Micromegas and Gassiplex telecope" ) ;


  int nchannel = 0;
  //Define bcIdAbsPrev, bcid_abs of first readout
  if (bcIdAbsPrev == -1)
  {
    nchannel = evt.GetNchannel();

    for (int j=0;j<nchannel;j++)
    {             
      MTChannel* channel = (MTChannel*)evt.GetChannels()->UncheckedAt(j);

      if (channel->GetChamberId() == 10)
      {
       return  channel->GetBcId_Abs();
      }
    }
  }


  /*m2 prototype time cut*/
  int dt_min = 4;
  int dt_max = 10;

  /*spark cut = Nhit inside 0-20 dt region*/
  int nhit_spark_cut = 1000;

  /*telescope energy cut*/
  int adc_cut = 25;

  int i = 0;
  int t = 0;
  int adc = 0;
  int xpos = 0;
  int ypos = 0;
  int zpos = 0;
  int zmax = 0;
  int digit = 0;
  int chbid = 0;
  int nhit_m2_cut = 0;
  double time = 0;
  TString name,title;
  Long64_t tprev,t0,t1,t2,t3,dt;

  int nhit[3];
  int xmax[3];
  int ymax[3];

  TH2I * hxy_tel_adc[3];


  for (int i=0;i<3;i++)
  {
    nhit[i] = 0;
    xmax[i] = 0;
    ymax[i] = 0;

    name = Form("hxy_tel_adc_chamber_%i",i+7);
    title = Form("Hit pattern (adc>%i) chamber %i;y (cm);x (cm)",adc_cut,i+7);
    hxy_tel_adc[i] = new TH2I(name,title,32,0,32,32,0,32);
  }

  //if t0 is defined start loop over events

  nchannel = evt.GetNchannel();

  if (nchannel!=0)
  {
    //determine time to previous readout
    for (int j=0;j<nchannel;j++)
    {
      MTChannel* channel = (MTChannel*)evt.GetChannels()->UncheckedAt(j);
      if (channel->GetChamberId() == 10)
      {
        t1 = channel->GetBcId_Abs();
        j = nchannel;
      }
    }

    time = t1 - bcIdAbsPrev;
    time *= 200e-9;
    bcIdAbsPrev = t1;


    //if readout occured after less than 3.35 s of previous readout
    //skip also first event as hit time info is wrong
    if (time <= 3.35)
    {
      //look for a track in telescope
      for (int j=0;j<nchannel;j++)
      {
        MTChannel* channel = (MTChannel*)evt.GetChannels()->UncheckedAt(j);
        chbid = channel->GetChamberId();

        //look in pad chambers
        if (chbid == 7 || chbid == 8 || chbid == 9)
        {
          adc = channel->GetAnalogValue();

          if (adc >= adc_cut)
          {
            nhit[chbid-7]++;
            hxy_tel_adc[chbid-7]->Fill(channel->GetY()*1e-4,channel->GetX()*1e-4,adc);
          }
        }
      }

      //require that 3 pad chambers have at least 1 hit
      if (nhit[0]>0 && nhit[1]>0 && nhit[2]>0)
      {
        //find position of ADC maxima

        for (int k=0;k<3;k++)
        {
          hxy_tel_adc[k]->GetMaximumBin(ymax[k],xmax[k],zmax);
        }

        //require aligned maxima

        if (ymax[0]==ymax[1] && ymax[0]==ymax[2] &&
            xmax[0]==xmax[1] && xmax[0]==xmax[2])
        {
          //reloop over hits
          nhit_m2_cut = 0;

          for (int j=0;j<nchannel;j++)
          {
            MTChannel* channel = (MTChannel*)evt.GetChannels()->UncheckedAt(j);

            //look in m2 prototype
            if (channel->GetChamberId() == 10)
            {                                 
              t2 = channel->GetBcId_Dif();
              t3 = channel->GetBcId_Hit();
              dt = t2 - t3;

              if (dt>=0 && dt<=20)
              {
                nhit_m2_cut++;
              }
            }
          }

          //reject high multiplicity event in time region of signal
          if (nhit_m2_cut < nhit_spark_cut)
          {
                  track = 1;    
            caloEvent.parameters().setValue("XTelescope", xmax[0] );
            caloEvent.parameters().setValue("YTelescope", ymax[0] );

            //reloop over hits
            LCCollectionVec* calVec = new LCCollectionVec( LCIO::CALORIMETERHIT )  ;
            for (int j=0;j<nchannel;j++)
            {
              MTChannel* channel = (MTChannel*)evt.GetChannels()->UncheckedAt(j);

              //look in m2 prototype
              if (channel->GetChamberId() == 10)
              {
                t3 = channel->GetBcId_Hit();
                dt = t2 - t3;

                if (dt>=dt_min && dt<=dt_max)
                {
                  digit = channel->GetDigitalValue();
                  xpos = (int)channel->GetX();
                  ypos = (int)channel->GetY();
                  zpos = (int)channel->GetZ();

                  xpos /= 10000;
                  ypos /= 10000;

                  t = (int)dt;

//                  CaloHit hit;
                        hResidualX.Fill(xmax[0]-xpos);
                        hResidualY.Fill(ymax[0]-ypos);
/*
                  hit.SetDeltaT(t);
                  hit.SetThreshold(digit);
*/

                  CalorimeterHitImpl* caloHit = new CalorimeterHitImpl() ;
                  //caloHit->parameters().setValue("DeltaT", t);
                  //caloHit->parameters().setValue("Threshold", digit);
                  try
                  {
                    // Create RawCalorimeterHit
                    caloHit->setEnergy( digit )  ;
                    caloHit->setTime( t )  ;
                    const float pos[3] = { xpos, ypos, 0};
                    caloHit->setPosition(pos);



                      LCFlagImpl chFlag(0) ;
                      //chFlag.setBit( LCIO::CHBIT_ID1 ) ;
                      chFlag.setBit( LCIO::RCHBIT_LONG )   ;  // long(1) - short(0) , incl./excl. position
                      chFlag.setBit( LCIO::RCHBIT_BARREL ) ;  // barrel(1) - endcap(0)
                      chFlag.setBit( LCIO::RCHBIT_TIME   ) ; // 1: time information stored
                      //chFlag.setBit( LCIO::RCHBIT_ID1   )  ;  // cellid1 stored
                      chFlag.setBit( LCIO::RCHBIT_NO_PTR ) ; // 1: pointer tag not added
                      //chFlag.setBit( LCIO::RCHBIT_ENERGY_ERROR )  ;   // 1: store energy error
                      calVec->setFlag( chFlag.getFlag()) ;

                    calVec->push_back( caloHit ) ;

                  } // end try block
                  catch (MicroException e)
                  {
                  }
                  //LCTOOLS::printRawCalorimeterHits(calVec);
                  // add all collections to the event

                }
              }
            }
            caloEvent.addCollection( calVec , "m2" ) ;
          }
          else { track = 0; }
        }
      } //end nhit[0]>0 && nhit[1]>0 && nhit[2]>0
    } // end time <= 3.35 && i != 0
  } // enf if nChannel != 0
    
  for (int i=0;i<3;i++)
  {
   delete hxy_tel_adc[i];
  }

  return bcIdAbsPrev;
}

---