Index: /trunk/Resolutions.cpp
===================================================================
--- /trunk/Resolutions.cpp	(revision 446)
+++ /trunk/Resolutions.cpp	(revision 447)
@@ -33,5 +33,5 @@
 
 /// \file Resolution.cpp
-/// \brief Prepares the resolution calculation
+/// \brief Resolution for CMS
 
 #include "TChain.h"
@@ -111,5 +111,5 @@
 
 
-void PairingJet(TLorentzVector &JETSm, const TLorentzVector &JET, const TClonesArray *branchJet, const float dR=0.25)
+void PairingJet(TLorentzVector &JETSm, const TLorentzVector &JET, const TClonesArray *branchJet)
 {
   JETSm.SetPxPyPzE(0,0,0,0);
@@ -125,5 +125,5 @@
         {
           deltaRtest = DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta());
-          if(deltaRtest < dR)
+          if(deltaRtest < 0.25)
              {
                JETSm = Att;
@@ -203,5 +203,5 @@
   if(argc != 3) {
       cout << " Usage: " << argv[0] << " input_file" << " output_file" << endl;
-      cout << " input_file - input file in Delphes-root format," << endl;
+      cout << " input_list - list of files in root format," << endl;
       cout << " output_file - output file." << endl;
       exit(1);
@@ -220,5 +220,5 @@
 
 
-  TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE");// Creates the file, but should be closed just after
+  TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
   outputFile->Close();
   
@@ -257,17 +257,9 @@
   
   RESOLution *DET = new RESOLution();
-  /*
-  string detectorcard = "data/DetectorCard_CMS.dat";
-  const float dR_jetpairing = 0.25;
-  const float jet_pt_cut = 1;
-  */
-  string detectorcard = "data/DetectorCard_ATLAS.dat";
-  const float dR_jetpairing = 0.2;
-  const float jet_pt_cut = 7;
-  DET->ReadDataCard(detectorcard);
+  DET->ReadDataCard("data/DetectorCard_CMS.dat");
 
 
   //Jet information
-  JetsUtil *JETRUN = new JetsUtil(detectorcard);
+  JetsUtil *JETRUN = new JetsUtil("data/DetectorCard_CMS.dat");
 
   TLorentzVector genMomentum(0,0,0,0);//TLorentzVector containing generator level information
@@ -277,5 +269,7 @@
   vector<fastjet::PseudoJet> input_particlesGEN;//for FastJet algorithm
   vector<fastjet::PseudoJet> sorted_jetsGEN;
+ 
   vector<int> NTrackJet;
+ 
   vector<TLorentzVector> towers;
  
@@ -315,6 +309,5 @@
 	  //Calculate ETMIS from generated particles
           if((pid == pNU1) || (pid == pNU2) || (pid == pNU3))PTmisGEN = PTmisGEN + genMomentum;
-
-	  //Electrons and muons	 
+	 
           if( (particle->Status == 1)    &&
             ((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
@@ -335,8 +328,8 @@
                  case pMU: // all muons with eta < DET->MAX_MU
                     PairingMuon(recoMomentum,genMomentum,branchMuon);
-                    if(recoMomentum.E()!=0){
+                    if(recoMomentum.E() !=0){
                      elementMuon = (RESOLMUON*) branchmuon->NewEntry();
-                     elementMuon->OverPT = 1./genMomentum.Pt();
-                     elementMuon->OverSmearedPT = 1./recoMomentum.Pt();}
+                     elementMuon->OverPT = (1/genMomentum.Pt());
+                     elementMuon->OverSmearedPT = (1/recoMomentum.Pt());}
                    break; // case pMU
                  default:
@@ -369,4 +362,5 @@
       elementEtmis->Ex = (-PTmisGEN).Px();
       elementEtmis->SEt = ScalarEt;
+
       elementEtmis->EtSmeare = (PTmisReco).Pt()-(PTmisGEN).Pt();
       elementEtmis->ExSmeare = (-PTmisReco).Px()-(PTmisGEN).Px();
@@ -376,4 +370,5 @@
 
       TSimpleArray<TRootGenParticle> TausHadr = TauHadr(branchGen);
+
       TLorentzVector JETreco(0,0,0,0);
       for (unsigned int i = 0; i < sorted_jetsGEN.size(); i++) {
@@ -381,12 +376,9 @@
         JETgen.SetPxPyPzE(sorted_jetsGEN[i].px(),sorted_jetsGEN[i].py(),sorted_jetsGEN[i].pz(),sorted_jetsGEN[i].E());
         PairingJet(JETreco,JETgen,branchJet);
-        if(JETreco.Pt()>jet_pt_cut)
+        if(JETreco.Pt()>1)
           {
             elementJet= (RESOLJET*) branchjet->NewEntry();
             elementJet->PT = JETgen.Et();
             elementJet->SmearedPT = JETreco.Et()/JETgen.Et();
-            elementJet->E = JETgen.E();
-            elementJet->dE = (JETreco.E()-JETgen.E())/JETgen.E() ;
-            elementJet->dE2 = pow( (JETreco.E()-JETgen.E())/JETgen.E() , 2.) ;
           }
       }
@@ -407,5 +399,5 @@
                   {
                     elementTaujet= (TAUHAD*) branchtaujet->NewEntry();
-                    elementTaujet->EnergieCen = EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E();
+                    elementTaujet->EnergieCen = (EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E());
                     elementTaujet->NumTrack = NumTracks(branchTracks,DET->TAU_track_pt,JETT.Eta(),JETT.Phi(),DET->TAU_track_scone);
                     if( (EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E()) > 0.95 
@@ -417,11 +409,16 @@
 	
       } // for itJet : loop on all jets
+//cout<<"i"<<endl;
       
       treeWriter->Fill();
     } // Loop over all events
   treeWriter->Write();
-
-  cout << detectorcard << " has been used.\n";
+float frac = numTauRec/numTau;  
+cout<<numTauRec<<endl;
+cout<<numTau<<endl;
+
   cout << "** Exiting..." << endl;
+  cout<<frac<<endl;
+
   
   delete treeWriter;
Index: /trunk/Resolutions_ATLAS.cpp
===================================================================
--- /trunk/Resolutions_ATLAS.cpp	(revision 447)
+++ /trunk/Resolutions_ATLAS.cpp	(revision 447)
@@ -0,0 +1,431 @@
+/***********************************************************************
+**                                                                    **
+**          /----------------------------------------------\          **
+**         |  Delphes, a framework for the fast simulation  |         **
+**         |       of a  generic collider experiment        |         **
+**          \-------------  arXiv:0903.2225v1  ------------/          **
+**                                                                    **
+**                                                                    **
+**   This package uses:                                               **
+**   ------------------                                               **
+**   ROOT: Nucl. Inst. & Meth. in Phys. Res. A389 (1997) 81-86        **
+**   FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210]      **
+**   Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2]       **
+**   FROG: [hep-ex/0901.2718v1]                                       **
+**   HepMC: Comput. Phys. Commun.134 (2001) 41                        **
+**                                                                    **
+** ------------------------------------------------------------------ **
+**                                                                    **
+**   Main authors:                                                    **
+**   -------------                                                    **
+**                                                                    **
+**                Severine Ovyn                Xavier Rouby           **
+**          severine.ovyn@uclouvain.be      xavier.rouby@cern         **
+**                                                                    **
+**         Center for Particle Physics and Phenomenology (CP3)        **
+**                Universite catholique de Louvain (UCL)              **
+**                     Louvain-la-Neuve, Belgium                      **
+**                                                                    **
+**                      Copyright (C) 2008-2009,                      **
+**                        All rights reserved.                        **
+**                                                                    **
+***********************************************************************/
+
+/// \file Resolution.cpp
+/// \brief Prepares the resolution calculation
+
+#include "TChain.h"
+#include "TApplication.h"
+#include "TFile.h"
+
+#include "ExRootTreeReader.h"
+#include "ExRootTreeWriter.h"
+#include "ExRootTreeBranch.h"
+#include "TreeClasses.h"
+
+#include "DataConverter.h"
+#include "HEPEVTConverter.h"
+#include "LHEFConverter.h"
+#include "STDHEPConverter.h"
+
+#include "SmearUtil.h"
+#include "JetsUtil.h"
+#include "BFieldProp.h"
+
+//#include "PseudoJet.hh"
+//#include "ClusterSequence.hh"
+
+#include<vector>
+#include<iostream>
+
+using namespace std;
+
+//------------------------------------------------------------------------------
+
+// //********************************** PYTHIA INFORMATION*********************************
+
+TSimpleArray<TRootGenParticle> TauHadr(const TClonesArray *GEN)
+ {
+   TIter it((TCollection*)GEN);
+   it.Reset();
+   TRootGenParticle *gen1;
+   TSimpleArray<TRootGenParticle> array,array2;
+
+   while((gen1 = (TRootGenParticle*) it.Next()))
+      {
+        array.Add(gen1);
+      }
+   it.Reset();
+   bool tauhad;
+   while((gen1 = (TRootGenParticle*) it.Next()))
+    {
+       tauhad=false;
+       if(abs(gen1->PID)==15)
+         {
+            int d1=gen1->D1;
+            int d2=gen1->D2;
+            if((d1 < array.GetEntries()) && (d1 > 0) && (d2 < array.GetEntries()) && (d2 > 0))
+               {
+                 tauhad=true;
+                 for(int d=d1; d < d2+1; d++)
+                    {
+                      if(abs(array[d]->PID)== pE || abs(array[d]->PID)== pMU)tauhad=false;
+                    }
+               }
+         }
+        if(tauhad)array2.Add(gen1);
+     }
+   return array2;
+ }
+
+double EnergySmallCone(const vector<TLorentzVector> &towers, const float eta, const float phi,float energy_scone,float JET_seed) {
+  double Energie=0;
+  for(unsigned int i=0; i < towers.size(); i++) {
+      if(towers[i].Pt() < JET_seed) continue;
+      if((DeltaR(phi,eta,towers[i].Phi(),towers[i].Eta()) < energy_scone)) {
+          Energie += towers[i].E();
+      }
+  }
+  return Energie;
+}
+
+
+void PairingJet(TLorentzVector &JETSm, const TLorentzVector &JET, const TClonesArray *branchJet, const float dR=0.25)
+{
+  JETSm.SetPxPyPzE(0,0,0,0);
+  float deltaRtest=5000;
+  TIter itJet((TCollection*)branchJet);
+  TRootJet *jet;
+  itJet.Reset();
+  while( (jet = (TRootJet*) itJet.Next()) )
+   {
+      TLorentzVector Att;
+      Att.SetPtEtaPhiE(jet->PT,jet->Eta,jet->Phi,jet->E);
+      if(DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
+        {
+          deltaRtest = DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta());
+          if(deltaRtest < dR)
+             {
+               JETSm = Att;
+             }
+        }
+    }
+}
+
+void PairingElec(TLorentzVector &ELECSm, const TLorentzVector &ELEC, const TClonesArray *branchElec)
+{
+  ELECSm.SetPxPyPzE(0,0,0,0);
+  float deltaRtest=5000;
+  TIter itElec((TCollection*)branchElec);
+  TRootElectron *elec;
+  itElec.Reset();
+  while( (elec = (TRootElectron*) itElec.Next()) )
+   {
+      TLorentzVector Att;
+      Att.SetPtEtaPhiE(elec->PT,elec->Eta,elec->Phi,elec->E);
+      if(DeltaR(ELEC.Phi(),ELEC.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
+        {
+          deltaRtest = DeltaR(ELEC.Phi(),ELEC.Eta(),Att.Phi(),Att.Eta());
+          if(deltaRtest < 0.025)
+             {
+               ELECSm = Att;
+             }
+        }
+    }
+}
+
+void PairingMuon(TLorentzVector &MUONSm, const TLorentzVector &MUON, const TClonesArray *branchMuon)
+{
+  MUONSm.SetPxPyPzE(0,0,0,0);
+  float deltaRtest=5000;
+  TIter itMuon((TCollection*)branchMuon);
+  TRootMuon *muon;
+  itMuon.Reset();
+  while( (muon = (TRootMuon*) itMuon.Next()) )
+   {
+      TLorentzVector Att;
+      Att.SetPxPyPzE(muon->Px,muon->Py,muon->Pz,muon->E);
+      if(DeltaR(MUON.Phi(),MUON.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
+        {
+          deltaRtest = DeltaR(MUON.Phi(),MUON.Eta(),Att.Phi(),Att.Eta());
+          if(deltaRtest < 0.025)
+             {
+               MUONSm = Att;
+             }
+        }
+    }
+}
+
+unsigned int NumTracks(const TClonesArray *branchTracks, const float pt_track, const float eta, const float phi,float track_scone) {
+  unsigned int numtrack=0;
+  TIter itTrack((TCollection*)branchTracks);
+  TRootTracks *track;
+  itTrack.Reset();
+  while( (track = (TRootTracks*) itTrack.Next()) )
+   {
+      if((track->PT < pt_track )||
+         (DeltaR(phi,eta,track->Phi,track->Eta) > track_scone)
+        )continue;
+      numtrack++;
+  }
+  return numtrack;
+}
+
+
+
+int main(int argc, char *argv[])
+{
+  int appargc = 2;
+  char *appName = "Resolution";
+  char *appargv[] = {appName, "-b"};
+  TApplication app(appName, &appargc, appargv);
+  
+  if(argc != 3) {
+      cout << " Usage: " << argv[0] << " input_file" << " output_file" << endl;
+      cout << " input_file - input file in Delphes-root format," << endl;
+      cout << " output_file - output file." << endl;
+      exit(1);
+  }
+
+  srand (time (NULL));         /* Initialisation du gÃ©nÃ©rateur */
+  
+  //read the input TROOT file
+  string inputfilename(argv[1]), outputfilename(argv[2]);
+
+  if(outputfilename.find(".root") > outputfilename.length() ) {
+    cout << "output_file should be a .root file!\n";
+    return -1;
+  }
+
+
+
+  TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE");// Creates the file, but should be closed just after
+  outputFile->Close();
+  
+  TChain chainGEN("GEN");
+  chainGEN.Add(inputfilename.c_str());
+  ExRootTreeReader *treeReaderGEN = new ExRootTreeReader(&chainGEN);
+  TChain chain("Analysis");
+  chain.Add(inputfilename.c_str());
+  ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
+  const TClonesArray *branchJet = treeReader->UseBranch("Jet");
+  const TClonesArray *branchElec = treeReader->UseBranch("Electron");
+  const TClonesArray *branchMuon = treeReader->UseBranch("Muon");
+  const TClonesArray *branchTracks = treeReader->UseBranch("Tracks");
+  const TClonesArray *branchTowers = treeReader->UseBranch("CaloTower");
+  const TClonesArray *branchGen = treeReaderGEN->UseBranch("Particle");
+  TIter itGen((TCollection*)branchGen);
+  
+  //write the output root file
+  ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputfilename, "Analysis");
+  ExRootTreeBranch *branchjet = treeWriter->NewBranch("JetPTResol", RESOLJET::Class());
+  ExRootTreeBranch *branchelec = treeWriter->NewBranch("ElecEResol", RESOLELEC::Class());
+  ExRootTreeBranch *branchmuon = treeWriter->NewBranch("MuonPTResol", RESOLMUON::Class());
+  ExRootTreeBranch *branchtaujet = treeWriter->NewBranch("TauJetPTResol", TAUHAD::Class());
+  ExRootTreeBranch *branchetmis = treeWriter->NewBranch("ETmisResol",ETMIS::Class());
+  
+  TRootGenParticle *particle;
+  
+  RESOLELEC * elementElec;
+  RESOLMUON *elementMuon;
+  RESOLJET *elementJet;
+  TAUHAD *elementTaujet;
+  ETMIS *elementEtmis;
+  
+  int numTau=0;
+  int numTauRec=0;
+  
+  RESOLution *DET = new RESOLution();
+  /*
+  string detectorcard = "data/DetectorCard_CMS.dat";
+  const float dR_jetpairing = 0.25;
+  const float jet_pt_cut = 1;
+  */
+  string detectorcard = "data/DetectorCard_ATLAS.dat";
+  const float dR_jetpairing = 0.2;
+  const float jet_pt_cut = 7;
+  DET->ReadDataCard(detectorcard);
+
+
+  //Jet information
+  JetsUtil *JETRUN = new JetsUtil(detectorcard);
+
+  TLorentzVector genMomentum(0,0,0,0);//TLorentzVector containing generator level information
+  TLorentzVector recoMomentum(0,0,0,0);//TLorentzVector containing generator level information
+  LorentzVector jetMomentum;
+  
+  vector<fastjet::PseudoJet> input_particlesGEN;//for FastJet algorithm
+  vector<fastjet::PseudoJet> sorted_jetsGEN;
+  vector<int> NTrackJet;
+  vector<TLorentzVector> towers;
+ 
+  // Loop over all events
+  Long64_t entry, allEntries = treeReader->GetEntries();
+  cout << "** Chain contains " << allEntries << " events" << endl;
+  for(entry = 0; entry < allEntries; ++entry)
+    {
+      TLorentzVector PTmisReco(0,0,0,0);
+      TLorentzVector PTmisGEN(0,0,0,0);
+      treeReader->ReadEntry(entry);
+      treeReaderGEN->ReadEntry(entry);
+      treeWriter->Clear();
+      if((entry % 100) == 0 && entry > 0 )  cout << "** Processing element # " << entry << endl;
+      
+      TSimpleArray<TRootGenParticle> bGen;
+      itGen.Reset();
+      TSimpleArray<TRootGenParticle> NFCentralQ;
+      
+      input_particlesGEN.clear();
+      towers.clear();
+      
+      // Loop over all particles in event
+      while( (particle = (TRootGenParticle*) itGen.Next()) )
+	{
+          genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
+	  
+	  int pid = abs(particle->PID);
+	  float eta = fabs(particle->Eta);
+	  
+          //input generator level particle for jet algorithm	  
+          if(particle->Status == 1 && eta < DET->CEN_max_calo_fwd)
+	    {
+	      input_particlesGEN.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
+	    }
+	  
+	  //Calculate ETMIS from generated particles
+          if((pid == pNU1) || (pid == pNU2) || (pid == pNU3))PTmisGEN = PTmisGEN + genMomentum;
+
+	  //Electrons and muons	 
+          if( (particle->Status == 1)    &&
+            ((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
+             (fabs(particle->Eta) < DET->CEN_max_calo_fwd)
+            )
+            {
+              eta=fabs(genMomentum.Eta());
+
+              switch(pid) {
+
+                 case pE: // all electrons with eta < DET->MAX_CALO_FWD
+                    PairingElec(recoMomentum,genMomentum,branchElec);
+                    if(recoMomentum.E()!=0){
+                     elementElec=(RESOLELEC*) branchelec->NewEntry();
+                     elementElec->E = genMomentum.E();
+                     elementElec->SmearedE = recoMomentum.E();}
+                   break; // case pE
+                 case pMU: // all muons with eta < DET->MAX_MU
+                    PairingMuon(recoMomentum,genMomentum,branchMuon);
+                    if(recoMomentum.E()!=0){
+                     elementMuon = (RESOLMUON*) branchmuon->NewEntry();
+                     elementMuon->OverPT = 1./genMomentum.Pt();
+                     elementMuon->OverSmearedPT = 1./recoMomentum.Pt();}
+                   break; // case pMU
+                 default:
+                   break;
+              } // switch (pid)
+           }
+
+	} // while
+     
+      //compute missing transverse energy from calo towers    
+      TIter itCalo((TCollection*)branchTowers);
+      TRootCalo *calo;
+      itCalo.Reset();
+      TLorentzVector Att(0.,0.,0.,0.);
+      float ScalarEt=0;
+      while( (calo = (TRootCalo*) itCalo.Next()) )
+           {
+             if(calo->E !=0){
+             Att.SetPtEtaPhiE(calo->getET(),calo->Eta,calo->Phi,calo->E);
+             towers.push_back(Att);
+             if(fabs(Att.Eta()) < DET->CEN_max_calo_fwd)
+               {
+                  ScalarEt = ScalarEt + calo->getET();
+                  PTmisReco = PTmisReco + Att;
+               }
+           }  
+      }
+      elementEtmis= (ETMIS*) branchetmis->NewEntry();
+      elementEtmis->Et = (PTmisGEN).Pt();
+      elementEtmis->Ex = (-PTmisGEN).Px();
+      elementEtmis->SEt = ScalarEt;
+      elementEtmis->EtSmeare = (PTmisReco).Pt()-(PTmisGEN).Pt();
+      elementEtmis->ExSmeare = (-PTmisReco).Px()-(PTmisGEN).Px();
+
+      //*****************************
+      sorted_jetsGEN=JETRUN->RunJetsResol(input_particlesGEN);
+
+      TSimpleArray<TRootGenParticle> TausHadr = TauHadr(branchGen);
+      TLorentzVector JETreco(0,0,0,0);
+      for (unsigned int i = 0; i < sorted_jetsGEN.size(); i++) {
+        TLorentzVector JETgen(0,0,0,0);
+        JETgen.SetPxPyPzE(sorted_jetsGEN[i].px(),sorted_jetsGEN[i].py(),sorted_jetsGEN[i].pz(),sorted_jetsGEN[i].E());
+        PairingJet(JETreco,JETgen,branchJet);
+        if(JETreco.Pt()>jet_pt_cut)
+          {
+            elementJet= (RESOLJET*) branchjet->NewEntry();
+            elementJet->PT = JETgen.Et();
+            elementJet->SmearedPT = JETreco.Et()/JETgen.Et();
+            elementJet->E = JETgen.E();
+            elementJet->dE = (JETreco.E()-JETgen.E())/JETgen.E() ;
+            elementJet->dE2 = pow( (JETreco.E()-JETgen.E())/JETgen.E() , 2.) ;
+          }
+      }
+      numTau = numTau+TausHadr.GetEntries();
+
+      TIter itJet((TCollection*)branchJet);
+      TRootJet *jet;
+      itJet.Reset();
+      while( (jet = (TRootJet*) itJet.Next()) )
+        {
+          TLorentzVector JETT(0,0,0,0);
+          JETT.SetPxPyPzE(jet->Px,jet->Py,jet->Pz,jet->E);
+          if(fabs(JETT.Eta()) < (DET->CEN_max_tracker - DET->TAU_track_scone)) 
+           {
+            for(Int_t i=0; i<TausHadr.GetEntries();i++)
+              {
+                if(DeltaR(TausHadr[i]->Phi,TausHadr[i]->Eta,JETT.Phi(),JETT.Eta())<0.1)
+                  {
+                    elementTaujet= (TAUHAD*) branchtaujet->NewEntry();
+                    elementTaujet->EnergieCen = EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E();
+                    elementTaujet->NumTrack = NumTracks(branchTracks,DET->TAU_track_pt,JETT.Eta(),JETT.Phi(),DET->TAU_track_scone);
+                    if( (EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E()) > 0.95 
+                     &&  (NumTracks(branchTracks,DET->TAU_track_pt,JETT.Eta(),JETT.Phi(),DET->TAU_track_scone))==1)numTauRec++;
+                 }
+              }
+           }
+         
+	
+      } // for itJet : loop on all jets
+      
+      treeWriter->Fill();
+    } // Loop over all events
+  treeWriter->Write();
+
+  cout << detectorcard << " has been used.\n";
+  cout << "** Exiting..." << endl;
+  
+  delete treeWriter;
+  delete treeReader;
+  delete DET;
+}
+
