/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/root/Tag_tracks.cpp

Go to the documentation of this file.

/* @version $Revision: 1566 $ * @modifiedby $Author: jacquem $ * @lastmodified $Date: 2012-03-08 14:00:10 +0100 (Thu, 08 Mar 2012) $ */
#include "tools/Log.hh"
#include "tools/Toolbox.hh"

#include "root/MTRun.hh"
#include "root/MTChannel.hh"
#include "root/MTEvent.hh"
#include "root/MTTrack.hh"
#include "MicroException.hh"

#include "mysql/Mysql.hh"


#include "mysql/Mysql.hh"

#include "TFile.h"
#include "TTree.h"
#include "TKey.h"
#include "TBranch.h"
#include "TRandom.h"
#include "TProfile.h"
#include "TH2I.h"
#include "TF1.h"


#include <sstream>
#include <sys/types.h>
#include <unistd.h>
#include <stdio.h>
#include <stdlib.h>


int errsv = errno;


using namespace std;

/*
 Modifying a existing TTree is often not possible because the tree is in a read-only
 file and you do not have permission to save the modified tree .
 Even if you do have the permission, you risk losing the
 original tree with an unsuccessful attempt to save  the modification.
*/



void Find_track(MTTrack *track,MTEvent& evt);

bool Valid_time(MTEvent *evt_prev, MTEvent *evt_curr);

void Fill_histos(int e,MTEvent* evt,TH1I *hdt,TH2I *hxy,TH2I *hxy_cut,TH1I *hnhit,TH1I *hnhit_cut,TH1I *hc_nhit,TH1I *hc_nhit_cut,TH1I *hnhit_ro,TH1I *hnhit_cut_ro,TH2F *hx_cor,TH2F *hy_cor,TH1F *hx_res,TH1F *hy_res,int chb);

void Calculate_offset(TH2F *hx_cor,TH2F *hy_cor,TH1F *hx_res,TH1F *hy_res,float* offset_x,float* offset_y,float* offset_x_err,float* offset_y_err);

void Average_pres_temp(TTree *tree,float* pres,float* temp,float* pres_rms,float* temp_rms,int npoints);



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

   if ( argc != 2  ) {
       FILE_LOG(logERROR)  << "usage:" << endl;
       FILE_LOG(logERROR)  << ":process_example fileToProcess" << endl;
       exit(1);
   }
  string inputfile = argv[1] ;

  // build temporaly fill path and name
  pid_t pid = getpid(); 
  std::ostringstream oss;
  oss << "/tmp/" << pid << ".root";
  string outputfile = oss.str();

  cout<<"Read   "<<inputfile<<endl;
  cout<<"Create "<<outputfile<<endl;



  // Open Root TFile in read only mode
  TFile infile(inputfile.c_str(),"READ");

  // Open Root TFile in write mode
  TFile oufile(outputfile.c_str(),"RECREATE");
  // Change current directory to "this" directory
  // Mandatory to use updateTree->GetCurrentFile() later
  oufile.cd();



  TH1I * hdt = new TH1I("hdt","Time to readout;bcid_dif - bcid_hit (200 ns)",100,0,100);

  TH2I * hxy = new TH2I("hxy","Hit position;y (pad number);x (pad number)",96,0,96,96,0,96);
  TH2I * hxy_cut = new TH2I("hxy_cut","Hit position after time cut;y (pad number);x (pad number)",96,0,96,96,0,96);

  TH1I * hnhit = new TH1I("hnhit","Number of hits;nhit",10000,0,10000);
  TH1I * hnhit_cut = new TH1I("hnhit_cut","Number of hits after time cut;nhit",10000,0,10000);

  TH1I * hc_nhit = new TH1I("hc_nhit","Number of hits;channel;nhit",9216,0,9216);
  TH1I * hc_nhit_cut = new TH1I("hc_nhit_cut","Number of hits after time cut;channel;nhit",9216,0,9216);

  TH2F * hx_cor = new TH2F("hx_cor","Correlation Prototype - telescope;Xtelescope (cm);Xm2 (cm)",100,0,100,200,0,100);
  TH2F * hy_cor = new TH2F("hy_cor","Correlation Prototype - telescope;Ytelescope (cm);Ym2 (cm)",100,0,100,200,0,100);

  TH1F * hx_res = new TH1F("hx_res","Correlation Prototype - telescope;Xm2 - Xtelescope (cm)",200,-100,100);
  TH1F * hy_res = new TH1F("hy_res","Correlation Prototype - telescope;Ym2 - Ytelescope (cm)",200,-100,100);



  //square events cuts
  UInt_t square_asu = 0;
  UInt_t square_bcid = 0;
  UInt_t square_cut = 1000;
  UInt_t m2_chb = 10;





  // Loop over original TTrees
  TIter nextkey(infile.GetListOfKeys());
  MTEvent *updateEvt = new MTEvent();

  TKey* key = (TKey*)nextkey();

  TTree *tree = (TTree*)key->ReadObj();

  // Create new TTree which will contain processed MTRun and MTEvent
  string treeName = tree->GetName();
  treeName = treeName + "_preprocessed";
  TTree *updateTree = new TTree(treeName.c_str(),"MicroMegas event");
  updateTree->Branch("MTEvent",&updateEvt,16000,2);

  MTRun * run  = (MTRun*)tree->GetUserInfo()->FindObject("MTRun");

  // Set Processed flag to MTRun
  run->SetProcessed(true);

  // Save MTRun information to new processed TTree
  updateTree->GetUserInfo()->Add(run);
  MTEvent *evt =  new MTEvent();
  TBranch *branch= tree->GetBranch("MTEvent");
  branch->SetAddress(&evt);


  //MTChannel* channel =NULL;
  int nbEntries = tree->GetEntries();
  MTEvent* prev_evt = NULL;

  cout<<endl;
  cout<<"NUMBER OF EVENTS = "<<nbEntries<<endl;
  cout<<endl;


  TH1I * hnhit_ro = new TH1I("hnhit_ro","Number of hits per readout;readout;nhit",nbEntries,0,nbEntries);
  TH1I * hnhit_cut_ro = new TH1I("hnhit_cut_ro","Number of hits per readout after time cut;readout;nhit",nbEntries,0,nbEntries);


  //for ( int evtNum = 0; evtNum < nbEntries ; evtNum++)
    for ( int evtNum = 0; evtNum < 500 ; evtNum++)
    {
      tree->GetEntry(evtNum);

      // Copy original MTEvent to new processed MTEvent
      *updateEvt = *evt;

      //check if time info is valid
      if ( prev_evt != NULL )
        {
          bool temp = Valid_time(prev_evt,updateEvt);
        }

      //look for square event in m2 prototype
      bool square = updateEvt->IsSquare(square_cut,square_asu,square_bcid,m2_chb);
      
      //Associate or not a track to this event
      //MTTrack* track = NULL;
      MTTrack* track = new MTTrack();
      Find_track(track,*updateEvt);

      //Write MTEvent into Branch
      cout  << "Write Event[" << evtNum << "] \r" << flush;
      updateTree->Fill();
      delete track;

      //Start basic analysis
      Fill_histos(evtNum,updateEvt,hdt,hxy,hxy_cut,hnhit,hnhit_cut,hc_nhit,hc_nhit_cut,hnhit_ro,hnhit_cut_ro,hx_cor,hy_cor,hx_res,hy_res,m2_chb);

      //keep event for time-between-readout calculation
      prev_evt = updateEvt;
    }



  //calculate offsets and errors
  float offset_x = 0;
  float offset_y = 0;

  float offset_x_err = 0;
  float offset_y_err = 0;

  Calculate_offset(hx_cor,hy_cor,hx_res,hy_res,&offset_x,&offset_y,&offset_x_err,&offset_y_err);
  


  //determine pressure and temperature
  float pres = 0;
  float temp = 0;
  float pres_rms = 0;
  float temp_rms = 0;
  int npoints = 10;
  Average_pres_temp(tree,&pres,&temp,&pres_rms,&temp_rms,npoints);


  //display profile infos

  float xmean = hxy_cut->GetMean(2);
  float ymean = hxy_cut->GetMean(1);
  float xrms = hxy_cut->GetRMS(2);
  float yrms = hxy_cut->GetRMS(1);
  float nhit_mean = hnhit->GetMean();
  float nhit_cut_mean = hnhit_cut->GetMean();
  int dt_max = hdt->GetMaximumBin()-1;
  float eff = 100. * (1 - 1.*hnhit_cut->GetBinContent(1)/hnhit_cut->GetEntries());

  cout<<"    Hit profile and RMS"<<endl;
  cout<<"Along X = ["<<xmean<<" +/- "<<xrms<<"] cm"<<endl;
  cout<<"Along Y = ["<<ymean<<" +/- "<<yrms<<"] cm"<<endl;
  cout<<endl;

  cout<<"    Nhits / readout"<<endl;
  cout<<"No time cut   = "<<nhit_mean<<endl;
  cout<<"With time cut = "<<nhit_cut_mean<<endl;
  cout<<endl;

  cout<<"    Timing"<<endl;
  cout<<"Delta_t_max = "<<dt_max<<" bcid"<<endl;
  cout<<endl;

  cout<<"    Efficiency = "<<eff<<" %"<<endl;
  cout<<endl;







  //Write offsets and other infos to database
  int runid = run->GetRunId();
  if (runid == 0)
    {
      cout<<"RunId not known, please enter: ";
      cin>>runid;
      cout<<endl;
      cout<<endl;
      //runid = 10227;
    }



  Mysql mysql;
  mysql.updateRunProcessInfo(runid, // runId
                             offset_x, // offsetX
                             offset_y, // offsetY 
                             pres, // pressure
                             temp, // temperature
                             dt_max);// deltaTMax 

  mysql.addRunChamberProcessInfo(runid,//runid
                                 m2_chb,  // chamberId,
                                 xmean, // xmean, 
                                 ymean, // ymean,
                                 eff,// inefficiency, 
                                 nhit_cut_mean,// nHitMean,
                                 dt_max);// deltaTMax;



  // Write Root Object to TFile
  TFile * curFile = updateTree->GetCurrentFile();
  updateTree->Write("", TObject::kOverwrite);
  updateTree->SetDirectory(0);

  // free heap memory
  delete evt;


  //Write histo to new TFile

  hdt->Write();
  hxy->Write();
  hxy_cut->Write();
  hc_nhit->Write();
  hc_nhit_cut->Write();
  hnhit->Write();
  hnhit_cut->Write();
  hnhit_ro->Write();
  hnhit_cut_ro->Write();
  hx_cor->Write();
  hy_cor->Write();
  hx_res->Write();
  hy_res->Write();



  infile.Close();
  oufile.Close();



  // Replace original root file by the new processed one
  string command = "cp -f ";
  command = command + outputfile + " "  + inputfile ;
  int result = system(command.c_str());
  if ( result == -1 ) { cout <<  strerror(errno) << endl; }     
  command = "rm -f " + outputfile;
  result = system(command.c_str());
  if ( result == -1 ) { cout <<  strerror(errno) << endl; }     


  cout << endl << "Done." << endl;


return 0;
}





/************************************************************/


void Average_pres_temp(TTree *tree,float* pres,float* temp,float* pres_rms,float* temp_rms,int npoints)

{

  int nbEntries = tree->GetEntries();

  MTEvent *evt =  new MTEvent();
  TBranch *branch= tree->GetBranch("MTEvent");
  branch->SetAddress(&evt);

  cout<<"    Average Pressure and Temperature ("<<npoints<<" points)"<<endl;

  int evtStep = floor(nbEntries/1./(npoints-1));

  float t = 0;
  float p = 0;

  for ( int evtNum = 0; evtNum < nbEntries ; evtNum += evtStep)
    {
      tree->GetEntry(evtNum);

      p = evt->GetPressure();
      t = evt->GetTemperature();

      *pres += p;
      *temp += t;

      *pres_rms += pow(p,2);
      *temp_rms += pow(t,2);
    }

  *pres /= npoints;
  *temp /= npoints;

  *pres_rms /= npoints;
  *temp_rms /= npoints;

  *pres_rms = sqrt(*pres_rms - pow(*pres,2));
  *temp_rms = sqrt(*temp_rms - pow(*temp,2));
  
  cout<<"P = ["<<*pres<<" +/- "<<*pres_rms<<"] mbar"<<endl;
  cout<<"T = ["<<*temp<<" +/- "<<*temp_rms<<"] deg C"<<endl;
  cout<<endl;

  delete evt;
  delete branch;

}










/************************************************************/



void Calculate_offset(TH2F *hx_cor,TH2F *hy_cor,TH1F *hx_res,TH1F *hy_res,float *offset_x,float *offset_y,float *offset_x_err,float *offset_y_err)
{

  TF1 * fx = new TF1("fx","[0]+[1]*x",0,100);
  TF1 * fy = new TF1("fy","[0]+[1]*x",0,100);

  fx->SetParameters(hx_res->GetMaximumBin() - 100,1);
  fy->SetParameters(hy_res->GetMaximumBin() - 100,1);

  TProfile * tpfx = hx_cor->ProfileX();
  TProfile * tpfy = hy_cor->ProfileX();

  tpfx->Fit(fx,"q"); 
  tpfy->Fit(fy,"q"); 

  *offset_x = fx->GetParameter(0);
  *offset_x_err = fx->GetParError(0);

  *offset_y = fy->GetParameter(0);
  *offset_y_err = fy->GetParError(0);

  tpfx->Delete();
  tpfy->Delete();

  fx->Delete();
  fy->Delete();

  cout<<endl;
  cout<<"    Offsets"<<endl;
  cout<<"Along X = ["<<*offset_x<<" +/- "<<*offset_x_err<<"] cm"<<endl;
  cout<<"Along Y = ["<<*offset_y<<" +/- "<<*offset_y_err<<"] cm"<<endl;
  cout<<endl;

}








/************************************************************/


void Fill_histos(int e,MTEvent* evt,TH1I *hdt,TH2I *hxy,TH2I *hxy_cut,TH1I *hnhit,TH1I *hnhit_cut,TH1I *hc_nhit,TH1I *hc_nhit_cut,TH1I *hnhit_ro,TH1I *hnhit_cut_ro,TH2F *hx_cor,TH2F *hy_cor,TH1F *hx_res,TH1F *hy_res,int chb)
{



  //look for a straight track to fill correlation histos
  MTTrack * trk = evt->GetTrack();
  bool correlation = 0;
  float bx = 0;
  float by = 0;

  if (trk!=NULL)
    {
      if (trk->GetFromPad() && trk->GetStraight())
        {
          correlation = 1;
          bx = trk->GetBx();
          by = trk->GetBy();
        }
    }

  //delete trk;


  int nhit = 0;
  int nhit_cut = 0;
  int nchannel = evt->GetNchannel();

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



      if (channel->GetChamberId() == chb)
        {
          int hardid = channel->GetHardId();
          int chipid = channel->GetChipId();

          int channelid = (chipid-1)*64 + hardid;

          float xpos = (int)(channel->GetX() / 10000);
          float ypos = (int)(channel->GetY() / 10000);

          int digit = channel->GetDigitalValue();

          hc_nhit->Fill(channelid);
          hxy->Fill(ypos,xpos);   
          nhit++;

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

          hdt->Fill(dt);

          if (dt>=6 && dt<=9)
            {
              hc_nhit_cut->Fill(channelid);
              hxy_cut->Fill(ypos,xpos);   
              nhit_cut++;

              if (correlation)
                {
                  //take m2 prototype dead area into account
                  //parallel to DIF side
                  if (xpos >= 32)
                    {
                      xpos += 0.5;

                      if (xpos >= 32*2)
                        {
                          xpos += 0.5;
                        }
                    }

                  //perpendicular to DIF side
                  if (ypos >= 48)
                    {
                      ypos += 0.5;
                    }

                  hx_res->Fill(xpos-bx);
                  hy_res->Fill(ypos-by);

                  hx_cor->Fill(bx,xpos);
                  hy_cor->Fill(by,ypos);
                }
            }
        }
    }

  hnhit->Fill(nhit);
  hnhit_cut->Fill(nhit_cut);

  hnhit_ro->Fill(e+1,nhit);
  hnhit_cut_ro->Fill(e+1,nhit_cut);

}








/***************************************************/






bool Valid_time(MTEvent *evt_prev, MTEvent* evt_curr)
{

  bool valid = 0;



  Long64_t tprevious = 0;
  Long64_t tcurrent  = 0;

  bool tprev_defined = 0;
  bool tcurr_defined = 0;
 
  // search time in previous event
  int nchannel = evt_prev->GetNchannel();
  for (int j=0;j<nchannel;j++)
    {             
      MTChannel* channel = (MTChannel*)evt_prev->GetChannels()->UncheckedAt(j);
      
      if (channel->GetChamberId() == 10)
        {
          tprevious = channel->GetBcId_Abs();
          tprev_defined = 1;
          j = nchannel;
        }
    }
  //if no m2 hits were found, no time info
  if (tprev_defined == 0){return 0;}

  // search time in current event
  nchannel = evt_curr->GetNchannel();
  for (int j=0;j<nchannel;j++)
    {             
      MTChannel* channel = (MTChannel*)evt_curr->GetChannels()->UncheckedAt(j);
      
      if (channel->GetChamberId() == 10)
        {
          tcurrent = channel->GetBcId_Abs();
          tcurr_defined = 1;
          j = nchannel;
        }
    }
  //if no m2 hits were found, no time info
  if (tcurr_defined == 0){return 0;}

  //calculate time difference

  if (tcurr_defined && tprev_defined)
    {
      Long64_t tdiff = tcurrent - tprevious;
      double dt = tdiff * 200e-9;

      if (dt <= 3.35)
        {
          valid = 1;
          evt_curr->SetTimeInfo(1);
        }
    }

  return valid;
}





/****************************************/




void Find_track(MTTrack *trk,MTEvent& evt)
{

  int adc_cut = 25;
  int tel_square_cut = 10;

  bool display = 0;

  TH2I * hxz = new TH2I("hxz","Hit pattern in XZ plane;z;x",40,0,40,16,0,16);
  TH2I * hyz = new TH2I("hyz","Hit pattern in YZ plane;z;y",40,0,40,16,0,16);
  TH2I * hxy[3];
  for (int i=0;i<3;i++){hxy[i] = new TH2I(Form("hxy_%i",i+7),Form("Hit pattern (adc>%i) chambeDr %i;y (cm);x (cm)",adc_cut,i+7),16,0,16,16,0,16);}

  float zmean[3] = {10,20,30};
                  
  float xavg = 0;
  float yavg = 0;
  float zavg = 20;
  float ax = 0;
  float bx = 0;
  float ay = 0;
  float by = 0;
  float ax1 = 0;
  float ax2 = 0;
  float ay1 = 0;
  float ay2 = 0;
  
  float xextrapol = 0;
  float yextrapol = 0;

  int content = 0;
  int nhit[3] = {0,0,0};
  int xmax[3] = {0,0,0};
  int ymax[3] = {0,0,0};
  float xmean[3] = {0,0,0};
  float ymean[3] = {0,0,0};
  int ztemp;

  bool track = 0;
  bool hot = 0;
  bool square_tel = 0;
  bool xstraight = 0;
  bool ystraight = 0;
  bool straight = 0;
  bool stiff = 0;




  //noisy channels
  bool noisy[3][96];
  for (int i=0;i<3;i++){
    for (int j=0;j<96;j++){
      noisy[i][j] = 0;}}
  
  /*HV pad*/
  noisy[0][77] = 1;
  noisy[1][77] = 1;
  noisy[2][77] = 1;
   
  /*floating pad*/
  noisy[0][9] = 1;
  noisy[2][33] = 1;

  int nhit_tel_cut[3] = {0,0,0};
  bool hot_tel[3] = {0,0,0};


  //fill XZ and YZ histos

  if (display){cout<<evt.GetNchannel()<<endl;}

  for (int j=0;j<evt.GetNchannel();j++)
    {
      MTChannel* channel = (MTChannel*)evt.GetChannels()->UncheckedAt(j);
      int chb = channel->GetChamberId() - 7;
      bool pad_telescope = chb >= 0 && chb < 3;

      if (pad_telescope)
        {
          int hardid = channel->GetHardId();

          if (noisy[chb][hardid] == 0)
            {
              int adc = channel->GetAnalogValue();

              if (adc >= adc_cut)
                {
                  nhit_tel_cut[chb]++;

                  int xpos = (int)(channel->GetX() / 10000);
                  int ypos = (int)(channel->GetY() / 10000);
                  int zpos = (int)(channel->GetZ() / 10000);

                  if (!hot_tel[chb]){hot_tel[chb] = 1;}

                  hxy[chb]->Fill(xpos,ypos,adc);
                  hxz->Fill((chb+1)*10,xpos,adc);
                  hyz->Fill((chb+1)*10,ypos,adc);
                }
            }
        }
    }
  
  //require 3 hot chambers
  track = 0;
  hot = hot_tel[0] && hot_tel[1] && hot_tel[2];

  if (hot)
    {
      //tag square events    
      square_tel = 0;    
      for (int k=0;k<3;k++)
        {
          if (nhit_tel_cut[k] >= tel_square_cut)
            {
              square_tel = 1;
            }
        }

      //find tracks
      if (!square_tel)
        {
          //get position with maximum charge
          for (int k=0;k<3;k++)
            {
              hxy[k]->GetMaximumBin(xmax[k],ymax[k],ztemp);
            }
          
          //request align maxima +/- 1 pad
          //define track from first and last chamber
          
          xextrapol = floor(0.5*(xmax[2] + xmax[0]));
          yextrapol = floor(0.5*(ymax[2] + ymax[0]));
          
          if ((fabs(xextrapol - xmax[1]) <= 1) && (fabs(yextrapol - ymax[1]) <= 1))
            {
              track = 1;
            }
          
          if (track)
            {
              /*************form clusters**************/      
                  
              for (int k=0;k<3;k++)
                {
                  if (nhit_tel_cut[k] != 1)
                    {
                      for (int xx=-1;xx<=1;xx++)
                        {
                          for (int yy=-1;yy<=1;yy++)
                            {
                              content = hxy[k]->GetBinContent(xmax[k] + xx,ymax[k] + yy);
                                    
                              if (content != 0)
                                {
                                  xmean[k] += (xmax[k] + xx - 0.5);
                                  ymean[k] += (ymax[k] + yy - 0.5);
                                  nhit[k]++;
                                }
                            }
                        }

                      if (nhit[k] != 1)
                        {
                          xmean[k] /= nhit[k];
                          ymean[k] /= nhit[k];
                        }
                    }
                  else
                    {
                      xmean[k] = xmax[k] - 0.5;
                      ymean[k] = ymax[k] - 0.5;
                    }
                }

              /************determine track parameters**********/
              
              //tracks in xz plane
              //if straight
              if (xmean[0]==xmean[1] && xmean[0]==xmean[2])
                {
                  ax = 0;
                  bx = xmean[0];
                  xstraight = 1;
                }
              //if inclined
              else
                {
                  xavg = 1/3.*(xmean[0] + xmean[1] + xmean[2]);
                    
                  for (int l=0;l<3;l++)
                    {
                      ax1 += ((zmean[l] - zavg) * (xmean[l] - xavg));
                      ax2 += pow(zmean[l] - zavg,2);
                    }
                      
                  ax = ax1/ax2;
                  bx = xavg - ax*zavg;
                }
                  
              //tracks in yz plane
              //if straight
              if (ymean[0]==ymean[1] && ymean[0]==ymean[2])
                {
                  ay = 0;
                  by = ymean[0];
                  ystraight = 1;
                }
              //if inclined
              else
                {
                  yavg = 1/3.*(ymean[0] + ymean[1] + ymean[2]);
                      
                  for (int l=0;l<3;l++)
                    {
                      ay1 += ((zmean[l] - zavg) * (ymean[l] - yavg));
                      ay2 += pow(zmean[l] - zavg,2);
                    }
                    
                  ay = ay1/ay2;
                  by = yavg - ay*zavg;
                }
              
              straight = xstraight && ystraight;
              stiff = straight && nhit_tel_cut[0]==1 && nhit_tel_cut[1]==1 && nhit_tel_cut[2]==1;

              trk->SetFromPad(1);
              trk->SetFromStrip(0);

              trk->SetStraight(straight);
              trk->SetSingleHit(stiff);

              trk->SetAx(ax);
              trk->SetAy(ay);

              trk->SetBx(bx);
              trk->SetBy(by);

              evt.SetTrack(trk);

            }
        }
    }


  if (display)
    {
      cout<<"Hot = "<<hot<<endl;
      if (hot)
        {
          cout<<"Track = "<<track<<endl;
          if (track)
            {
              cout<<"XZ plane : "<<ax<<" "<<bx<<endl;
              cout<<"YZ plane : "<<ay<<" "<<by<<endl;
            }
        }
    }


  hxz->Delete();
  hyz->Delete();

  for (int k=0;k<3;k++)
    {
      hxy[k]->Delete();
    }


}

---