/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/Iro/profile4layers_simu.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) $ */
#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 <TRandom.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 <TF1.h>
#include "tools/Arguments.hh"



using namespace std;


const int maxHits = 10000;
const int maxEvents = 10000;


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

  Arguments arg(argc,argv," inputrootfilename outputrootfilename hitcutstart");

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



  int Nsample = 3;
  int z_sample[3] = {14,24,45};
  float calib[3]  = {1,1,1};

  TString infilename = argv[1];
  TString outfilename = argv[2];

  int hitcutstart = atoi(argv[3]);



  float calib_layer[50];
  float layer_response = 1.0;
  float layer_response_rel_rms = 0.25;



  int Nlayer = 50;

  int nhit_calo_cut = 60;

  int nshower = 0;
  int nmuon = 0;
  int npenetrating_muon = 0;

  int nstart = 0;

  int nhit_cumul = 0;

  int zstart = -1;

  int nevt_tot = 0;






  TH1F * hdensity = new TH1F("hdensity","",100,0,50);
  TH2F * hnhit_density = new TH2F("hnhit_density","",1000,0,2000,200,0,100);

  TH1F * hrhadron = new TH1F("hrhadron","",1000,0,10);
  TH2F * hnhit_rhadron = new TH2F("hnhit_rhadron","",1000,0,2000,200,0,10);

  TProfile * hlong1 = new TProfile("hlong1","",50,0,50,0,150,"");
  TProfile * hlong2 = new TProfile("hlong2","",50,0,50,0,150,"");
  TProfile * hlong3 = new TProfile("hlong3","",50,0,50,0,150,"");

  TGraphErrors * tgmip = new TGraphErrors();
  TGraphErrors * tgshower = new TGraphErrors();

  tgmip->SetName("tgmip");
  tgshower->SetName("tgshower");

  TH1F * hlayer_mult = new TH1F("hlayer_mult","",50,0,50);

  TH1I * hnhit_calo = new TH1I("hnhit_calo","",500,0,5000);
  TH1I * hnhit_evt  = new TH1I("hnhit_evt","all showers",500,0,5000);
  TH1I * hnhit_evt2 = new TH1I("hnhit_evt2","all showers starting wthin 5 first planes",500,0,5000);
  TH1I * hnhit_evt3 = new TH1I("hnhit_evt3","all showers without MIP part",500,0,5000);
  TH1I * hnhit_evt4 = new TH1I("hnhit_evt4","all showers without MIP part starting wthin 5 first planes",500,0,5000);
  TH1I * hnhit_evt5  = new TH1I("hnhit_evt5","all showers with multiplicity per layer",500,0,5000);
  TH1I * hnhit_evt6 = new TH1I("hnhit_evt6","all showers without MIP part starting wthin 5 first planes with multiplicity per layer",500,0,5000);

  TH1I * hstart = new TH1I("hstart","",Nlayer,0,Nlayer);

  TH1I * hnhit_layer = new TH1I("hnhit_layer","",200,0,200);

  TH1I * hz_nhit = new TH1I("hz_nhit","",50,0,50);
  TH1I * hz_nhit_cumul = new TH1I("hz_nhit_cumul","",50,0,50);

  TH2I * hnhit_zstart = new TH2I("hnhit_zstart","Nhit for pions vs shower start",50,0,50,400,0,2000);
  TH2I * hnhit2_zstart = new TH2I("hnhit2_zstart","Nhit after removal of MIP hits for pions vs shower start",50,0,50,400,0,2000);


  TString name,title;
  TH2F * hxy[Nlayer];
  TH1F * hx[Nlayer];
  TH1F * hy[Nlayer];
  TH1I * hnhit_mip[Nlayer];
  TH1I * hnhit_shower[Nlayer];
  TH1I * hnhit_shower2[Nlayer];
  TH1I * hnhit_noise[Nlayer];
  TH1I * hdt[Nlayer];
  int nhit[Nlayer];
  int nnoise[Nlayer];

  for (int i=0;i<Nlayer;i++)
    {
      nhit[i] = 0;
      nnoise[i] = 0;

      name = Form("hnhit_mip_%i",i);
      title = Form("Number of CaloHit - chb %i",i);
      hnhit_mip[i]= new TH1I(name,title,100,0,100);

      name = Form("hnhit_shower_%i",i);
      title = Form("Number of CaloHit - chb %i",i);
      hnhit_shower[i]= new TH1I(name,title,1000,0,1000);

      name = Form("hnhit_shower2_%i",i);
      title = Form("Number of CaloHit - chb %i",i);
      hnhit_shower2[i]= new TH1I(name,title,1000,0,1000);

      name = Form("hnhit_noise_%i",i);
      title = Form("Number of NoiseHit - chb %i",i);
      hnhit_noise[i]= new TH1I(name,title,100,0,100);

      name = Form("hdt_%i",i);
      title = Form("Delta_t - chb %i",i);
      hdt[i]= new TH1I(name,title,21,-10,10);

      name = Form("hxy_%i",i);
      title = Form("CaloHit position - chb %i",i);

      if (i<48)
        {
          hxy[i] = new TH2F(name,title,96,0,96*1.04,96,0,96*1.04);

          name = Form("hx_%i",i);
          title = Form("X - chb %i",i);
          hx[i] = new TH1F(name,title,96,0,96*1.04);

          name = Form("hy_%i",i);
          title = Form("Y - chb %i",i);
          hy[i] = new TH1F(name,title,96,0,96*1.04);
        }
      else
        {
          hxy[i] = new TH2F(name,title,194,0,97,193,0,96.5);

          name = Form("hx_%i",i);
          title = Form("X - chb %i",i);
          hx[i] = new TH1F(name,title,194,0,97);

          name = Form("hy_%i",i);
          title = Form("Y - chb %i",i);
          hy[i] = new TH1F(name,title,193,0,96.5);
        }
    }




  for (int i=0 ; i<50 ; i++)
    {
      calib_layer[i] = layer_response * fabs(gRandom->Gaus(1,layer_response * layer_response_rel_rms));
      hlayer_mult->Fill(i,calib_layer[i]);
    }
  






  TFile f(infilename);


  TIter nextkey(f.GetListOfKeys());
  TKey *key;

  int NbTrees = 0;

  while (key = (TKey*)nextkey()) 
    {
      TString keyname = key->GetName();
      TTree *tree = (TTree*)key->ReadObj();                
      cout << "---------------------------- NEW TTREE -----------------"<< endl;

      int nEvts = tree->GetEntries();

      int nHits;
      int HitX[maxHits];
      int HitY[maxHits];
      int HitZ[maxHits];
      float HitEng[maxHits];
      tree->SetBranchAddress("nHits", &nHits);
      tree->SetBranchAddress("HitX", &HitX);
      tree->SetBranchAddress("HitY", &HitY);
      tree->SetBranchAddress("HitZ", &HitZ);
      tree->SetBranchAddress("HitEng", &HitEng);


      // Loop over events                                                                                                                                                                                                

      for (int evt = 0; evt < nEvts; evt++)
        {
          if (evt >= maxEvents) continue;

          tree->GetEntry(evt);

          int nhit_seperate1 = 0;
          int nhit_seperate2 = 0;
          float rhadron = 0;

          int nhit_evt = 0;
          int nhit_evt2 = 0;
          int nhit_evt3 = 0;
          int nhit_evt4 = 0;
          float nhit_evt5 = 0;
          float nhit_evt6 = 0;

          int nhotlayer = 0;
          float density = 0;

          int nhot_rear  = 0;
          int nhot_front = 0;

          for(int i=0 ; i < Nlayer ; i++)
            {
              nhit[i] = 0;
            }
                
          zstart = -1;


          cout<<"evt "<<evt<<" -> nhit = "<<nHits<<" \r"<<flush;

          hnhit_calo->Fill(nHits);

          


          /********** read of hits ************/
                                                                                                                                                                            
          for (int hit = 0; hit < nHits; hit++)
            {
              int xpos = HitX[hit] + 48;
              int ypos = HitY[hit] + 48;
              int zSlot = HitZ[hit];
              float Edep = HitEng[hit];

              if (Edep !=0 )
                {
                  hxy[zSlot]->Fill(xpos,ypos);
                        
                  nhit[zSlot]++;
                }
            }   




          /********** sum of hits ************/

          for(int i=0 ; i < Nlayer ; i++)
            {
              nhit_evt += nhit[i];
              nhit_evt5 += (nhit[i] * calib_layer[i]);
              
              if (nhit[i] != 0)
                {
                  density += nhit[i];
                  nhotlayer++;

                  if (i<5){nhit_seperate1++;}
                  if (i>=15){nhit_seperate2++;}

                  if (i<10){nhot_front++;}
                  if (i>=19 && i<30){nhot_rear++;}
                }
            }






          /************ calculate density ****************/

          if (nhotlayer != 0)
            {
              density /= (1.*nhotlayer);
              hdensity->Fill(density);
            }

          if (nhit_seperate2 != 0)
            {
              rhadron = 1. * nhit_seperate1 / nhit_seperate2;
            }




          /******** select muons *******/

          if (density !=0 && density < 2)
            {
              /********** select penetrating ************/

              nmuon++;

              if (nhot_front >= 6 && nhot_rear >= 6)
                {
                  
                  npenetrating_muon++;
                  
                  for(int i=0 ; i < Nlayer ; i++)
                    {
                      hnhit_mip[i]->Fill(nhit[i]);
                    }
                }
            }




          /******** select showers *******/

          if (density >= 5 && nhit_evt >= nhit_calo_cut)
            {          
              
              hrhadron->Fill(rhadron);
              hnhit_rhadron->Fill(nhit_evt,rhadron);

              nshower++;
              
              /******** find shower start *********/
              
              for(int i=0 ; i < Nlayer-2 ; i++)
                {
                  
                  /******** rising number of hits in 3 consecutive layer *********/
                  
                  if ( (nhit[i+2] > nhit[i+1]) && (nhit[i+1] > nhit[i]) )
                    {
                      for (int k=0 ; k<3 ; k++)
                        {                         

                          /******* shower begins when nhit >= 4 ***********/

                          if (nhit[i+k] >= hitcutstart)
                            {
                              nstart++;
                              zstart = i+k;
                              
                              hstart->Fill(zstart);
                              k = 3;
                              i = Nlayer;
                            }
                        }
                    }
                }
              
              
              
              /********* fill profiles *********/
              
              if (zstart != -1)
                {

                  hnhit_zstart->Fill(zstart,nHits);
                  
                  for(int i=0 ; i < Nlayer ; i++)
                    {
                      hnhit_shower[i]->Fill(nhit[i]);
                      
                      hz_nhit->Fill(i,nhit[i]);
                      hlong1->Fill(i,nhit[i]);

                      //all showers starting wthin 5 first planes
                      if (zstart < 5)
                        {
                          nhit_evt2 += nhit[i];
                        }
                      
                      //all showers without MIP part
                      if (i >= zstart && i <= 45)
                        {
                          nhit_evt3 += nhit[i];

                          //all showers without MIP part starting wthin 5 first planes
                          if (zstart < 5)
                            {
                              nhit_evt4 += nhit[i];
                              nhit_evt6 += (nhit[i] * calib_layer[i]);
                            }
                        }
                    }                                               


                  //all showers
                  hnhit_evt->Fill(nhit_evt);
                  hnhit_evt5->Fill(nhit_evt5);

                  //all showers starting wthin 5 first planes
                  if (nhit_evt2 != 0){hnhit_evt2->Fill(nhit_evt2);}

                  //all showers without MIP part
                  hnhit_evt3->Fill(nhit_evt3);

                  //all showers without MIP part starting wthin 5 first planes
                  if (nhit_evt4 != 0){hnhit_evt4->Fill(nhit_evt4);hnhit_evt6->Fill(nhit_evt6);}

                  


                  for (int k=0 ; k<Nsample ; k++)
                    {
                      if (zstart <= z_sample[k])
                        {
                          hnhit_shower2[z_sample[k] - zstart]->Fill(nhit[z_sample[k]]);
                          hlong2->Fill(z_sample[k] - zstart,nhit[z_sample[k]]);
                          hlong3->Fill(z_sample[k] - zstart,nhit[z_sample[k]] * calib[k]);
                        }
                    }
                }
            }
        }
    }



  cout<<endl;
  cout<<endl;
  cout<<"Nevt_tot = "<<nevt_tot<<endl;
  cout<<endl;
  cout<<"Nshower = "<<nshower<<endl;
  cout<<"Nstart = "<<nstart<<endl;
  cout<<endl;
  cout<<"Nmuon = "<<nmuon<<endl;
  cout<<"Npenetrating_muon = "<<npenetrating_muon<<endl;
  cout<<endl;





  cout<<endl;
  cout<<"Write to TFile "<<outfilename<<endl;
  cout<<endl;
  
  TFile * tf = new TFile(outfilename,"RECREATE");
  
  for (int i=0;i<Nlayer;i++)
    {
      tgmip->SetPoint(i,i,hnhit_mip[i]->GetMean());
      tgmip->SetPointError(i,0,hnhit_mip[i]->GetRMS());
      
      tgshower->SetPoint(i,i,hnhit_shower[i]->GetMean());
      tgshower->SetPointError(i,0,hnhit_shower[i]->GetRMS());
      
      hxy[i]->Write();
      hdt[i]->Write();
      hnhit_mip[i]->Write();
      hnhit_shower[i]->Write();
      hnhit_shower2[i]->Write();
    }
  
  tgmip->SetMarkerStyle(20);
  tgmip->SetMarkerColor(2);
  
  tgshower->SetMarkerStyle(20);
  tgshower->SetMarkerColor(3);
  
  tgmip->Write();
  tgshower->Write();

  hlayer_mult->Write();

  hrhadron->Write();
  hnhit_rhadron->Write();
  hnhit_zstart->Write();
  hdensity->Write();
  hnhit_density->Write();
  hstart->Write();
  hlong1->Write();
  hlong2->Write();
  hlong3->Write();
  hnhit_evt->Write();
  hnhit_evt2->Write();
  hnhit_evt3->Write();
  hnhit_evt4->Write();
  hnhit_evt5->Write();
  hnhit_evt6->Write();
  hnhit_layer->Write();
  hz_nhit->Write();
  hz_nhit_cumul->Write();
  hnhit_calo->Write();
  
  tf->Close();

  
}

---