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

#include "IO/LCWriter.h"
#include "IMPL/LCRunHeaderImpl.h"
#include "IMPL/LCEventImpl.h"
#include "IMPL/LCTOOLS.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 <string>
#include "mysql/Mysql.hh"
#include <fstream>
#include "TH2I.h"

Include dependency graph for createLcioCalorimeterHit.cc:

Go to the source code of this file.

Functions
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)
int	main (int argc, char **argv)

---

Function Documentation

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

Long64_t buildCaloEvent	(	const MTEvent &	evt,
		LCEventImpl &	,
		Long64_t	bcidAbsPrev,
		bool &	track,
		TH1F &	hResidualX,
		TH1F &	hResidualY
	)

Definition at line 192 of file createLcioCalorimeterHit.cc.

  {

  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;
}

void usage

(

void

)

Definition at line 45 of file createLcioCalorimeterHit.cc.

{
    cout  << "     : buildCaloHits -f filename outputfileName" << endl;
}

int main	(	int	argc,
		char **	argv
	)

Definition at line 50 of file createLcioCalorimeterHit.cc.

                                {

  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;
}

---