/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/root/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 <TROOT.h>
//#include <TRint.h>
#include "TH2I.h"
//#include "TNtuple.h"


#include "root/CaloRun.hh"
#include "root/CaloEvent.hh"
#include "root/MTRun.hh"
#include "root/MTEvent.hh"
#include "root/MTChannel.hh"
#include "root/MTTrack.hh"

void selectAndCopyEvent(const string& inputFileName,const string& outputFileName,int deltaTMin, int deltaTMax);
void buildCaloEvent(MTEvent& evt, CaloEvent&, int deltaTMin, int deltaTMax, float offset_x, float offset_y);

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

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

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


  int deltaTMin = -1;
  int deltaTMax = -1;





  // Get inputFile name in database by timestamp limit
  if ( strcmp (argv[1] ,"-t" )== 0 )
  {
    if ( argc < 5 ) { usage(); return -1; } 

    if ( argc == 7 ) 
    {
      deltaTMin = atoi ( argv[5] );
      deltaTMax = atoi ( argv[6] );
    }
    
    string outputFileName = argv[4];
    // RECREATE Root file to clear its contains if it alread exist
    TFile caloFile(outputFileName.c_str(), "RECREATE","Micromegas calorimeter data");
    caloFile.Write();
    caloFile.Close();

    Mysql sql;
    string start = argv[2];
    string end = argv[3];
    
    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,deltaTMin,deltaTMax);
      }
      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];

    if ( argc == 6 ) 
    {
      deltaTMin = atoi ( argv[4] );
      deltaTMax = atoi ( argv[5] );
    }

    



    // 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 ) 
        {
          cout<<"Read ROOT file : "<<inputRootFile<<endl;
          selectAndCopyEvent(inputRootFile,outputFileName,deltaTMin,deltaTMax);
        }
      }
    }
    myReadFile.close();
    return 0;
  }
  else 
  {
    usage();
    return -1;
  }
  return 0;
}










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

void selectAndCopyEvent(const string& inputFileName,const string& outputFileName,int deltaTMin, int deltaTMax)
{

 
  TFile readFile(inputFileName.c_str(),"READONLY");
  TFile caloFile(outputFileName.c_str(), "UPDATE","Micromegas calorimeter data");

  ui16 runId = 0;


  TIter nextkey(readFile.GetListOfKeys());
  TKey *key = NULL;
  CaloRun run;

  key= (TKey*)nextkey();
  

  //Get tree to read
  TTree *readTree = (TTree*)key->ReadObj();
  MTEvent *readEvt =  new MTEvent();
  TBranch *branch= readTree->GetBranch("MTEvent");
  branch->SetAddress(&readEvt);


  //Get RunId
  MTRun * readRun  = (MTRun*)readTree->GetUserInfo()->FindObject("MTRun");
  int runid = readRun->GetRunId();
  if (runid==0)
    {
      //cout<<"RunId not known, please enter: ";
      //cin>>runid;
      //cout<<endl;
      //cout<<endl;
      runid = 10227;
    }

  //Get telescope/prototype offset from database
  Mysql mysql;
  float offset_x = mysql.getRunOffsetX(runid);
  float offset_y = mysql.getRunOffsetY(runid);

  cout<<"OFFSETS = "<<offset_x<<" "<<offset_y<<endl;



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

  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;
  cout<<"Tree entries = "<<nbEntries<<endl;

  MTTrack * trk = NULL;


  //for ( evtNum = 0; evtNum < 100 ; evtNum++)
  for ( evtNum = 0; evtNum < nbEntries ; evtNum++)
    {
      cout << "Done for entry " << evtNum + 1 << " / " << nbEntries << "\r" << flush;

      // Get input Event as readEvt
      readTree->GetEntry(evtNum);

      bool valid_time = readEvt->GetTimeInfo();
cout << "readEvt->GetTimeInfo["<< readEvt->GetTimeInfo() << "]"  << endl;
      bool square_m2 = readEvt->GetSquareEvt();

      if (valid_time && !square_m2)
        {

          trk = readEvt->GetTrack();
cout << "readEvt->GetTrack()[" << readEvt->GetTrack() << "]" << endl;

          if (trk != NULL)
            {
              buildCaloEvent(*readEvt,*evt,deltaTMin,deltaTMax,offset_x,offset_y);
              tree.Fill();
              evt->Clear();
            }
        }
    }


  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 the dt distribution between 5 and 10

void buildCaloEvent(MTEvent& evt, CaloEvent& caloEvent,int deltaTMin, int deltaTMax, float offset_x, float offset_y)
  {

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

    if (deltaTMin != -1) { deltaTMin = dt_min; }
    if (deltaTMax != -1) { deltaTMax = dt_max; }


    //Distance from first telescope chamber to m2 prototype
    float dist_z = 156.;


    //Get track
    MTTrack * trk = evt.GetTrack();
    
    int trk_quality = 0;

    if (trk->GetFromPad())
      {
        trk_quality = 10;
      }

    if (trk->GetFromStrip())
      {
        trk_quality = 20;
      }

    if (trk->GetStraight())
      {
        trk_quality++;

        if (trk->GetSingleHit())
          {
            trk_quality++;
          }
      }

    
    caloEvent.SetTrkQuality(trk_quality);

    float extrapol_x = offset_x;
    float extrapol_y = offset_y;
    
    //Extrapolate at prototype
    
    //straight tracks
    if (trk_quality%10 != 0)
      {
        extrapol_x += trk->GetBx();
        extrapol_y += trk->GetBy();     
      }

    //inclined tracks
    else
      {
        extrapol_x += trk->GetAx() * dist_z + trk->GetBx();
        extrapol_y += trk->GetAy() * dist_z + trk->GetBy();
      }


    caloEvent.SetXExtrapol(extrapol_x);
    caloEvent.SetYExtrapol(extrapol_y);


  

    //loop over hits

    int nchannel = evt.GetNchannel();
    int 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)
          {
            Long64_t dt = channel->GetBcIdDif() - channel->GetBcIdHit();
              
            int zpos = (int)channel->GetZ() / 10000;

            if (dt>=dt_min && dt<=dt_max)
              {
                int digit = channel->GetDigitalValue();
                int xpos = (int)channel->GetX() / 10000;
                int ypos = (int)channel->GetY() / 10000;
                  
                int t = (int)dt;
                  
                CaloHit hit;

                hit.SetX(xpos);
                hit.SetY(ypos);           
                hit.SetZ(zpos);
                hit.SetDeltaT(t);
                hit.SetThreshold(digit);

                caloEvent.AddHit(hit);  
              }
          }
      }
  }

---