/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/lcio/buildCaloHits.cpp

Go to the documentation of this file.

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

#include <Riostream.h>
#include <sstream>
#include <string>
#include "mysql/Mysql.hh"
#include "MicroException.hh"
#include <iostream>
#include <fstream>
#include "root/CaloHit.hh"
#include "root/MTRun.hh"
#include "root/MTEvent.hh"
#include "root/MTChannel.hh"
#include "root/MTMicrorocChip.hh"


#include "TH2I.h"


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

void usage(void)
{
    cout  << "usage: buildCaloHits -t timestampMin timestampMax outputfileName" << endl;
    cout  << "     : buildCaloHits -f filename outputfileName" << endl;
    cout << "1312577280 1312590660" << endl;
}

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

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


  // Get inputFile name in database by timestamp limit
  if ( strcmp (argv[1] ,"-t" )== 0 )
  {
    if ( argc != 5 ) { usage(); return -1; } 
    string outputFileName = argv[4];
    Mysql sql;
    string start = argv[2];
    string end = argv[3];
    
    // RECREATE Root file to clear its contains if it alread exist
//    TFile caloFile(outputFileName.c_str(), "RECREATE","Micromegas calorimeter data");

    try 
    {
      // Get list of TFile to filter from Database
      const vector<string>& ids = sql.getMergedFilesByTimestamp(start,end);
      vector<string>::const_iterator cii;

      // Loop over input TFile
      for(cii=ids.begin(); cii!=ids.end(); cii++)
      {

        string inputFileName = *cii;
        selectAndCopyEvent(inputFileName,outputFileName);
      }
      delete &ids;
    }
    catch (MicroException e) { }
  }

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

     // RECREATE Root file to clear its contains if it alread exist
     TFile caloFile(outputFileName.c_str(), "RECREATE","Micromegas calorimeter data");
     caloFile.Write();
     caloFile.Close();

     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;

  TFile caloFile(outputFileName.c_str(), "UPDATE","Micromegas calorimeter data");
  TIter nextkey(readFile.GetListOfKeys());
  TKey *key = NULL;
  CaloRun run;
  while (key= (TKey*)nextkey() )
  {
    TTree *readTree = (TTree*)key->ReadObj();

    MTEvent *readEvt =  new MTEvent();
    TBranch *branch= readTree->GetBranch("MTEvent");
    branch->SetAddress(&readEvt);

    // Create a new TTree which will contain  filtered events
    stringstream label;
    label<< inputFileName << "_" << runId ; 
    TTree tree(label.str().c_str(),"Purged MicroMegas event");
    cout << "work on :" << label.str() << endl;
    tree.SetMaxTreeSize(4000000000); // 4 Mb

    Double32_t temp=0;
    Double32_t pres=0;

    Long64_t bcIdAbsPrev = -1;

    CaloEvent *evt = new CaloEvent();
    // create a branch with run 
    tree.Branch("CaloEvent",&evt,16000,2);

    run.SetRunId(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++)
    {
      // 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,*evt,bcIdAbsPrev, track, hResidualX,hResidualY);
      if( track )
      {
        tree.Fill();  // fill the tree with the current event
        evt->Clear();
      } 
    }
    pres = pres / evtNum;
    temp = temp / evtNum;
    run.SetTemperature(temp);
    run.SetPressure(pres);
    
    run.SetOffsetX(-101 + hResidualX.GetMaximumBin()); 
    run.SetOffsetY(-101 + hResidualY.GetMaximumBin()); 

    tree.GetUserInfo()->Add(&run); 
    TFile *curFile =tree.GetCurrentFile() ;
    //curFile->Write("", TObject::kOverwrite);
    tree.Write();
  
    delete readEvt;
    delete evt;
  }
  caloFile.Close();
  readFile.Close();
 // delete run;
}





//
// 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, CaloEvent& caloEvent, Long64_t bcIdAbsPrev ,bool& track,TH1F& hResidualX, TH1F& hResidualY)
  {

  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.SetXTelescope(xmax[0]);
            caloEvent.SetYTelescope(ymax[0]);
            //reloop over hits
            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;
                  hit.SetX(xpos);
                  hit.SetY(ypos);
                        hResidualX.Fill(xmax[0]-xpos);
                        hResidualY.Fill(ymax[0]-ypos);
                  hit.SetZ(0);
                  hit.SetDeltaT(t);
                  hit.SetThreshold(digit);
                  caloEvent.AddHit(hit);        
                }
              }
            }
          }
          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;
}

---