rsAnalc0p01_o.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/XANAVertex/interface/XANAVertex.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)
{
   
// version TChain
   TChain *T = new TChain("ESD");
   T->Add("/data/scratch/collard/Production_xana151/grav/oscar/oRSee_1500_c001_1.root");
   T->Add("/data/scratch/collard/Production_xana151/grav/oscar/oRSee_1500_c001_2.root");
   T->Add("/data/scratch/collard/Production_xana151/grav/oscar/oRSee_1500_c001_3.root");
   T->Add("/data/scratch/collard/Production_xana151/grav/oscar/oRSee_1500_c001_4.root");
   T->Add("/data/scratch/collard/Production_xana151/grav/oscar/oRSee_1500_c001_5.root");
   T->Add("/data/scratch/collard/Production_xana151/grav/oscar/oRSee_1500_c001_6.root");
   T->Add("/data/scratch/collard/Production_xana151/grav/oscar/oRSee_1500_c001_7.root");
   T->Add("/data/scratch/collard/Production_xana151/grav/oscar/oRSee_1500_c001_8.root");
   T->Add("/data/scratch/collard/Production_xana151/grav/oscar/oRSee_1500_c001_9.root");
   T->Add("/data/scratch/collard/Production_xana151/grav/oscar/oRSee_1500_c001_10.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);
  
  // 0TeV75 
  TH1F *histo2000 = new TH1F("histo2000","M RSgrav gen ee ",50, 0.55,0.85);
  TH1F *histo2001 = new TH1F("histo2001","M RSgrav gen    ",50, 0.55,0.85);
  TH1F *histo2002 = new TH1F("histo2002","M RSgrav rec    ",50, 0.55,0.85);
  TH1F *histo2003 = new TH1F("histo2003","M RSgrav gen    ",100, 0.55,0.85);
  TH1F *histo2004 = new TH1F("histo2004","M RSgrav rec    ",100, 0.55,0.85);
  TH1F *histo2005 = new TH1F("histo2005","M RSgrav gen    ",200,0.55,0.85);
  TH1F *histo2006 = new TH1F("histo2006","M RSgrav rec    ",200,0.55,0.85);
  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.);

  // 1TeV 
  TH1F *histo2100 = new TH1F("histo2100","M RSgrav gen ee ",50, 0.8,1.1);
  TH1F *histo2101 = new TH1F("histo2101","M RSgrav gen    ",50, 0.8,1.1);
  TH1F *histo2102 = new TH1F("histo2102","M RSgrav rec    ",50, 0.8,1.1);
  TH1F *histo2103 = new TH1F("histo2103","M RSgrav gen    ",100,0.8,1.1);
  TH1F *histo2104 = new TH1F("histo2104","M RSgrav rec    ",100,0.8,1.1);
  TH1F *histo2105 = new TH1F("histo2105","M RSgrav gen    ",200,0.8,1.1);
  TH1F *histo2106 = new TH1F("histo2106","M RSgrav rec    ",200,0.8,1.1);
  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.);

  // 1TeV25 
  TH1F *histo2200 = new TH1F("histo2200","M RSgrav gen ee ",50, 1.05,1.35);
  TH1F *histo2201 = new TH1F("histo2201","M RSgrav gen    ",50, 1.05,1.35);
  TH1F *histo2202 = new TH1F("histo2202","M RSgrav rec    ",50, 1.05,1.35);
  TH1F *histo2203 = new TH1F("histo2203","M RSgrav gen    ",100,1.05,1.35);
  TH1F *histo2204 = new TH1F("histo2204","M RSgrav rec    ",100,1.05,1.35);
  TH1F *histo2205 = new TH1F("histo2205","M RSgrav gen    ",200,1.05,1.35);
  TH1F *histo2206 = new TH1F("histo2206","M RSgrav rec    ",200,1.05,1.35);
  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.);

  // 1TeV5 
  TH1F *histo2300 = new TH1F("histo2300","M RSgrav gen ee ",50, 1.3,1.6);
  TH1F *histo2301 = new TH1F("histo2301","M RSgrav gen    ",50, 1.3,1.6);
  TH1F *histo2302 = new TH1F("histo2302","M RSgrav rec    ",50, 1.3,1.6);
  TH1F *histo2303 = new TH1F("histo2303","M RSgrav gen    ",100,1.3,1.6);
  TH1F *histo2304 = new TH1F("histo2304","M RSgrav rec    ",100,1.3,1.6);
  TH1F *histo2305 = new TH1F("histo2305","M RSgrav gen    ",200,1.3,1.6);
  TH1F *histo2306 = new TH1F("histo2306","M RSgrav rec    ",200,1.3,1.6);
  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.);

  // 1TeV75 
  TH1F *histo2400 = new TH1F("histo2400","M RSgrav gen ee ",50, 1.55,1.85);
  TH1F *histo2401 = new TH1F("histo2401","M RSgrav gen    ",50, 1.55,1.85);
  TH1F *histo2402 = new TH1F("histo2402","M RSgrav rec    ",50, 1.55,1.85);
  TH1F *histo2403 = new TH1F("histo2403","M RSgrav gen    ",100,1.55,1.85);
  TH1F *histo2404 = new TH1F("histo2404","M RSgrav rec    ",100,1.55,1.85);
  TH1F *histo2405 = new TH1F("histo2405","M RSgrav gen    ",200,1.55,1.85);
  TH1F *histo2406 = new TH1F("histo2406","M RSgrav rec    ",200,1.55,1.85);
  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.);

  // 2TeV 
  TH1F *histo2500 = new TH1F("histo2500","M RSgrav gen ee ",50, 1.8,2.1);
  TH1F *histo2501 = new TH1F("histo2501","M RSgrav gen    ",50, 1.8,2.1);
  TH1F *histo2502 = new TH1F("histo2502","M RSgrav rec    ",50, 1.8,2.1);
  TH1F *histo2503 = new TH1F("histo2503","M RSgrav gen    ",100,1.8,2.1);
  TH1F *histo2504 = new TH1F("histo2504","M RSgrav rec    ",100,1.8,2.1);
  TH1F *histo2505 = new TH1F("histo2505","M RSgrav gen    ",200,1.8,2.1);
  TH1F *histo2506 = new TH1F("histo2506","M RSgrav rec    ",200,1.8,2.1);
  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.);

  // 2TeV5 
  TH1F *histo2600 = new TH1F("histo2600","M RSgrav gen ee ",50, 2.3,2.6);
  TH1F *histo2601 = new TH1F("histo2601","M RSgrav gen    ",50, 2.3,2.6);
  TH1F *histo2602 = new TH1F("histo2602","M RSgrav rec    ",50, 2.3,2.6);
  TH1F *histo2603 = new TH1F("histo2603","M RSgrav gen    ",100,2.3,2.6);
  TH1F *histo2604 = new TH1F("histo2604","M RSgrav rec    ",100,2.3,2.6);
  TH1F *histo2605 = new TH1F("histo2605","M RSgrav gen    ",200,2.3,2.6);
  TH1F *histo2606 = new TH1F("histo2606","M RSgrav rec    ",200,2.3,2.6);
  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.);


  //   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.);

    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.);
    

    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();


      // z vertex 
    float zvertex=0;
    if (nv>0)
    {
     XANAVertex *vertexlink = (XANAVertex*)((*vertices)[0]);
     zvertex=vertexlink->getPosition().z();
    }
    
   
    
    // 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=false;
               float xcor1=RSCorrectionFunc(energy1,  eta1,  FAMOS);
               float xcor2=RSCorrectionFunc(energy2,  eta2,  FAMOS);
               xcor1=1./xcor1;
               xcor2=1./xcor2;
               energy1*=xcor1;
               energy2*=xcor2;

               // Correction de la position du vertex
               RSZvertexCorrection(zvertex,eta1);
               RSZvertexCorrection(zvertex,eta2);
               
               // 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. );
                        
                        histo2402->Fill( eeMassRec /1000. );
                        histo2404->Fill( eeMassRec /1000. );
                        histo2406->Fill( eeMassRec /1000. );
                        
                        histo2502->Fill( eeMassRec /1000. );
                        histo2504->Fill( eeMassRec /1000. );
                        histo2506->Fill( eeMassRec /1000. );
                        
                        histo2602->Fill( eeMassRec /1000. );
                        histo2604->Fill( eeMassRec /1000. );
                        histo2606->Fill( eeMassRec /1000. );
                        
                        if (725<eeMassRec && eeMassRec<775)
                        {
                          histo2007->Fill( costhetaGen );
                          histo2008->Fill( costhetaRec );
                          histo2009->Fill( costhetaGen );
                          histo2010->Fill( costhetaRec );
                        }
                        else if (975<eeMassRec && eeMassRec<1025)
                        {
                          histo2107->Fill( costhetaGen );
                          histo2108->Fill( costhetaRec );
                          histo2109->Fill( costhetaGen );
                          histo2110->Fill( costhetaRec );
                        }
                        else if (1225<eeMassRec && eeMassRec<1275)
                        {
                          histo2207->Fill( costhetaGen );
                          histo2208->Fill( costhetaRec );
                          histo2209->Fill( costhetaGen );
                          histo2210->Fill( costhetaRec );
                        }
                        else if (1475<eeMassRec && eeMassRec<1525)
                        {
                          histo2307->Fill( costhetaGen );
                          histo2308->Fill( costhetaRec );
                          histo2309->Fill( costhetaGen );
                          histo2310->Fill( costhetaRec );
                        }
                        else if (1725<eeMassRec && eeMassRec<1775)
                        {
                          histo2407->Fill( costhetaGen );
                          histo2408->Fill( costhetaRec );
                          histo2409->Fill( costhetaGen );
                          histo2410->Fill( costhetaRec );
                        }
                        else if (1975<eeMassRec && eeMassRec<2025)
                        {
                          histo2507->Fill( costhetaGen );
                          histo2508->Fill( costhetaRec );
                          histo2509->Fill( costhetaGen );
                          histo2510->Fill( costhetaRec );
                        }
                        else if (2475<eeMassRec && eeMassRec<2525)
                        {
                          histo2607->Fill( costhetaGen );
                          histo2608->Fill( costhetaRec );
                          histo2609->Fill( costhetaGen );
                          histo2610->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();

     
     filehisto->Close();

  
}

---