/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/lcio/max_lcio.cc

Go to the documentation of this file.

#include "lcio.h"

#include "IO/LCWriter.h"
#include "EVENT/LCIO.h"
#include <EVENT/CalorimeterHit.h>
#include "EVENT/LCGenericObject.h"
#include "IMPL/LCCollectionVec.h"
#include "IMPL/CalorimeterHitImpl.h"

#include "TNtuple.h"
#include "TH1F.h"
#include "TFile.h"
#include "TGraph.h"


#include <cstdlib>
#include <iostream>
#include <sstream>


using namespace std ;
using namespace lcio ;



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


  if (argc==3){

  static string lcioFile = argv[1] ;
  TString rootfilename = argv[2];

  IO::LCReader* lcReader = IOIMPL::LCFactory::getInstance()->createLCReader() ;
  lcReader->open( lcioFile ) ;

  bool display = 0;

  int n = 1e9;
  int i = 0;

  LCEvent* evt = NULL;

  int nthEvent = 500;
  lcReader->skipNEvents(nthEvent) ;

  TFile * tf = new TFile(rootfilename.Data(),"RECREATE");

  TNtuple * nthit = new TNtuple("nthit","AHCAL Calo Hits","x:y:z:e");
  TNtuple * ntdif = new TNtuple("ntdif","AHCAL synchro DIF trigger number","trg");

  TGraph * tg_dif_trg_nevent = new TGraph();
  tg_dif_trg_nevent->SetName("tg_dif_trg_nevent");
  tg_dif_trg_nevent->SetTitle("Synchro DIF AHCAL trigger number");
  tg_dif_trg_nevent->GetXaxis()->SetTitle("AHCAL event number");
  tg_dif_trg_nevent->GetYaxis()->SetTitle("DIF trigger number");

  TH1F * hen = new TH1F("hen","Calo Hit energy;energy (GeV)",500,0,5);
  TH1F * hx = new TH1F("hx","Calo Hit position X;x (mm)",400,-600,600);
  TH1F * hy = new TH1F("hy","Calo Hit position Y;y (mm)",400,-600,600);
  TH1F * hz = new TH1F("hz","Calo Hit position Z;z (mm)",1200,0,1200);


  float energy = 0;
  float xpos = 0;
  float ypos = 0;
  float zpos = 0;
  int dif_trg = 0;

  TString collec_name;
  TString collec_type;

  cout<<endl;
  cout<<"Read file "<<lcioFile<<endl;
  cout<<endl;


 
  MTEvent *ahcalEvt =  new MTEvent();
  TTree mt(lcioFile,"AHCAL MT format");
  mt.SetMaxTreeSize(4000000000); // 4 Mb
  mt.Branch("MTEvent",&ahaclEvt,16000,2);// second argument is class of evt!!!






  while( (evt = lcReader->readNextEvent()) != 0 ){
      
    i++;
    // Fill MT Event
    ahcalEvt->SetGlobalTriggerCounter(i);

    if (i%100==0){cout<<"Read event "<<i<<" \r"<<flush;}

    //LCTOOLS::dumpEventDetailed( evt );

    const vector<string>* collNames = evt->getCollectionNames();
        
    vector<string>::const_iterator it;
        
    for ( it = collNames->begin(); it!=collNames->end(); ++it) {

      LCCollection*  collec = evt->getCollection(*it);
      if (display){cout<<endl;cout << " collection "<<*it<<" of type " << collec->getTypeName() << endl;}

      collec_name = *it;
      collec_type = collec->getTypeName();


      /*Calo Hits*/

      if (collec_type == "CalorimeterHit"){

              for ( int element = 0; element < collec->getNumberOfElements(); element++){

                CalorimeterHit *hit = dynamic_cast<CalorimeterHit*>(collec->getElementAt(element));

          if (hit != NULL){

                  if(display){cout<<"Hit with energy "<<hit->getEnergy()<<" GeV at position "<<hit->getPosition()[0]<<" "<<hit->getPosition()[1]<<" "<<hit->getPosition()[2]<<endl;}

                  energy = hit->getEnergy();
                  xpos = hit->getPosition()[0];
                  ypos = hit->getPosition()[1];
                  zpos = hit->getPosition()[2];
       
                  nthit->Fill(xpos,ypos,zpos,energy);
       
       
                  hen->Fill(energy);
                  hx->Fill(xpos);
                  hy->Fill(ypos);
                  hz->Fill(zpos);}

            // Fill MT Event Hits
            MTChannel mtHit();
            mtHit.SetX(xpos);
            mtHit.SetY(ypos);
            mtHit.SetZ(zpos);
            mtHit.SetAnalogValue( int(energy));
            ahcalEvt->AddChannel(mtHit);
            
      }
   } 


      /*DIF trigger*/

      if (collec_name == "CALDAQ2_DifTrigger"){

        for ( int element = 0; element < collec->getNumberOfElements(); element++){

          LCGenericObject *obj = dynamic_cast<LCGenericObject*>(collec->getElementAt(element));

          if (obj != NULL){     

            dif_trg = obj->getIntVal(0);

            ntdif->Fill(dif_trg);
            tg_dif_trg_nevent->SetPoint(tg_dif_trg_nevent->GetN(),i,dif_trg);

      // Fill MT Event dif synchro
      ahcalEvt->SetDifSynchro(dif_trg);

            if (display){cout<<"  dif trigger number = "<<dif_trg<<endl;}}}}}


    

    if (i==n){break;}
    
    mt->Fill();  
    ahcalEvt->SetNchannel(0); 
  } // end while event



  lcReader->close() ;
  delete lcReader ;



  tg_dif_trg_nevent->Write();


  tf->Write();
  tf->Close();}


  else{cout<<"argument 1 is the slcio file - arg 2 is the final rootfile with position, energy and trigger number"<<endl;}


  return 0;}

---