/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/Iro/myfirstprogam.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 <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;


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

  Arguments arg(argc,argv," run1 [run2] max [run20] -o outputrootfilename");

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



  int Nsample = 3;

  int z_sample[3] = {14,24,39};

  float calib[3] = {1.8/1.641,1.8/1.587,1.8/1.611};

  int nhitcutstart = 9;

  int nhit_calo_cut = 220;  


  int NbOfRuns = arg.getNbOfArgs();

  TString runnumber[20];
  TString outfilename = "temp.root";

  for (int nr=0 ; nr < NbOfRuns ; nr ++ )
    {
      runnumber[nr] = arg.getArgument(nr+1); 
    }

  if ( arg.getOption("-o").size() != 0 )
    {
      outfilename  = arg.getOption("-o");
    }

  if ( arg.getOption("-c").size() != 0 )
    {
      nhitcutstart  = atoi(arg.getOption("-c").c_str());

      TString temp = arg.getOption("-c");
      for (int i=0;i<5;i++){outfilename.Chop();}
      outfilename += "_cut" + temp + ".root";
    }


  if ( arg.getOption("-n").size() != 0 )
    {
      nhit_calo_cut = atoi(arg.getOption("-n").c_str());
      cout<<"Nhit_calo_cut = "<<nhit_calo_cut<<endl;
    }

   


  int Nlayer = 50;

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

  int nstart = 0;

  int nhit_evt = 0;
  int nhit_cumul = 0;

  int zstart = -1;

  int nevt_tot = 0;

  int ndouble = 0;

  int nbHits_prev = -1;





  TH1F * hdensity = new TH1F("hdensity","",200,0,100);
  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);

  TH2I * hbug = new TH2I("hbug","",100,0,1000,100,0,1000);

  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();
  TGraphErrors * tgnoise = new TGraphErrors();

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

  TH1I * hnhit_calo = new TH1I("hnhit_calo","",1000,0,10000);
  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 * 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);

  TH1F * hz_rmsx = new TH1F("hz_rmsx","",50,0,50);
  TH1F * hz_rmsy = new TH1F("hz_rmsy","",50,0,50);

  TH1F * hx_evt = new TH1F("hx_evt","",96,0,96*1.04);
  TH1F * hy_evt = new TH1F("hy_evt","",96,0,96*1.04);

  TH2I * hnhit_zstart = new TH2I("hnhit_zstart","",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);
        }
    }







  TChain chain("T");

  for ( int nr=0; nr< NbOfRuns; nr++){

    for (int i=0 ; i<4 ; i++)
      {
        TString SubFile = Form("_I%i",i);
        TString fname = Form("/lapp_data/LC/data/TB2012/SPS_H2_may_2012/Production/CaloHits/CaloHit_DHCAL_" + runnumber[nr] + SubFile + "_0.slcio.root");
        chain.Add(fname);
        cout << fname << endl;
      }
  }

  TObjArray *fileElements=chain.GetListOfFiles();
  TIter next(fileElements);
  TChainElement *chEl=0;
  int NchainFiles=0;

  while (( chEl=(TChainElement*)next() ))
    {
      TFile f(chEl->GetTitle());
      TIter nextkey(f.GetListOfKeys());
      TKey *key;

      int NbTrees = 0;

      while (key = (TKey*)nextkey()) 
        {

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

            TTree *tree = (TTree*)key->ReadObj();                
            cout << "---------------------------- NEW TTREE -----------------"<< endl;
            int nbEvt =0;
            int nbHit = 0;
            NbTrees++;
            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 ;
                  }
              }

            catch ( ... ) {}

            int RunId=1;
            CaloEvent *evt =  new CaloEvent();
            TBranch *branch= tree->GetBranch("CaloEvent");
            branch->SetAddress(&evt);

            CaloHit* hit =NULL;
            int nbEntries = tree->GetEntries();




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

                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;
                
                int nbHits = evt->GetNHits();
                int maxDeltat = evt->GetDeltaTmax();

                int nhotlayer = 0;
                float density = 0;

                int nhot_rear  = 0;
                int nhot_front = 0;

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


                if (evtNum!=0){hbug->Fill(nbHits,nbHits_prev);}

                nevt_tot++;

                    
                    //cout<<"Evt "<<evtNum<<" Nhit = "<<nbHits<<" - Nhit_prev = "<<nbHits_prev<<endl;
                    cout<<"evt "<<evtNum<<" -> nhit = "<<nbHits<<" \r"<<flush;

                    hnhit_calo->Fill(nbHits);



                    /******** time cut ******/

                    for(int i=0 ; i < nbHits ; i++)
                      {
                        hit = (CaloHit*)evt->GetHits()->UncheckedAt(i);
                        
                        int deltat = hit->GetDeltaT();
                        int zSlot = hit->GetLayer();
                        
                        float xpos =  hit->GetX()*1e-4;
                        float ypos =  hit->GetY()*1e-4;
                        
                        hxy[zSlot]->Fill(xpos,ypos);
                        
                        hx[zSlot]->Fill(xpos);
                        hy[zSlot]->Fill(ypos);
                        
                        hx_evt->Fill(xpos);
                        hy_evt->Fill(ypos);

                        hdt[zSlot]->Fill(maxDeltat-deltat);
                        
                        if (zSlot >= 48)
                          {
                            if  (fabs(maxDeltat-deltat) <= 2)
                              {
                                nhit[zSlot]++;
                              }
                            else
                              {
                                nnoise[zSlot]++;
                              }
                          }
                        else
                          {
                            if  (fabs(maxDeltat-deltat) <= 1)
                              {
                                nhit[zSlot]++;
                              }
                            else
                              {
                                nnoise[zSlot]++;
                              }
                          }
                      }

                


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


                    for(int i=0 ; i < Nlayer ; i++)
                      {
                        hnhit_noise[i]->Fill(nnoise[i]);

                        nhit_evt += nhit[i];
                        
                        if (nhit[i] != 0)
                          {
                            hnhit_layer->Fill(nhit[i]);

                            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);
                        hnhit_density->Fill(nhit_evt,density);
                      }

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




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

                    if (density !=0 && density < 5)
                      {
                        /********** 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++;

                        /********** select pions - reject electrons ************/

                        if (rhadron <= 0.4)
                          {                     
                            nhshower++;
                        
                            /******** 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] >= nhitcutstart)
                                          {
                                            nstart++;
                                            zstart = i+k;
                                            
                                            hstart->Fill(zstart);
                                            k = 3;
                                            i = Nlayer;
                                          }
                                      }
                                  }
                              }
                            


                            /********* fill profiles *********/
                            
                            if (zstart != -1)
                              {
                                
                                hnhit_zstart->Fill(zstart,nhit_evt);
                                
                                
                                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];
                                          }
                                      }
                                  }
                                
                                
                                //all showers
                                hnhit_evt->Fill(nhit_evt);
                                
                                //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);}

                  

                                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]);
                                      }
                                  }
                              }
                          }
                      }

                    if (nbHits_prev == nbHits){ndouble++;}

                nbHits_prev = nbHits;
              }          
        }
    }



  cout<<endl;
  cout<<endl;
  cout<<"Nevt_tot = "<<nevt_tot<<endl;
  cout<<"Ndouble = "<<ndouble<<endl;
  cout<<endl;
  cout<<"Nshower = "<<nshower<<endl;
  cout<<"Nhadronshower = "<<nhshower<<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());
      
      tgnoise->SetPoint(i,i,hnhit_noise[i]->GetMean());
      tgnoise->SetPointError(i,0,hnhit_noise[i]->GetRMS());
      
      hxy[i]->Write();
      hdt[i]->Write();
      hnhit_mip[i]->Write();
      hnhit_shower[i]->Write();
      hnhit_shower2[i]->Write();
      hnhit_noise[i]->Write();
    }
  
  tgmip->SetMarkerStyle(20);
  tgmip->SetMarkerColor(2);
  
  tgshower->SetMarkerStyle(20);
  tgshower->SetMarkerColor(3);
  
  tgnoise->SetMarkerStyle(20);
  tgnoise->SetMarkerColor(4);
  
  tgmip->Write();
  tgshower->Write();
  tgnoise->Write();

  hrhadron->Write();
  hnhit_rhadron->Write();
  hnhit_zstart->Write();
  hbug->Write();
  hdensity->Write();
  hnhit_density->Write();
  hstart->Write();
  hlong1->Write();
  hlong2->Write();
  hlong3->Write();
  hx_evt->Write();
  hy_evt->Write();
  hnhit_evt->Write();
  hnhit_evt2->Write();
  hnhit_evt3->Write();
  hnhit_evt4->Write();
  hnhit_layer->Write();
  hz_nhit->Write();
  hz_nhit_cumul->Write();
  hz_rmsx->Write();
  hz_rmsy->Write();
  hnhit_calo->Write();
  
  tf->Close();

  
}

---