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

Go to the documentation of this file.

/* @version $Revision: 1376 $ * @modifiedby $Author: jacquem $ * @lastmodified $Date: 2011-11-18 16:11:02 +0100 (Fri, 18 Nov 2011) $ */
// This program extracts the nb of hits for Muons and Pions  to compare with MC and do linearity or resolution
//Yannis 17 Aug 2012
#include <map>
#include "TKey.h"
#include <TH2F.h>
#include <TProfile.h>
#include <TROOT.h>
#include <TRint.h>
#include <TKey.h>
#include <TList.h>
#include <TTree.h>
#include <TFile.h>
#include <TChain.h>
#include "TChainElement.h"
#include <TMath.h>
#include <stdlib.h>
#include <TSystem.h>
#include <iostream>
#include <sstream>
#include "Log.hh"


#include "MicroException.hh"
#include <string>

#include "mysql/Mysql.hh"


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

#include "TGraphErrors.h"
#include "TGraphAsymmErrors.h"
#include <TF1.h>
#include "tools/Arguments.hh"

using namespace std;
const int maxEvents = 200000; // How many events to read
const int maxTrees = 100;  // How many trees to read
const int maxHits = 5000; // Max nb of considered hits
const float AverageMult = 1.7; // Average multiplicity per track. Nb of Hits is normalized to this value

int main(int argc, char**argv){
  Arguments arg(argc,argv," getNbOfHits run1 [run2] max [run20] -p particleType -r opt");

  if ( arg.getNbOfArgs() < 1) { 
    arg.usage();
    exit(1);
  }

  // Analyse the command line to find what to do
  
  int NbOfRuns = arg.getNbOfArgs();

  string rootName[20];
  for (int nr=0 ; nr < NbOfRuns ; nr ++ ){rootName[nr].assign(arg.getArgument(nr+1)); }

  string  ParticleType = "Pions";
  if ( arg.getOption("-p").size() != 0 )
    {
      ParticleType  = arg.getOption("-p");
    }

  string  RootFileType = "EM";  // With Efficiency and Multiplicity Corrections
  if ( arg.getOption("-r").size() != 0 )
    {
      RootFileType  = arg.getOption("-r");
    }

  TString Period = "SPS_aug_2012";
  if ( arg.getOption("-t").size() != 0 )
    {
      Period  = arg.getOption("-t");
    }

  cout << " Starting Analysis Run " << rootName[0] << endl;
  int search = 1;

  int NbHitLayer[49];
  int IsPenetratingMIP;

  TH1I* hnhit_roi = new TH1I("hnhit_roi","Number of hits in tight DeltaT",200,0,2000.);
  TH1F* hDensity = new TH1F("hDensity","Hits Density ",200,0,50.);
  TH1I* hHitsTot48 = new TH1I("hHitsTot48", "Nb of Hits Thresh 1 Chamber 48", 200, 0, 200);  
  TH1I* hHitsTot49 = new TH1I("hHitsTot49", "Nb of Hits Thresh 1 Chamber 49", 200, 0, 200);  
  TH1I * hnHitForw = new TH1I("hnhitForw","MIP Number of hits Forw",50,0,50);
  TH2I* hxy = new TH2I("hxy","Hit Map on first layer",96, 0, 96, 96, 0, 96);  
  TH2I* hxyOrg = new TH2I("hxyOrg","Hit Map on Central Pad",96, 0, 96, 96, 0, 96);  
  TH2I* hPadCentr = new TH2I("hPadCentr","Central Pad",96, 0, 96, 96, 0, 96);  
  TH2F* hXYExtr = new TH2F("hXYExtr","XY Extrapolation",96, 0, 96, 96, 0, 96); 
  TH2F*   hS3S9 = new TH2F("hS3S9","NB of Hits in 3x3/10x10 Nhits", 400,0.,2000., 100, 0., 1.); 
  TProfile * hRegMult = new TProfile("hRegMult", "Nb of Hits in n+1xn+1 cumulative matrice", 50, 0, 50., 0.9, 1000.);
  TProfile * hRegMultC = new TProfile("hRegMultC", "Nb of Hits Total-center matri", 50, 0, 50., 0., 1000.);      

  TH2F*   hRMSx = new TH2F("hRMSx","NB of Hits VS RMSX", 400,0.,2000.,100,0,25);
  TH2F*   hRMSy = new TH2F("hRMSy","NB of Hits VS RMSY", 400,0.,2000.,100,0,25);
  TH2I* hNhitDeb = new TH2I("hNhitDeb", " Deb", 400,0.,2000.,400,0.,2000.);
  TH1F* hNhitDebrap = new TH1F("hNhitDebrap", " Debaa", 100, 0., 5.);
  TH2F* hNhitRapHits = new TH2F("hNhitRapHits", " Debaa", 400,0.,2000.,100, 0., 5.);
  TH1F*  hRG1G2 = new TH1F("hRG1G2", "G1/G2", 100, 0., 10.); 
  TH2F*  hnHitsRG1G2 = new TH2F("hnHitsRG1G2", "Nhits vs G1/G3", 400,0.,2000.,100, 0., 10.);
  //  if(rootName[0]=="715615" ||  rootName[0]=="715692" || rootName[0]=="715699" || rootName[0]=="715747"){xhlim = 700;}
  TH1I* hNbHitsPions = new TH1I("hhNbHitsPions","Number of hits for Pions",520,0,2600.);
  TH1I* hNbHitsPions1 = new TH1I("hhNbHitsPions1","Number of hits for Pions",520,0,2600.);
  TH1I* hNbHitsPions2 = new TH1I("hhNbHitsPions2","Number of hits for Pions",200,0,1000.);
  TH1I* hNbHitsPions3 = new TH1I("hhNbHitsPions3","Number of hits for Pions",100,0,500.);
  TH1I* hNbHitsPionsSC = new TH1I("hhNbHitsPionsSC","Number of hits for Pions",520,0,2600.);
  TProfile * hNoiseLay = new TProfile("hNoiseLay","Noise Per layer",50, 0., 50., 0., 100.,"S");
  TH1F* hShCont1 = new TH1F("hShCont1","Rear to Forward Nb hits ratio", 100, 0., 1.);
  TH1F* hNbHitsRear = new TH1F("hNbHitsRear","Nb Hits Rear", 100, 0., 100.);
 TH1F* hNbHitsBoard = new TH1F("hNbHitsBoard","Nb Hits Board", 200, 0., 500.);

  TH2I * hxyMap[50]; 
  TH2I * hxyMap1[50];
  TH2I * hxyMap2[50];
  TH2I * hxyMap3[50];

  TString name,title;
  int NbinsEffx=12;
  int NbinsEffy=12;
  int hbin=96;
  for (int i=0;i<50;i++)
    {
      name = Form("hxyMap_%i",i);
      title = Form("Hit distribution in Layer %i",i);
      hxyMap[i]=  new TH2I(name,title,96,0,hbin,96,0,hbin);
      hxyMap1[i]= new TH2I(name+" Th1",title,96,0,hbin,96,0,hbin);
      hxyMap2[i]= new TH2I(name+" Th2",title,96,0,hbin,96,0,hbin);
      hxyMap3[i]= new TH2I(name+" Th3",title,96,0,hbin,96,0,hbin);
    }
  TH2I * hxyTot = new TH2I("hxyTot","",96,0,hbin,96,0,hbin);

  TString OutFilename = "/lapp_data/LC/Detecteurs/Sid/karyotak/DataTB2012/rootfiles/"+ParticleType+rootName[0]+RootFileType+".root";
  TFile tf(OutFilename,"RECREATE");


  cout << " Writting file : " << OutFilename << endl;
  
  int nbEvtTot=0;
  int NbTrees = 0;
  int RunId=1;
  int nhwrt=0;
  int nhit_roi_old=1;
  int nhit_roi;
  float RG1G3 = -100.;
  int maxDeltatOld = -100;
  int maxDeltat    = -100;
  UInt_t  BcIdAbsOld =0;
  UInt_t  BcIdAbs  =0;
  string Ebeam; double BeamEnergy ;
  double  nHitsTot; double  nHitsTot1;  double  nHitsTot2;       double  nHitsTot3;
  double ShCont1;  double NbHitsRear; double NbHitsBoard; int zstart;

  TString PionsOutputFileName="/lapp_data/LC/Detecteurs/Sid/karyotak/DataTB2012/rootfiles/DSTPions"+rootName[0]+RootFileType+".root";
  TFile* PionsRootFile = new TFile(PionsOutputFileName, "recreate");
  TTree* Ptree = new TTree("Pions","");
  Ptree->Branch("Ebeam",   &BeamEnergy,    "BeamEnergy/D");
  Ptree->Branch("Nhtot",   &nHitsTot,  "nHitsTot/D");
  Ptree->Branch("Nhtot1",  &nHitsTot1, "nHitsTot/D");
  Ptree->Branch("Nhtot2",  &nHitsTot2, "nHitsTot/D");
  Ptree->Branch("Nhtot3",  &nHitsTot3, "nHitsTot/D");
  Ptree->Branch("ShCont1", &ShCont1,   "ShCont1/D");
  Ptree->Branch("NbHitsRear", &NbHitsRear,   "NbHitsRear/D");
  Ptree->Branch("NbHitsBoard", &NbHitsBoard, "NbHitsBoard/D");
  Ptree->Branch("Zstart",  &zstart, "zstart/I");

  // Define Beam spot limits

  Int_t YPads[10] =  {52, 48, 63, 43, 63, 63, 44, 45, 45, 50};
  Int_t XPads[10] =  {58, 50, 57, 48, 58, 61, 45, 61, 55, 42};
  Int_t YPadsL[10] = {10, 10, 10, 10, 25, 25, 10, 10, 10, 5};
  Int_t YPadsH[10] = {10, 10, 10, 10, 10, 10, 10, 25, 10, 5};
  Int_t XPadsL[10] = {10, 10, 10, 10, 25, 25, 10, 25, 10, 5};
  Int_t XPadsH[10] = {10, 10, 10, 10, 10, 10, 10, 10, 10, 5};
  int jBeam = 9;
  TString Runs[9]={"714489","714525","714531","714546", "714553","714556","714559","714562","714573"};
  for (int r=0; r<9 && ParticleType=="Pions" ; r++){
    cout << r << " " << Runs[r] << " " << " Check run for Beam spot " << rootName[0] << endl;
    if( rootName[0]==Runs[r] ){jBeam = r; cout << "JBeam " << jBeam << " Beam spot at " << YPads[r] << " " << XPads[r] << endl; break; }
}
  cout << " Found Beam spot at " << YPads[jBeam] << " " << XPads[jBeam] << endl;
  
  TChain chain("T");
  int nbHit = 0;
  Double_t Zlayer[maxHits]; // Before it was Float
  Double_t Xlayer[maxHits];
  Double_t Ylayer[maxHits];
  Double_t ErPos[maxHits];

  cout << " Reading Files : " << endl;
  for ( int nr=0; nr< NbOfRuns; nr++){
    for (int i =0; i<5 ; i++){
      TString SubFile = Form("_I%i",i);
      TString fname = Form("/lapp_data/LC/data/TB2012/"+Period+"/Production/CaloHits/CaloHit_DHCAL_"+rootName[nr]+SubFile+"_0.slcio.root");
      chain.Add(fname);
      cout << fname << endl;
    }
  }
  Float_t EffMap[49]; Float_t MultMap[49];
  for (int m=0; m<=49 ; m++){EffMap[m]=1.; MultMap[m] = AverageMult;}
  if(RootFileType =="EM" || RootFileType == "E" ){ // Load files only if needed
    // Get efficiencies and Multiplicities

    TString EffFilename="/gpfs/LAPP-DATA/LC/Detecteurs/Sid/karyotak/DataTB2012/rootfiles/Muons"+rootName[0]+"Eff.root"; // Gets the efficiency of the 1st run 
    //    TString EffFilename="/lapp_data/LC/Detecteurs/Sid/karyotak/DataTB2012/rootfiles/MuonsEff"+rootName[0]+".root"; // Gets the efficiency of the 1st run 
    TFile ChEff(EffFilename);
    TH1F* hGlEffChamber = (TH1F*) ChEff.Get("hGlEffChamber");
    TProfile * hLayerMult = (TProfile*) ChEff.Get("hLayerMult");
    cout << "Reading efficiencies and Multiplicities " << endl;
    for (int l = 0; l<=49; l++){
      if(RootFileType=="E" || RootFileType=="EM"){EffMap[l] = 0.01*hGlEffChamber->GetBinContent(l+1);}
      if(RootFileType=="M" || RootFileType=="EM")MultMap[l]=hLayerMult->GetBinContent(l+1);
      cout << "Layer " << l << " Eff " << EffMap[l] << " Multiplicity " << MultMap[l] << endl;
    }
  }

  TObjArray *fileElements=chain.GetListOfFiles();
  TIter next(fileElements);
  TChainElement *chEl=0;
  int NchainFiles=0;
  while (( chEl=(TChainElement*)next() )) 
    {
      TFile f(chEl->GetTitle());
      //  TFile f(rootName.c_str());
      TIter nextkey(f.GetListOfKeys());
      TKey *key;
      NchainFiles++;

      int nbEvt = 0;

      cout << " Reading a new file in the chain " << NchainFiles << endl;
      while (key = (TKey*)nextkey()) 
        {

          TString keyname = key->GetName();
          if ((keyname == "hResidualX") || (keyname == "hResidualY")){}
          else
            {

              TTree *tree = (TTree*)key->ReadObj();                
              cout << "---------------------------- NEW TTREE -----------------"<< endl;
              nbHit = 0;
              NbTrees++;
              if( NbTrees > maxTrees ) break;
              CaloRun *run=NULL;

              try {
                CaloRun *run=dynamic_cast<CaloRun *>(tree->GetUserInfo()->First());
                if( run != NULL ){cout << "Run[" << run->GetRunId() <<  "], temperature[" << run->GetTemperature() <<  "] pression[" << run->GetPressure() << "] " << endl ; }

              
                // Get Beam energy form runId
                Mysql mysql;
                try {
                  Ebeam  = mysql.getRunEnergy(run->GetRunId());
                }
                catch (MicroException e) {}
                int le = Ebeam.find(" "); string EbeamN = Ebeam.substr(0,le);  BeamEnergy = strtod (Ebeam.c_str(),NULL); 
                cout << "Beam[" << Ebeam << "]  " << BeamEnergy << endl;
              }
              catch ( ... ) {}
              CaloEvent *evt =  new CaloEvent();
              TBranch *branch= tree->GetBranch("CaloEvent");
              branch->SetAddress(&evt);

              CaloHit* hit =NULL;
              int nbEntries = tree->GetEntries();
              cout << " Nb of events to read " << nbEntries << endl;
    
              for ( int evtNum = 0; evtNum < nbEntries ; evtNum++)
                {
                  if(nbEvt>=maxEvents) break;
                  nbEvt++; nbEvtTot++;

                  hxyTot->Reset();
                  for (int i=0; i<50; i++){hxyMap[i]->Reset();hxyMap1[i]->Reset();hxyMap2[i]->Reset();hxyMap3[i]->Reset();}
                  tree->GetEntry(evtNum);

                  int nbHits = evt->GetNHits();

                  //number of hits in region of interest
                  nHitsTot=0.;  nHitsTot1=0.;  nHitsTot2=0.;  nHitsTot3=0.;
                  int nhit_fit = 0;
                  float xpos = 0;
                  float ypos = 0;
                  float zpos = 0;
                  int ncol = 0;
                  int nrow = 0;
                  int zSlot = 0;
                  int deltat=0;
                  int Thresh = 0;
                  nhit_roi_old = nhit_roi;
                  maxDeltatOld = maxDeltat;
                  BcIdAbsOld = BcIdAbs;

                  nhit_roi = 0;
                  float Rad1lay=10000.;
                  int XOrg=1000; int YOrg=1000;
                  maxDeltat = evt->GetDeltaTmax();
                  BcIdAbs = evt->GetBcIdAbs();

                  for (int layer=0 ; layer <=49 ; layer++){
                    NbHitLayer[layer]=0;
                  }
                  for (int i=0 ; i<maxHits ;i++){
                    Zlayer[i]=0;
                    Xlayer[i]=0;
                    Ylayer[i]=0;
                    double theta = (0.0136/100.)*sqrt((i+1)*2/1.757)*(1+0.038*log((i+1)*2/1.757));
                    ErPos[i] = sqrt(pow(10./sqrt(12.),2)+pow(2*28.*theta,2));
                  }

                  int NbHitLayer33=0;
                  int NbHitLayer77=0;

                  for(int i=0 ; i < nbHits ; i++)
                    {
                      hit = (CaloHit*)evt->GetHits()->UncheckedAt(i);

                      zSlot = hit->GetLayer() ;
                      deltat=hit->GetDeltaT();
                      xpos =  hit->GetX()/1000. ; // unit = mm
                      ypos =  hit->GetY()/1000. ;
                      zpos =  hit->GetZ()/1000.;
                      ncol = hit->GetColInChamber();
                      nrow = hit->GetRowInChamber();
                      Thresh = hit->GetThreshold();
                      // cout << Thresh << " " << endl;
                  if(ParticleType=="Noise"){
                      if(zSlot <48 && abs(deltat-maxDeltat)<=2) continue; // This is for RPCs
                      if(zSlot >=48 && abs(deltat-maxDeltat)<=5) continue; // For Micromegas
                  }
                  else{
                      if(zSlot <48 && abs(deltat-maxDeltat)>2) continue; // This is for RPCs
                      if(zSlot >=48 && abs(deltat-maxDeltat)>5) continue; // For Micromegas
                  }

                      //                      if(zSlot >=48 ) continue; // Ignore  Micromegas


                      nhit_roi++; // Nb of good hits
                      hxyMap[zSlot]->Fill(ncol,nrow,1);
                      if(Thresh == 1) {hxyMap1[zSlot]->Fill(ncol,nrow,1);}
                      if(Thresh == 2) {hxyMap2[zSlot]->Fill(ncol,nrow,1);}
                      if(Thresh == 3) {hxyMap3[zSlot]->Fill(ncol,nrow,1);}
                      hxyTot->Fill(ncol,nrow);

                      if(zSlot==0){
                        //  hxy->Fill(ncol,nrow);
                        float rad = sqrt( (nrow-58.)*(nrow-58.)+(ncol-50.)*(ncol-50.) );
                        if(rad < Rad1lay) {XOrg=ncol; YOrg=nrow; Rad1lay=rad;}
                      }
                      //    cout << XOrg << " Y org" << YOrg << endl;
                      NbHitLayer[zSlot]++; 
                      if(abs(ncol-XOrg)<=2 && abs(nrow-YOrg)<=2){NbHitLayer33++;}
                      if(abs(ncol-XOrg)<=10 && abs(nrow-YOrg)<=10){NbHitLayer77++;}


                      // Prepare for fitting a muon
                      // Choose layers to make the fit
                      if(zSlot>=48) continue; // Do not include MicroMegas in the fit
                      if(nhit_fit==0){
                        Zlayer[nhit_fit] = zpos; // Convert coordinates to mm
                        Xlayer[nhit_fit] = xpos;
                        Ylayer[nhit_fit] = ypos;
                        nhit_fit++;
                      }
                      else{
                        if (nhit_fit < maxHits){
                          if(NbHitLayer[zSlot] >1 ) continue ; // keep only 1st hit in layers
                          // Select hits that have almost the same slope as the previous one            
                          Double_t slopy = TMath::ATan2( (ypos-Ylayer[nhit_fit-1]) , (zpos-Zlayer[nhit_fit-1]) );
                          Double_t slopx = TMath::ATan2( (xpos-Xlayer[nhit_fit-1]) , (zpos-Zlayer[nhit_fit-1]) );
                          if ( TMath::Abs(slopx)> 5. || TMath::Abs(slopy)>5. ) continue;
                          Zlayer[nhit_fit] = zpos; 
                          Xlayer[nhit_fit] = xpos;
                          Ylayer[nhit_fit] = ypos;
                          nhit_fit++;
                        }
                      }
                    } // end hit

                  // Ca c'est une rustine provisoire pour eviter les doublons d'evenements
                  hNhitDeb->Fill(nhit_roi,nhit_roi_old); hNhitDebrap->Fill((float)nhit_roi/nhit_roi_old);
                  hNhitRapHits->Fill((float)nhit_roi,(float)nhit_roi/nhit_roi_old);
                  //      if(fabs( (float)nhit_roi/nhit_roi_old - 1.)<=0.12 ) continue;

                  YOrg=hxyTot->ProjectionX()->GetMaximumBin()-1;
                  XOrg=hxyTot->ProjectionY()->GetMaximumBin()-1;
                  hxyOrg->Fill(YOrg,XOrg);
                  float RmsY = hxyTot->GetRMS(1);
                  float RmsX = hxyTot->GetRMS(2);
                  hRMSx->Fill(nhit_roi,RmsX);
                  hRMSy->Fill(nhit_roi,RmsY);

                  hnhit_roi->Fill(nhit_roi);
                
                  // Compute Hit Density to find MIP events
                  int NbLayers=0;
                  for (int layer=0 ; layer < 50 ; layer++){
                    if(NbHitLayer[layer]>0)
                      {
                        NbLayers++;
                      }
                  }

                  float HitDensity = 1. * nhit_roi/NbLayers;
                  hDensity->Fill(HitDensity);

                  //Look for penetrating MIPS-----------------------------

                  int NbLayersForw=0;
                  int NbLayersBack=0;

                  IsPenetratingMIP = 0;
                  if(HitDensity>1 && HitDensity < 3){

                    for (int layer=0 ; layer < 10 ; layer++){
                      if(0 < NbHitLayer[layer] && NbHitLayer[layer] <= 3){NbLayersForw++;}
                    }
                    for (int layer=38 ; layer < 48 ; layer++){
                      if(0 < NbHitLayer[layer]  && NbHitLayer[layer] <= 3){NbLayersBack++;}
                    }

                    hnHitForw->Fill(NbLayersForw);
                    hnHitForw->Fill(NbLayersBack+20);

                    //Define penetrating MIP
                    if(NbLayersForw>6 && NbLayersBack > 6) { IsPenetratingMIP = 1;}
                  } // end of definition of Penetrating MIPS

                  hS3S9->Fill( (float)nhit_roi, (float)NbHitLayer33/NbHitLayer77);
                  // Looking for Noise !------------------------------
                  if(ParticleType=="Noise"){
                      for (int nly=0 ; nly<=47 ; nly++){
                        nHitsTot = nHitsTot +  hxyMap[nly]->GetEntries() ;
                        hNoiseLay->Fill(nly, hxyMap[nly]->GetEntries());
                        hNbHitsPions->Fill( (int) nHitsTot );                 }
                  }
                  // Looking for Pions ---------------------------------------------------
                  //                      if( nhit_roi>150 && HitDensity>3  && ParticleType == "Pions"){
                  if( nhit_roi>80 && HitDensity>4 && ParticleType == "Pions" ){
                    hxy->Fill(YOrg,XOrg);
                    
                    // if(nhwrt<15) {tf.cd();nhwrt++; hxyTot->Write();f.cd();} 
                    int lG1=4 ; int lG2=15 ; int lG30=30;
                    int NbHitsG1 = 0;       int NbHitsG2 = 0; int NbHitsG3 = 0;  double NbHits30=0; NbHitsRear=0; NbHitsBoard=0;
                    int zstartFound=0; zstart=0;
                    for (int kl=0 ; kl<=47 ; kl++){
                      if(kl <=lG1 ) { NbHitsG1 = NbHitsG1 + (int)hxyMap[kl]->GetEntries();}
                      if(kl > lG1 && kl< lG2 ) { NbHitsG2 = NbHitsG2+ (int)hxyMap[kl]->GetEntries();}
                      if(kl > lG2 ) { NbHitsG3 = NbHitsG3 + (int)hxyMap[kl]->GetEntries();}
                      if(kl <= lG30 ) { NbHits30 = NbHits30 + hxyMap[kl]->GetEntries();}
                      if(kl >= 41 && kl< 47 ) { NbHitsRear = NbHitsRear + hxyMap[kl]->GetEntries();}

                      /******** rising number of hits in 3 consecutive layer *********/
                                
                      if ( zstartFound==0 && ((int)hxyMap[kl+2]->GetEntries() > (int)hxyMap[kl+1]->GetEntries()) && 
                           ((int)hxyMap[kl+1]->GetEntries() > (int)hxyMap[kl]->GetEntries()) ){
                        for (int k=0 ; k<3 ; k++){
                          /******* shower begins when nhit >= 4 ***********/
                          if ((int)hxyMap[kl+k]->GetEntries() >=4){
                            zstartFound=1;
                            zstart = kl+k;
                            break;
                          }
                        }
                      }
                    }
                    RG1G3 = (float)NbHitsG1 / NbHitsG3;
                    hRG1G2->Fill( (float)NbHitsG1/NbHitsG3 );
                    hnHitsRG1G2->Fill(nhit_roi, (float)NbHitsG1/NbHitsG3 );
                    ShCont1 = NbHitsRear / NbHits30;


                    //--------------Define Pions---------------
                    if(RG1G3 < 0.5 ) {
                      if(   // Fiducial cuts   
                         (XPads[jBeam]-XPadsL[jBeam] < XOrg && XOrg < XPads[jBeam]+XPadsH[jBeam]) && 
                         (YPads[jBeam]-YPadsL[jBeam] < YOrg && YOrg < YPads[jBeam]+YPadsH[jBeam]) ) {
                        //      cout << " Hidsto" <<ShCont1 << " " << NbHitsRear << endl;
        hShCont1->Fill(ShCont1); hNbHitsRear->Fill(NbHitsRear);
                        int nHits48 = (int)hxyMap[48]->GetEntries();
                        int nHits49 = (int)hxyMap[49]->GetEntries();
                        hHitsTot48->Fill(nHits48); 
                        hHitsTot49->Fill(nHits49);
                        for (int nly=0 ; nly<=47 ; nly++){
                          if( EffMap[nly] > 0) {nHitsTot  = nHitsTot +  (AverageMult)/MultMap[nly]  *  hxyMap[nly]->GetEntries() / EffMap[nly] ;}
                          if( EffMap[nly] > 0) {nHitsTot1 = nHitsTot1 +  (AverageMult)/MultMap[nly] * hxyMap1[nly]->GetEntries() / EffMap[nly] ;}
                          if( EffMap[nly] > 0) {nHitsTot2 = nHitsTot2 +  (AverageMult)/MultMap[nly] * hxyMap2[nly]->GetEntries() / EffMap[nly] ;}
                          if( EffMap[nly] > 0) {nHitsTot3 = nHitsTot3 +  (AverageMult)/MultMap[nly] * hxyMap3[nly]->GetEntries() / EffMap[nly] ;}
                          NbHitsBoard = NbHitsBoard + (hxyMap[nly]->GetEntries() - hxyMap[nly]->Integral(9,88,9,88));
                        }
                        //  if(nHitsTot<180 && nhwrt<50) {tf.cd();nhwrt++; hxyTot->Write();f.cd();} 
                        hNbHitsBoard->Fill(NbHitsBoard);
                        hNbHitsPions->Fill( (int) nHitsTot );
                        hNbHitsPions1->Fill( (int) nHitsTot1 );
                        hNbHitsPions2->Fill( (int) nHitsTot2 );
                        hNbHitsPions3->Fill( (int) nHitsTot3 );
                        if(ShCont1 <0.1  || NbHitsRear < 6) {hNbHitsPionsSC->Fill( (int) nHitsTot );}
                        //      cout << "Tree " << ShCont1 << " " << NbHitsRear << endl;
                        Ptree->Fill(); 
                      }//end of fiducial cuts
                    } // end of selection on g1/g3
                  } // End of Pions

                  // Now I select penetrating Muons
                  //      cout << IsPenetratingMIP << "  " << nhit_roi << endl;
                  if(IsPenetratingMIP==1 && nhit_roi<150 && nhit_roi>50 && ParticleType == "Muons" ){
                    float zmax = 1372.;
                    TGraphErrors *grx = new TGraphErrors(nhit_fit,Zlayer,Xlayer,0,ErPos);
                    TF1 *f1 = new TF1("f1","pol1",0,zmax);
                    grx->Fit("f1","Q");
                    double betaX = f1->GetParameter(0);
                    double alfaX = f1->GetParameter(1);
                    double chi2X = f1->GetChisquare();
                    delete f1; delete grx;
                    TGraphErrors *gry = new TGraphErrors(nhit_fit,Zlayer,Ylayer,0,ErPos);
                    TF1 *f2 = new TF1("f2","pol1",0,zmax);
                    gry->Fit("f2","Q");
                    double betaY = f2->GetParameter(0);
                    double alfaY = f2->GetParameter(1);
                    double chi2Y = f2->GetChisquare();
                    delete f2; delete gry;
                   //  if(nhwrt<32 & nhit_roi>8 ) { // Now check for doublons
//                    tf.cd();
//                    nhwrt++; 
//                    cout << " num Evt local " << nbEvt << " new/old  Nb of hits "  <<  nhit_roi << " / " << nhit_roi_old  << " new/old delatMax "  << maxDeltatOld << " "  <<  maxDeltat << " BcId " << BcIdAbs << " / " << BcIdAbsOld << endl ; 
//                    hxyTot->Write();
//                    f.cd();
//                  } 
                    for(int i=48; i<=49 && ( chi2X/(nhit_fit-2) < 6. && chi2Y/(nhit_fit-2) < 6 ); i++){
                      float Xextrapol = alfaX*i*28. + betaX ; // X extrapolate to chamber i in mm
                      float Yextrapol = alfaY*i*28. + betaY ; // Y extrapolate to chamber i in mm
                      hXYExtr->Fill(Yextrapol/10.,Xextrapol/10.);
                      // Find Central pad from computed x,y position
                      Int_t xpadCentr = -1000;
                      Int_t ypadCentr = -1000;
                      if( i>47 && ( fabs(Yextrapol-322.5) < 20 || fabs(Yextrapol-650.) < 20 ) ) continue;
                      if( i>47 && ( fabs(Xextrapol-485.0) < 20  ) ) continue;
                      if(i<48){
                        xpadCentr = (int)(Xextrapol/10.4);
                        ypadCentr = (int)(Yextrapol/10.4);
                      }
                      else{
                        if(Xextrapol<=480.) {xpadCentr = (int)(Xextrapol/10.);}
                        else {xpadCentr = (int)( (Xextrapol-5)/10 );}
                        if(Yextrapol<=320.) {ypadCentr = (int)( (Yextrapol/10.) );}
                        if(Yextrapol>325 && Yextrapol<=645.) {ypadCentr = (int)( (Yextrapol-5.)/10.);}
                        if(Yextrapol>=646) {ypadCentr = (int)((Yextrapol-10.)/10.);}
                      }
                      hPadCentr->Fill(ypadCentr,xpadCentr);                 
                      if (xpadCentr >8 && ypadCentr >8 && xpadCentr <88 && ypadCentr < 88 ){
                        //                    if (xpadCentr >0 && ypadCentr >0 ){
                        int nHits48 = (int)hxyMap[48]->GetEntries();
                        int nHits49 = (int)hxyMap[49]->GetEntries();
                        if(i==48){hHitsTot48->Fill(nHits48); } // cout << " Nb hits 48 " << nHits48 << endl;
                        if(i==49){hHitsTot49->Fill(nHits49);}

                        

                      }  // end checking the pad center
                    } // end loop on all layers 
                  } // end penetrating muon
         
                  if (nbEvtTot%10000 == 0 ) {
                    cout << "Processing event: " << nbEvtTot
                         << " Local event " << nbEvt << endl;  
                  }  
                  tf.cd(0);
                }  // end event 
              cout << " Finished with events " << nbEvt << " in this TREE " << endl;
            } // else 
        } // end of tree or run
      cout << " Ending  file in chain " << endl;
    } // end of chain
  cout << " Going to finsh the work after" << nbEvtTot << " events "<< endl;
  tf.cd();
  hnhit_roi->Write();
  hDensity->Write();
  hHitsTot48->Write(); 
  hHitsTot49->Write();
  Double_t b148=hHitsTot48->GetBinContent(1); Double_t ent48=hHitsTot48->GetEntries(); cout << " Innefficiency 48 " << b148/ent48 << " 1st bin " << b148 << " All " << ent48 << endl;
  Double_t b149=hHitsTot49->GetBinContent(1); Double_t ent49=hHitsTot49->GetEntries(); cout << " Innefficiency 49 " << b149/ent49 << endl;
  hnHitForw->Write();
  hxy->Write();
  hRegMult->Write();
  hRegMultC->Write();
  hPadCentr->Write();
  hXYExtr->Write();
  hS3S9->Write();
  hRMSx->Write();
  hRMSy->Write();
  hNhitDeb->Write();
  hNhitDebrap->Write();hNhitRapHits->Write();
  hRG1G2->Write() ; hnHitsRG1G2->Write();
  hNbHitsPions->Write(); hNbHitsPions1->Write(); hNbHitsPions2->Write(); hNbHitsPions3->Write(); 
  hNbHitsPionsSC->Write(); hNoiseLay->Write();
  hxyOrg->Write();
  hShCont1->Write();hNbHitsRear->Write();hNbHitsBoard->Write();

  PionsRootFile->cd();
  Ptree->Write();

  }

---