/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 AddCaloHit(CaloEvent& caloEvent,MTChannel* channel);


// 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");



  TIter nextkey(readFile.GetListOfKeys());
  TKey *key = NULL;
  CaloRun run;
  ui32 skipCandidate = 0;
  ui32 nbHits = 0;

  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++)
    {
      UInt_t lastCandidate = 0;
  //cout << "----lastCandidate[" << lastCandidate << "]"  << endl;

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


      if ( readEvt->IsValid() == true  )
      {
        bool findCaloEvent = false;
        int nbChannel = readEvt->GetNchannel();

        // Une fois les hits tries par absolutute time et les dtCanditade tries:
        // 1/  on prend le 1er  dtCandidate
        // 2/  on boucle sur les hits d'un MTEvent.
        //      tant que le hit est a +ou- acceptation du dtCanditate, on ajoute
        //        le CalHit correspondant au CaloEvent en cours 
        //     Des qu'on trouve un hit en dehors de l acceptation, on enregistre en CaloEvent
        //     on passe sur le dtCanditate suivant

        bool endOfCandidate = false;
        const std::vector<UInt_t> &dtCanditade = readEvt->GetDtCandidate();
        if ( dtCanditade.size() > 0 )
        {
          vector<UInt_t>::const_iterator iter =  dtCanditade.begin() ; 
          UInt_t candidate = *iter;
          lastCandidate = candidate;
  //cout << "----Candidate[" << candidate << "]"  << endl;
  //cout << "----lastCandidate[" << lastCandidate << "]"  << endl;


          for(int i=0;i<nbChannel && endOfCandidate == false ;i++)
          {
            MTChannel* channel = (MTChannel*)readEvt->GetChannels()->UncheckedAt(i);
            Long64_t dt = channel->GetBcIdDif() - channel->GetBcIdHit();
  //cout << "Channel[" << i << "/" << nbChannel <<  "], dt["<< dt << "] bcIdAbs[" << channel->GetBcIdAbs() << "]" <<  endl; 

            if (dt > 0 && fabs(dt-candidate) <= acceptation)
            { 
              AddCaloHit(*evt,channel);
              nbHits++;
              findCaloEvent = true;
  //cout << "Add CaloHit" << endl;
            }
            if (  i == nbChannel-1 )
            { // last channel
              evt->SetEventId(eventId++);
              evt->SetDeltaTmax(candidate);
              evt->SetBcIdAbs(channel->GetBcIdAbs());
  //cout <<  " lst event channel set CalEvent[" << eventId  << "]" << endl;
              tree.Fill();
              evt->Clear();
              findCaloEvent = false;
            }
            else if ( (fabs(dt-candidate) > acceptation && findCaloEvent == true) ) 
            {
              evt->SetEventId(eventId++);
              evt->SetDeltaTmax(candidate);
              evt->SetBcIdAbs(channel->GetBcIdAbs());
  //cout <<  "set CalEvent[" << eventId <<  "]  candidate[" << candidate  << "]" << endl;

              tree.Fill();
              evt->Clear();
              findCaloEvent = false;
              iter++;
              if ( iter == dtCanditade.end() ) { endOfCandidate = true; }
              while ( iter != dtCanditade.end()  &&  *iter <= (lastCandidate + acceptation  ) ) 
              {
  //cout << "----Skip Candidate[" << *iter << "] <=  (lastCandidate + acceptation )[" << (lastCandidate +  acceptation)  <<  "]"  << endl;

                 skipCandidate++;
                 iter++;
                 if ( iter == dtCanditade.end() ) { endOfCandidate = true; }
              }
              if ( ! endOfCandidate )
              {
                candidate = *iter;
                lastCandidate = candidate;
  //cout << "----Candidate[" << candidate << "]"  << endl;
  //cout << "----lastCandidate[" << lastCandidate << "]"  << endl;
                i--;  // on premet au channel courant d'etre ajoute au dt candidat suivant
              }
            } 
            
          } // end loop MTChannel
        }
      }
      cout << "Done  entry " << evtNum + 1 << " / " << nbEntries << "\r" <<  endl; // flush;;
    }


    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." <<  endl;
  cout << "nb skipped events "  << skipCandidate <<  endl;
  cout << "nb hits "  << nbHits <<  endl;
}



// Construct new CaloEvent from MTEvent for dtCandidate
bool AddCaloHit(CaloEvent& caloEvent,MTChannel* channel)
{
  //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->GetBcIdDif() - channel->GetBcIdHit();

    //if ( dt > 0 && 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;
}

---