/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/root/buildCaliceCaloHits.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 "tools/Arguments.hh"

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

void selectAndCopyEvent(const string& inputFileName,const string& outputFileName,UInt_t acceptation,UInt_t limit);
bool buildCaliceEvent(MTEvent& evt, CaloEvent& caloEvent,UInt_t dtCanditade,UInt_t acceptation);

// algoritm for calice cubic meter in ram full mode

void usage(void)
{
    cout  << "     : buildCaloHits inputfilename outputfileName acceptation -l eventNum" << endl;
}

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

  
  Arguments arg(argc,argv,"buildCaloHits inputfilename outputfileName acceptation -l eventNum");
  if (arg.getNbOfArgs() != 3 || arg.getNbOfOptions() > 1 )
  {
    arg.usage();
    return -1;
  }

  // Get input file name in a input file
  string outputFileName = arg.getArgument(2);

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

  string inputRootFile = arg.getArgument(1);
  cout<<"Read ROOT file : "<<inputRootFile<<endl;
  UInt_t acceptation = atoi(arg.getArgument(3).c_str());

  unsigned int limit = 0;
  if ( arg.getOption("-l").size() != 0 )
  {
    limit = atoi(arg.getOption("-l").c_str());
  } 

  cout << inputRootFile << endl;
  cout << outputFileName << endl;

  selectAndCopyEvent(inputRootFile,outputFileName,acceptation,limit);
  return 0;
}


/*-----------------------------------------------------------------------------------------------*/
void selectAndCopyEvent(const string& inputFileName,const string& outputFileName,UInt_t acceptation,UInt_t limit)
{
 
  TFile readFile(inputFileName.c_str(),"READONLY");
  TFile caloFile(outputFileName.c_str(), "UPDATE","Micromegas calorimeter data");

  ui32 nbEmptyEvent = 0;


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

  UInt_t eventId = 0;
  while (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();
    cout << "read MTRun id[ " << readRun->GetRunId() << "]" << endl;
    float offset_x = 0;
    float offset_y = 0;

    // 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(40000000000); // 40 GB



    CaloEvent *evt = new CaloEvent();
    // create a branch with run 
    tree.Branch("CaloEvent",&evt,16000,2);
    run.SetRunId(runId);
    cout << "write CaloRun id[ " << run.GetRunId() << "]" << endl;

    // Loop over event in current input TTree
    int nbEntries = readTree->GetEntries();
    ui16 evtNum = 0;

    if ( limit == 0 ) { limit = nbEntries; }

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

      // Get input Event as readEvt
      readTree->GetEntry(evtNum);
      if ( readEvt->IsValid() == true )
      {
        const std::vector<UInt_t> &dtCanditade = readEvt->GetDtCandidate();
        for ( vector<UInt_t>::const_iterator iter =  dtCanditade.begin() ; iter != dtCanditade.end() ; iter++ )
        { 
          UInt_t candidate = *iter;

          evt->SetDeltaTmax(candidate);

          buildCaliceEvent(*readEvt,*evt,candidate,acceptation);
        //  cout << "evt->GetNHits()[" << evt->GetNHits() <<  "]"  << endl;
          if ( evt->GetNHits() == 0 )
          {
            nbEmptyEvent++;
          }
          else 
          {
            evt->SetEventId(eventId++);
            tree.Fill();
          }
          evt->Clear();
        }
      }
    }


    tree.GetUserInfo()->Add(&run); 
    TFile *curFile =tree.GetCurrentFile() ;
    //curFile->Write("", TObject::kOverwrite);
    tree.Write("", TObject::kOverwrite);
    cout << "Write TTree" <<endl;
    //tree.Write();
    delete readEvt;
    delete evt;

  }
  caloFile.Close();
  readFile.Close();
  cout << "Done: " << eventId << " events reconstructed. Empty event:"  << nbEmptyEvent <<  endl;
}



// Construct new CaloEvent from MTEvent for dtCandidate
bool buildCaliceEvent(MTEvent& evt, CaloEvent& caloEvent, UInt_t dtCanditade, UInt_t acceptation)
{
  //loop over hits
  int nchannel = evt.GetNchannel();

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

    Long64_t dt = channel->GetBcId_Dif() - channel->GetBcId_Hit();

    if (fabs(dt-dtCanditade) <= acceptation)
    {
      int digit = channel->GetDigitalValue();
      int xpos = (int)channel->GetX();
      int ypos = (int)channel->GetY();
      int zpos = (int)channel->GetZ();
      int layer = channel->GetChamberId();

      const MTChannelSoftId &soft = channel->GetSoftId();

      CaloHit hit;

      hit.SetX(xpos);
      hit.SetY(ypos);
      hit.SetZ(zpos);
      hit.SetRowInChamber(channel->GetRowInChamber());
      hit.SetColInChamber(channel->GetColInChamber());
      hit.SetDeltaT(dt);
      hit.SetThreshold(digit);
      hit.SetLayer(layer);
      hit.SetSoftId(soft);


      caloEvent.AddHit(hit);
    }
  }
  return true;
}



// 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 != NULL )
    {
      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);
  }//  end if track exist


  

    //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->GetBcId_Dif() - channel->GetBcId_Hit();
              
            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);  
   }
  }
 }
}

---