rsAnalc0p1.cpp

Go to the documentation of this file.

#include <XANADOO/domain/interface/XANAEsdEvent.h>
#include <XANADOO/XANAClusters/interface/XANACluster.h>
#include <XANADOO/XANAClusters/interface/XANACaloRecHit.h>
#include <XANADOO/XANAClusters/interface/XANAEmCluster.h>
#include <XANADOO/XANAClusters/interface/XANASuperCluster.h>
#include <XANADOO/XANAClusters/interface/XANAHadCluster.h>
#include <XANADOO/XANATracks/interface/XANATrack.h>
#include <XANADOO/XANAMcInfo/interface/XANAGeneratorEvent.h>
#include <XANADOO/XANAMcInfo/interface/XANAGeantEvent.h>
#include <XANADOO/XANAMcInfo/interface/XANAGeneratorParticle.h>
#include <XANADOO/XANAElectronCandidate/interface/XANAElectronCandidate.h>
#include <XANADOO/XANAElectronCandidate/interface/XANAElectronTrack.h>
#include <XANADOO/XANAElectronCandidate/interface/XANAElectronSeed.h>
#include <XANADOO/XANAMuons/interface/XANAMuonCandidate.h>
#include <XANADOO/XANATriggerInfo/interface/XANATriggerInfo.h>
#include <XANADOO/XANAJets/interface/XANAJet.h>

#include <XANADOO/XANAGeometry/interface/XANAEcalGeometry.h>

#include <XANADOO/XANAnalysisRS/interface/RSElecSelector.h>
#include <XANADOO/XANAnalysisRS/interface/RSCorrectionFunc.h>
#include <XANADOO/XANAnalysisRS/interface/RSSaturation.h>
#include <XANADOO/XANAnalysisRS/interface/RSPhotonAssoc.h>
#include <XANADOO/XANAnalysisRS/interface/ClassRSElec.h>
#include <XANADOO/XANAnalysisRS/interface/RSVertexCorr.h>
#include <XANADOO/XANAnalysisRS/interface/RSDistribAng.h>
#include <XANADOO/XANAnalysisRS/interface/RSMassEstim.h>

#include "TTree.h"
#include "TBits.h"
#include "TObjectTable.h"
#include "TFile.h"
#include "TBranch.h"
#include "TClonesArray.h"
#include "TObjArray.h"
#include "TH1.h"
#include "TH2.h"
#include "TSystem.h"
#include "TChain.h"

#include <string>

using namespace std;

Short_t ne_candidatelectron=0,ns_supercluster=0;
TObjArray *genpart = 0;  
TObjArray *clusters = 0;
TObjArray *hadClusters = 0;
TObjArray *superClusters = 0;
TObjArray *preshClusters = 0;
TRefArray *hits = 0;
TObjArray *tracks = 0;
TObjArray *vertices = 0;
TObjArray *electrons = 0;
TObjArray *electronTracks = 0;
TObjArray *muons = 0;


  
int main(int argc, char** argv)
{
  //   Connect to file generated
  TFile *f;
    //  if (argc==1) f = new TFile("SE04_E25_Pall_PT30.root");
  if (argc==1) f = new TFile("/data/scratch/collard/Production_xana151/zprim/famos/xZprimALRM_m1tev.root");

  else f = new TFile(*++argv);
 
    //   Read Tree named "T" in memory. Tree pointer is assigned the same name
  TTree *T = (TTree*)f->Get("ESD");

   
// version TChain
//   TChain *T = new TChain("ESD");
//   T->Add("file1.root");
//   T->Add("file2.root");
 
// example of how to suppress reading of branches
// has to be put before setting the branch addresses!
// T->SetBranchStatus("GeantEvent*",0);

// Reconstructed event
  XANAEsdEvent *event = new XANAEsdEvent();  
  T->SetBranchAddress("Event.",&event);
  // MC information (generator and GEANT event)
  XANAGeneratorEvent *generatorEvent = new XANAGeneratorEvent();
  T->SetBranchAddress("GeneratorEvent.",&generatorEvent);
    
  XANAGeantEvent *geantEvent = new XANAGeantEvent();
  T->SetBranchAddress("GeantEvent.",&geantEvent);
  


 
  // Histo
  // sans contrainte de masse
  TH1F *histo1000 = new TH1F("histo1000","M RSgrav gen ee ",180,0.,4.5);
  TH1F *histo1001 = new TH1F("histo1001","M RSgrav gen",180,0.,4.5);
  TH1F *histo1002 = new TH1F("histo1002","M RSgrav rec",180,0.,4.5);
  TH1F *histo1003 = new TH1F("histo1003","Delta M     ",100,-50.,50.);
  TH1F *histo1004 = new TH1F("histo1004","CosTheta gen",20,-1.,1.);
  TH1F *histo1005 = new TH1F("histo1005","CosTheta rec",20,-1.,1.);
  TH1F *histo1006 = new TH1F("histo1006","CosTheta gen",10,-1.,1.);
  TH1F *histo1007 = new TH1F("histo1007","CosTheta rec",10,-1.,1.);
  TH1F *histo1008 = new TH1F("histo1008","Delta Cos th",100,-0.5,0.5);
  TH1F *histo1009 = new TH1F("histo1009","Delta Cos th",100,-0.2,0.2);
  
  // 2TeV 
  TH1F *histo2000 = new TH1F("histo2000","M RSgrav gen ee ",50, 1.5,2.5);
  TH1F *histo2001 = new TH1F("histo2001","M RSgrav gen    ",50, 1.5,2.5);
  TH1F *histo2002 = new TH1F("histo2002","M RSgrav rec    ",50, 1.5,2.5);
  TH1F *histo2003 = new TH1F("histo2003","M RSgrav gen    ",100, 1.5,2.5);
  TH1F *histo2004 = new TH1F("histo2004","M RSgrav rec    ",100, 1.5,2.5);
  TH1F *histo2005 = new TH1F("histo2005","M RSgrav gen    ",1000,1.5,2.5);
  TH1F *histo2006 = new TH1F("histo2006","M RSgrav rec    ",1000,1.5,2.5);
  TH1F *histo2007 = new TH1F("histo2007","CosTheta gen    ",20,-1.,1.);
  TH1F *histo2008 = new TH1F("histo2008","CosTheta rec    ",20,-1.,1.);
  TH1F *histo2009 = new TH1F("histo2009","CosTheta gen    ",10,-1.,1.);
  TH1F *histo2010 = new TH1F("histo2010","CosTheta rec    ",10,-1.,1.);

  // 2TeV5 
  TH1F *histo2100 = new TH1F("histo2100","M RSgrav gen ee ",50, 2.0,3.0);
  TH1F *histo2101 = new TH1F("histo2101","M RSgrav gen    ",50, 2.0,3.0);
  TH1F *histo2102 = new TH1F("histo2102","M RSgrav rec    ",50, 2.0,3.0);
  TH1F *histo2103 = new TH1F("histo2103","M RSgrav gen    ",100, 2.0,3.0);
  TH1F *histo2104 = new TH1F("histo2104","M RSgrav rec    ",100, 2.9,3.0);
  TH1F *histo2105 = new TH1F("histo2105","M RSgrav gen    ",1000,2.0,3.0);
  TH1F *histo2106 = new TH1F("histo2106","M RSgrav rec    ",1000,2.0,3.0);
  TH1F *histo2107 = new TH1F("histo2107","CosTheta gen    ",20,-1.,1.);
  TH1F *histo2108 = new TH1F("histo2108","CosTheta rec    ",20,-1.,1.);
  TH1F *histo2109 = new TH1F("histo2109","CosTheta gen    ",10,-1.,1.);
  TH1F *histo2110 = new TH1F("histo2110","CosTheta rec    ",10,-1.,1.);

  // 3TeV 
  TH1F *histo2200 = new TH1F("histo2200","M RSgrav gen ee ",50, 2.5,3.5);
  TH1F *histo2201 = new TH1F("histo2201","M RSgrav gen    ",50, 2.5,3.5);
  TH1F *histo2202 = new TH1F("histo2202","M RSgrav rec    ",50, 2.5,3.5);
  TH1F *histo2203 = new TH1F("histo2203","M RSgrav gen    ",100, 2.5,3.5);
  TH1F *histo2204 = new TH1F("histo2204","M RSgrav rec    ",100, 2.5,3.5);
  TH1F *histo2205 = new TH1F("histo2205","M RSgrav gen    ",1000,2.5,3.5);
  TH1F *histo2206 = new TH1F("histo2206","M RSgrav rec    ",1000,2.5,3.5);
  TH1F *histo2207 = new TH1F("histo2207","CosTheta gen    ",20,-1.,1.);
  TH1F *histo2208 = new TH1F("histo2208","CosTheta rec    ",20,-1.,1.);
  TH1F *histo2209 = new TH1F("histo2209","CosTheta gen    ",10,-1.,1.);
  TH1F *histo2210 = new TH1F("histo2210","CosTheta rec    ",10,-1.,1.);

  // 3TeV5
  TH1F *histo2300 = new TH1F("histo2300","M RSgrav gen ee ",50, 3.0,4.0);
  TH1F *histo2301 = new TH1F("histo2301","M RSgrav gen    ",50, 3.0,4.0);
  TH1F *histo2302 = new TH1F("histo2302","M RSgrav rec    ",50, 3.0,4.0);
  TH1F *histo2303 = new TH1F("histo2303","M RSgrav gen    ",100, 3.0,4.0);
  TH1F *histo2304 = new TH1F("histo2304","M RSgrav rec    ",100, 3.0,4.0);
  TH1F *histo2305 = new TH1F("histo2305","M RSgrav gen    ",1000,3.0,4.0);
  TH1F *histo2306 = new TH1F("histo2306","M RSgrav rec    ",1000,3.0,4.0);
  TH1F *histo2307 = new TH1F("histo2307","CosTheta gen    ",20,-1.,1.);
  TH1F *histo2308 = new TH1F("histo2308","CosTheta rec    ",20,-1.,1.);
  TH1F *histo2309 = new TH1F("histo2309","CosTheta gen    ",10,-1.,1.);
  TH1F *histo2310 = new TH1F("histo2310","CosTheta rec    ",10,-1.,1.);

  // 3TeV5 sans saturation
  TH1F *histo2400 = new TH1F("histo2400","M RSgrav gen ee ",50, 3.0,4.0);
  TH1F *histo2401 = new TH1F("histo2401","M RSgrav gen    ",50, 3.0,4.0);
  TH1F *histo2402 = new TH1F("histo2402","M RSgrav rec    ",50, 3.0,4.0);
  TH1F *histo2403 = new TH1F("histo2403","M RSgrav gen    ",100, 3.0,4.0);
  TH1F *histo2404 = new TH1F("histo2404","M RSgrav rec    ",100, 3.0,4.0);
  TH1F *histo2405 = new TH1F("histo2405","M RSgrav gen    ",1000,3.0,4.0);
  TH1F *histo2406 = new TH1F("histo2406","M RSgrav rec    ",1000,3.0,4.0);
  TH1F *histo2407 = new TH1F("histo2407","CosTheta gen    ",20,-1.,1.);
  TH1F *histo2408 = new TH1F("histo2408","CosTheta rec    ",20,-1.,1.);
  TH1F *histo2409 = new TH1F("histo2409","CosTheta gen    ",10,-1.,1.);
  TH1F *histo2410 = new TH1F("histo2410","CosTheta rec    ",10,-1.,1.);
  
  // 4TeV
  TH1F *histo2500 = new TH1F("histo2500","M RSgrav gen ee ",50, 3.5,4.5);
  TH1F *histo2501 = new TH1F("histo2501","M RSgrav gen    ",50, 3.5,4.5);
  TH1F *histo2502 = new TH1F("histo2502","M RSgrav rec    ",50, 3.5,4.5);
  TH1F *histo2503 = new TH1F("histo2503","M RSgrav gen    ",100, 3.5,4.5);
  TH1F *histo2504 = new TH1F("histo2504","M RSgrav rec    ",100, 3.5,4.5);
  TH1F *histo2505 = new TH1F("histo2505","M RSgrav gen    ",1000,3.5,4.5);
  TH1F *histo2506 = new TH1F("histo2506","M RSgrav rec    ",1000,3.5,4.5);
  TH1F *histo2507 = new TH1F("histo2507","CosTheta gen    ",20,-1.,1.);
  TH1F *histo2508 = new TH1F("histo2508","CosTheta rec    ",20,-1.,1.);
  TH1F *histo2509 = new TH1F("histo2509","CosTheta gen    ",10,-1.,1.);
  TH1F *histo2510 = new TH1F("histo2510","CosTheta rec    ",10,-1.,1.);

  // 4TeV sans saturation
  TH1F *histo2600 = new TH1F("histo2600","M RSgrav gen ee ",50, 3.5,4.5);
  TH1F *histo2601 = new TH1F("histo2601","M RSgrav gen    ",50, 3.5,4.5);
  TH1F *histo2602 = new TH1F("histo2602","M RSgrav rec    ",50, 3.5,4.5);
  TH1F *histo2603 = new TH1F("histo2603","M RSgrav gen    ",100, 3.5,4.5);
  TH1F *histo2604 = new TH1F("histo2604","M RSgrav rec    ",100, 3.5,4.5);
  TH1F *histo2605 = new TH1F("histo2605","M RSgrav gen    ",1000,3.5,4.5);
  TH1F *histo2606 = new TH1F("histo2606","M RSgrav rec    ",1000,3.5,4.5);
  TH1F *histo2607 = new TH1F("histo2607","CosTheta gen    ",20,-1.,1.);
  TH1F *histo2608 = new TH1F("histo2608","CosTheta rec    ",20,-1.,1.);
  TH1F *histo2609 = new TH1F("histo2609","CosTheta gen    ",10,-1.,1.);
  TH1F *histo2610 = new TH1F("histo2610","CosTheta rec    ",10,-1.,1.);

  // 4TeV5
  TH1F *histo2700 = new TH1F("histo2700","M RSgrav gen ee ",50, 4.0,5.0);
  TH1F *histo2701 = new TH1F("histo2701","M RSgrav gen    ",50, 4.0,5.0);
  TH1F *histo2702 = new TH1F("histo2702","M RSgrav rec    ",50, 4.0,5.0);
  TH1F *histo2703 = new TH1F("histo2703","M RSgrav gen    ",100, 4.0,5.0);
  TH1F *histo2704 = new TH1F("histo2704","M RSgrav rec    ",100, 4.0,5.0);
  TH1F *histo2705 = new TH1F("histo2705","M RSgrav gen    ",1000,4.0,5.0);
  TH1F *histo2706 = new TH1F("histo2706","M RSgrav rec    ",1000,4.0,5.0);
  TH1F *histo2707 = new TH1F("histo2707","CosTheta gen    ",20,-1.,1.);
  TH1F *histo2708 = new TH1F("histo2708","CosTheta rec    ",20,-1.,1.);
  TH1F *histo2709 = new TH1F("histo2709","CosTheta gen    ",10,-1.,1.);
  TH1F *histo2710 = new TH1F("histo2710","CosTheta rec    ",10,-1.,1.);

  // 4TeV5 sans saturation
  TH1F *histo2800 = new TH1F("histo2800","M RSgrav gen ee ",50, 4.0,5.0);
  TH1F *histo2801 = new TH1F("histo2801","M RSgrav gen    ",50, 4.0,5.0);
  TH1F *histo2802 = new TH1F("histo2802","M RSgrav rec    ",50, 4.0,5.0);
  TH1F *histo2803 = new TH1F("histo2803","M RSgrav gen    ",100, 4.0,5.0);
  TH1F *histo2804 = new TH1F("histo2804","M RSgrav rec    ",100, 4.0,5.0);
  TH1F *histo2805 = new TH1F("histo2805","M RSgrav gen    ",1000,4.0,5.0);
  TH1F *histo2806 = new TH1F("histo2806","M RSgrav rec    ",1000,4.0,5.0);
  TH1F *histo2807 = new TH1F("histo2807","CosTheta gen    ",20,-1.,1.);
  TH1F *histo2808 = new TH1F("histo2808","CosTheta rec    ",20,-1.,1.);
  TH1F *histo2809 = new TH1F("histo2809","CosTheta gen    ",10,-1.,1.);
  TH1F *histo2810 = new TH1F("histo2810","CosTheta rec    ",10,-1.,1.);

  //   Start main loop on all events
  Int_t nb = 0; 
  Int_t nevent = (Int_t)T->GetEntries();

//  nevent =10;

  
  cout <<  "Nombre d entrees a traiter : " << nevent << endl; 
  for (Int_t i=0; i<nevent; i++) {
    nb += T->GetEntry(i); 


       // Monte Carlo information
   
      // generator event
    Short_t mc_np = 0, mc_nsp = 0, mc_nssp = 0, mc_npp = 0, mc_nspp = 0 ;
    mc_np   = generatorEvent->getNumberOfParticles();
    mc_nsp  = generatorEvent->getNumberOfSignalParticles();
    mc_nssp = generatorEvent->getNumberOfStableSignalParticles();
    mc_npp  = generatorEvent->getNumberOfPileupParticles();
    mc_nspp = generatorEvent->getNumberOfStablePileupParticles();
    genpart=generatorEvent->getParticles();
    int npart=0,iRSgrav=11,ilink1=0,ilink2=0;
    for (int gp=0;gp<mc_np;gp++) {
         XANAGeneratorParticle *gpart = (XANAGeneratorParticle *)((*genpart)[gp]);
         if ((gpart->getPartId()==11 || gpart->getPartId()==-11) &&  gpart->getMo1()==12 )
         { 
               if (npart==0)
               {
                    // elec 1
                    npart++;
                    ilink1=gp;
               }
               else if (npart==1)
               {
                    //elec 2
                   npart++;
                   ilink2=gp;
               }
               else 
               {
                   cout << " pas possible d'avoir 3 electrons rattaches au RSgrav " << endl;
               }
         } 
    }
    XANAGeneratorParticle *gRSgrav = (XANAGeneratorParticle *)((*genpart)[iRSgrav]);
    XANAGeneratorParticle *gElec1 = (XANAGeneratorParticle *)((*genpart)[ilink1]);
    XANAGeneratorParticle *gElec2 = (XANAGeneratorParticle *)((*genpart)[ilink2]);

    float trueMassGen=gRSgrav->getMass();
    float eeMassGen=RSMassGen(gElec1,gElec2);
    float costhetaGen = RSdistribAngGen(gElec1,gElec2);

    histo1000->Fill( eeMassGen/1000.);
    histo2000->Fill( eeMassGen/1000.);
    histo2100->Fill( eeMassGen/1000.);
    histo2200->Fill( eeMassGen/1000.);
    histo2300->Fill( eeMassGen/1000.);
    histo2400->Fill( eeMassGen/1000.);
    histo2500->Fill( eeMassGen/1000.);
    histo2600->Fill( eeMassGen/1000.);
    histo2700->Fill( eeMassGen/1000.);
    histo2800->Fill( eeMassGen/1000.);

    histo1001->Fill( trueMassGen/1000.);

    histo2001->Fill( trueMassGen/1000.);
    histo2003->Fill( trueMassGen/1000.);
    histo2005->Fill( trueMassGen/1000.);

    histo2101->Fill( trueMassGen/1000.);
    histo2103->Fill( trueMassGen/1000.);
    histo2105->Fill( trueMassGen/1000.);

    histo2201->Fill( trueMassGen/1000.);
    histo2203->Fill( trueMassGen/1000.);
    histo2205->Fill( trueMassGen/1000.);

    histo2301->Fill( trueMassGen/1000.);
    histo2303->Fill( trueMassGen/1000.);
    histo2305->Fill( trueMassGen/1000.);

    histo2401->Fill( trueMassGen/1000.);
    histo2403->Fill( trueMassGen/1000.);
    histo2405->Fill( trueMassGen/1000.);

    histo2501->Fill( trueMassGen/1000.);
    histo2503->Fill( trueMassGen/1000.);
    histo2505->Fill( trueMassGen/1000.);

    histo2601->Fill( trueMassGen/1000.);
    histo2603->Fill( trueMassGen/1000.);
    histo2605->Fill( trueMassGen/1000.);

    histo2701->Fill( trueMassGen/1000.);
    histo2703->Fill( trueMassGen/1000.);
    histo2705->Fill( trueMassGen/1000.);
    
    histo2801->Fill( trueMassGen/1000.);
    histo2803->Fill( trueMassGen/1000.);
    histo2805->Fill( trueMassGen/1000.);

    histo1004->Fill( costhetaGen );
    histo1006->Fill( costhetaGen );

    // Reconstructed event
    
    // Definitions
    Short_t nc = 0, nh = 0, nhc = 0, npc = 0,  nt = 0, nv = 0, nvt = 0;
    Short_t nm = 0, nelt = 0 ;    
  

    nc  = event->getNumberOfEmClusters();         clusters = event->getEmClusters(); 
    nh  = event->getNumberOfTrackHits();          hits = event->getEmClusterHits();   
    nhc = event->getNumberOfHadClusters();        hadClusters = event->getHadClusters();   
    ns_supercluster  = event->getNumberOfSuperClusters(); 
                                                  superClusters = event->getSuperClusters(); 
    npc = event->getNumberOfPreshClusters();      preshClusters = event->getPreshClusters(); 
    nt  = event->getNumberOfTracks();             tracks = event->getTracks();
    nv  = event->getNumberOfVertices();           vertices = event->getVertices();   
    nvt = event->getNumberOfVertexTracks(); 
    ne_candidatelectron  = event->getNumberOfElectronCandidates(); 
                                                  electrons = event->getElectronCandidates(); 
    nm  = event->getNumberOfMuonCandidates();     muons = event->getMuonCandidates(); 
    nelt    = event->getNumberOfElectronTracks(); electronTracks = event->getElectronTracks();


    
    // electrons
    Int_t indexel1=-1, indexel2=-1;
    
    // Selection des 2 candidats electrons les plus energetiques

    RSSelection(indexel1, indexel2);

    // On ne travaille qu'avec les 2 Candidats Electrons les + energetiques
    
    if (indexel1>-1 && indexel2>-1)
    {
               XANAElectronCandidate *elec1 = (XANAElectronCandidate*)((*electrons)[indexel1]);
               XANAElectronCandidate *elec2 = (XANAElectronCandidate*)((*electrons)[indexel2]);
              
               float energy1 =elec1->getSuperClusterEnergy()/elec1->getSuperCluster()->getEnergyScaleFactor();
               float eta1    =elec1->getSuperCluster()->getPosition().eta();
               float phi1    =elec1->getSuperCluster()->getPosition().phi();
               
               float energy2 =elec2->getSuperClusterEnergy()/elec2->getSuperCluster()->getEnergyScaleFactor();
               float eta2    =elec2->getSuperCluster()->getPosition().eta();
               float phi2    =elec2->getSuperCluster()->getPosition().phi();
              
               // Traitement de la saturation
               int iflag_satur=0;
               RSSaturation(elec1, energy1, iflag_satur);
               RSSaturation(elec2, energy2, iflag_satur);              
               
               // Association des photons
               float EtIsol1=0.;
               float EtIsol2=0.;
               RSAssocPhoton(elec1, energy1, EtIsol1);
               RSAssocPhoton(elec2, energy2, EtIsol2);
               
               // Correction de l'energie pour eta
               
               bool FAMOS=true;
               float xcor1=RSCorrectionFunc(energy1,  eta1,  FAMOS);
               float xcor2=RSCorrectionFunc(energy2,  eta2,  FAMOS);
               xcor1=1./xcor1;
               xcor2=1./xcor2;
               energy1*=xcor1;
               energy2*=xcor2;
               
               // Candidat RSElec
               
               RSElec *cElec1 = new RSElec(elec1,energy1,eta1,phi1,EtIsol1);
               RSElec *cElec2 = new RSElec(elec2,energy2,eta2,phi2,EtIsol2);

                // Masse Reconstruite
                float eeMassRec=RSMassRec(cElec1,cElec2);
                
                // Cut en iso (H/E pas encore dispo dans FAMOS 1_2_0)
                if ( cElec1->getElecIsol()< 0.02* cElec1->getElecEt() &&
                     cElec2->getElecIsol()< 0.02* cElec2->getElecEt() )
                    {
                        float costhetaRec = RSdistribAngRec(cElec1,cElec2);

                        histo1002->Fill( eeMassRec /1000. );
                        histo1003->Fill( trueMassGen-eeMassRec);
                        
                        histo1005->Fill(costhetaRec);
                        histo1007->Fill(costhetaRec);
                        histo1008->Fill( costhetaGen-costhetaRec);
                        histo1009->Fill( costhetaGen-costhetaRec);

                        histo2002->Fill( eeMassRec /1000. );
                        histo2004->Fill( eeMassRec /1000. );
                        histo2006->Fill( eeMassRec /1000. );
                        
                        histo2102->Fill( eeMassRec /1000. );
                        histo2104->Fill( eeMassRec /1000. );
                        histo2106->Fill( eeMassRec /1000. );
                        
                        histo2202->Fill( eeMassRec /1000. );
                        histo2204->Fill( eeMassRec /1000. );
                        histo2206->Fill( eeMassRec /1000. );
                        
                        histo2302->Fill( eeMassRec /1000. );
                        histo2304->Fill( eeMassRec /1000. );
                        histo2306->Fill( eeMassRec /1000. );
                        
                        histo2502->Fill( eeMassRec /1000. );
                        histo2504->Fill( eeMassRec /1000. );
                        histo2506->Fill( eeMassRec /1000. );
                        
                        histo2702->Fill( eeMassRec /1000. );
                        histo2704->Fill( eeMassRec /1000. );
                        histo2706->Fill( eeMassRec /1000. );
                        
                        if (iflag_satur==0)
                        {
                           histo2402->Fill( eeMassRec /1000. );
                           histo2404->Fill( eeMassRec /1000. );
                           histo2406->Fill( eeMassRec /1000. );
                        
                           histo2602->Fill( eeMassRec /1000. );
                           histo2604->Fill( eeMassRec /1000. );
                           histo2606->Fill( eeMassRec /1000. );

                           histo2802->Fill( eeMassRec /1000. );
                           histo2804->Fill( eeMassRec /1000. );
                           histo2806->Fill( eeMassRec /1000. );
                        }

                        if (1900<eeMassRec && eeMassRec<2100)
                        {
                          histo2007->Fill( costhetaGen );
                          histo2008->Fill( costhetaRec );
                          histo2009->Fill( costhetaGen );
                          histo2010->Fill( costhetaRec );
                        }
                        else if (2400<eeMassRec && eeMassRec<2600)
                        {
                          histo2107->Fill( costhetaGen );
                          histo2108->Fill( costhetaRec );
                          histo2109->Fill( costhetaGen );
                          histo2110->Fill( costhetaRec );
                        }
                        else if (2900<eeMassRec && eeMassRec<3100)
                        {
                          histo2207->Fill( costhetaGen );
                          histo2208->Fill( costhetaRec );
                          histo2209->Fill( costhetaGen );
                          histo2210->Fill( costhetaRec );
                        }
                        else if (3400<eeMassRec && eeMassRec<3600)
                        {
                          histo2307->Fill( costhetaGen );
                          histo2308->Fill( costhetaRec );
                          histo2309->Fill( costhetaGen );
                          histo2310->Fill( costhetaRec );
                          if (iflag_satur==0)
                          {
                            histo2407->Fill( costhetaGen );
                            histo2408->Fill( costhetaRec );
                            histo2409->Fill( costhetaGen );
                            histo2410->Fill( costhetaRec );
                          }
                        }
                        else if (3900<eeMassRec && eeMassRec<4100)
                        {
                          histo2507->Fill( costhetaGen );
                          histo2508->Fill( costhetaRec );
                          histo2509->Fill( costhetaGen );
                          histo2510->Fill( costhetaRec );
                          if (iflag_satur==0)
                          {
                            histo2607->Fill( costhetaGen );
                            histo2608->Fill( costhetaRec );
                            histo2609->Fill( costhetaGen );
                            histo2610->Fill( costhetaRec );
                          }
                        }
                        else if (4400<eeMassRec && eeMassRec<4600)
                        {
                          histo2707->Fill( costhetaGen );
                          histo2708->Fill( costhetaRec );
                          histo2709->Fill( costhetaGen );
                          histo2710->Fill( costhetaRec );
                          if (iflag_satur==0)
                          {
                            histo2807->Fill( costhetaGen );
                            histo2808->Fill( costhetaRec );
                            histo2809->Fill( costhetaGen );
                            histo2810->Fill( costhetaRec );
                          }
                        }



                    }



    }


    event->clear(); 
    generatorEvent->clear();
    geantEvent->clear();
  }
  
//  f.Close();
  //
     TFile *filehisto = new TFile("Mass_0.root","RECREATE");

     histo1000->Write();
     histo1001->Write();
     histo1002->Write();
     histo1003->Write();
     histo1004->Write();
     histo1005->Write();
     histo1006->Write();
     histo1007->Write();
     histo1008->Write();
     histo1009->Write();


     histo2000->Write();
     histo2001->Write();
     histo2002->Write();
     histo2003->Write();
     histo2004->Write();
     histo2005->Write();
     histo2006->Write();
     histo2007->Write();
     histo2008->Write();
     histo2009->Write();
     histo2010->Write();

     histo2100->Write();
     histo2101->Write();
     histo2102->Write();
     histo2103->Write();
     histo2104->Write();
     histo2105->Write();
     histo2106->Write();
     histo2107->Write();
     histo2108->Write();
     histo2109->Write();
     histo2110->Write();

     histo2200->Write();
     histo2201->Write();
     histo2202->Write();
     histo2203->Write();
     histo2204->Write();
     histo2205->Write();
     histo2206->Write();
     histo2207->Write();
     histo2208->Write();
     histo2209->Write();
     histo2210->Write();

     histo2300->Write();
     histo2301->Write();
     histo2302->Write();
     histo2303->Write();
     histo2304->Write();
     histo2305->Write();
     histo2306->Write();
     histo2307->Write();
     histo2308->Write();
     histo2309->Write();
     histo2310->Write();

     histo2400->Write();
     histo2401->Write();
     histo2402->Write();
     histo2403->Write();
     histo2404->Write();
     histo2405->Write();
     histo2406->Write();
     histo2407->Write();
     histo2408->Write();
     histo2409->Write();
     histo2410->Write();

     histo2500->Write();
     histo2501->Write();
     histo2502->Write();
     histo2503->Write();
     histo2504->Write();
     histo2505->Write();
     histo2506->Write();
     histo2507->Write();
     histo2508->Write();
     histo2509->Write();
     histo2510->Write();

     histo2600->Write();
     histo2601->Write();
     histo2602->Write();
     histo2603->Write();
     histo2604->Write();
     histo2605->Write();
     histo2606->Write();
     histo2607->Write();
     histo2608->Write();
     histo2609->Write();
     histo2610->Write();

     histo2700->Write();
     histo2701->Write();
     histo2702->Write();
     histo2703->Write();
     histo2704->Write();
     histo2705->Write();
     histo2706->Write();
     histo2707->Write();
     histo2708->Write();
     histo2709->Write();
     histo2710->Write();

     histo2800->Write();
     histo2801->Write();
     histo2802->Write();
     histo2803->Write();
     histo2804->Write();
     histo2805->Write();
     histo2806->Write();
     histo2807->Write();
     histo2808->Write();
     histo2809->Write();
     histo2810->Write();
     
     filehisto->Close();

  
}

---