Index: trunk/Delphes.cpp
===================================================================
--- trunk/Delphes.cpp	(revision 260)
+++ trunk/Delphes.cpp	(revision 264)
@@ -31,5 +31,5 @@
 
 /// \file Delphes.cpp
-/// \brief executable for the Delphes
+/// \brief Executable for Delphes
 
 #include "TChain.h"
@@ -41,4 +41,5 @@
 #include "ExRootTreeWriter.h"
 #include "ExRootTreeBranch.h"
+#include "ExRootProgressBar.h"
 
 #include "DataConverter.h"
@@ -48,4 +49,5 @@
 
 #include "SmearUtil.h"
+#include "CaloUtil.h"
 #include "BFieldProp.h"
 #include "TriggerUtil.h"
@@ -56,6 +58,4 @@
 #include <vector>
 #include <iostream>
-
-#include "Utilities/ExRootAnalysis/interface/ExRootProgressBar.h"
 
 using namespace std;
@@ -77,4 +77,6 @@
 int main(int argc, char *argv[])
 {
+unsigned int nnnn=0, mmmm=0;
+
   int appargc = 2;
   char *appName= new char[20];
@@ -86,5 +88,4 @@
   delete [] appName;
   delete [] appOpt;
-
 
   if(argc != 3 && argc != 4 && argc != 5) {
@@ -242,13 +243,10 @@
   //Propagation of tracks in the B field
   TrackPropagation *TRACP = new TrackPropagation(DET); 
-  //TrackPropagation *TRACP = new TrackPropagation(DetDatacard); 
   
   //Jet information
   JetsUtil *JETRUN = new JetsUtil(DET);
-  //JetsUtil *JETRUN = new JetsUtil(DetDatacard);
   
   //VFD information
   VeryForward * VFD = new VeryForward(DET);
-  //VeryForward * VFD = new VeryForward(DetDatacard);
  
   // data converters
@@ -299,5 +297,6 @@
   ExRootTreeBranch *branchMuon = treeWriter->NewBranch("Muon", TRootMuon::Class());
   ExRootTreeBranch *branchPhoton = treeWriter->NewBranch("Photon", TRootPhoton::Class());
-  ExRootTreeBranch *branchTracks = treeWriter->NewBranch("Tracks", TRootTracks::Class());
+  //ExRootTreeBranch *branchTracks = treeWriter->NewBranch("Tracks", TRootTracks::Class());
+  ExRootTreeBranch *branchTrack = treeWriter->NewBranch("Tracks", D_Track::Class());
   ExRootTreeBranch *branchETmis = treeWriter->NewBranch("ETmis", TRootETmis::Class());
   ExRootTreeBranch *branchCalo = treeWriter->NewBranch("CaloTower", TRootCalo::Class());
@@ -306,10 +305,9 @@
   ExRootTreeBranch *branchFP420 = treeWriter->NewBranch("FP420hits", TRootRomanPotHits::Class());
   
-  TRootGenParticle *particle; 
   TRootETmis *elementEtmis;
   TRootElectron *elementElec;
   TRootMuon *elementMu;
   TRootPhoton *elementPhoton;
-  TRootTracks *elementTracks;
+  D_Track * elementTrack;
   TRootCalo *elementCalo;
   
@@ -323,11 +321,32 @@
   vector<TLorentzVector> TrackCentral;  
   vector<PhysicsTower> towers;
-  vector<ParticleUtil> electron;
-  vector<ParticleUtil> muon;
-  vector<ParticleUtil> gamma;
+  //vector<ParticleUtil> electron;
+  vector<D_Particle> electron;
+  //vector<ParticleUtil> muon;
+  vector<D_Particle> muon;
+  //vector<ParticleUtil> gamma;
+  vector<D_Particle> gamma;
 
   TSimpleArray<TRootGenParticle> NFCentralQ;
-  float iPhi=0,iEta=0;
-
+
+D_CaloList list_of_calorimeters;
+D_CaloElement CentralCalo("centralcalo", 
+	-DET->CEN_max_calo_cen, DET->CEN_max_calo_cen,
+	 DET->ELG_Ccen,  DET->ELG_Ncen,  DET->ELG_Scen,
+         DET->HAD_Chcal, DET->HAD_Nhcal, DET->HAD_Shcal);
+D_CaloElement ForwardCalo("forwardcalo",
+	 DET->CEN_max_calo_cen, DET->CEN_max_calo_fwd, 
+         DET->ELG_Cfwd,  DET->ELG_Nfwd,  DET->ELG_Sfwd,
+         DET->HAD_Chf,   DET->HAD_Nhf,   DET->HAD_Shf  );
+D_CaloElement BackwardCalo("backwardcalo",
+	-DET->CEN_max_calo_fwd, -DET->CEN_max_calo_cen,
+         DET->ELG_Cfwd,  DET->ELG_Nfwd,  DET->ELG_Sfwd,
+         DET->HAD_Chf,   DET->HAD_Nhf,   DET->HAD_Shf  );
+//D_CaloElement CastorCalo("castor",5.5,6.6,1,0,0,1,0,0);
+list_of_calorimeters.addElement(CentralCalo);
+list_of_calorimeters.addElement(ForwardCalo);
+list_of_calorimeters.addElement(BackwardCalo);
+//list_of_calorimeters.addElement(CastorCalo);
+list_of_calorimeters.sortElements();
 
 
@@ -360,9 +379,13 @@
       towers.clear();
       input_particles.clear();
-      
+
+
+/*
+if(0) { // OLD SMEARING ALGORITHM without energy flow
       // 2.1 Loop over all particles in event
       itGen.Reset();
       while( (particle = (TRootGenParticle*) itGen.Next()) ) {
           int pid = abs(particle->PID);
+          float iPhi=0,iEta=0;
 
 	
@@ -381,5 +404,5 @@
 	  // the ordering of conditions have been optimised for speed : put first the STATUS condition
           if( (particle->Status == 1)    &&
-	      ((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
+	      (pid != pNU1) && (pid != pNU2) && (pid != pNU3) &&
               (fabs(particle->Eta) < DET->CEN_max_calo_fwd)
 	      ) 
@@ -431,9 +454,9 @@
 	      if(genMomentumBfield.E() !=0 && pid != pMU) {
 		// in case the Bfield is not simulated, checks that charged particles have enough pt to reach the calos
-		if ( !DET->FLAG_bfield && charge!=0 && genMomentumBfield.Pt() <= DET->TRACK_ptmin ) { /* particules do not reach calos */ }
+		if ( !DET->FLAG_bfield && charge!=0 && genMomentumBfield.Pt() <= DET->TRACK_ptmin ) {}// particules do not reach calos 
 		else { // particles reach calos	
 	  	  // applies the calo segmentation and returns iEta & iPhi
 		  DET->BinEtaPhi(genMomentumBfield.Phi(), genMomentumBfield.Eta(), iPhi, iEta); 
-                  if(iEta != -100 && iPhi != -100) { 
+                  if(iEta != UNDEFINED && iPhi != UNDEFINED) { 
                       momentumCaloSegmentation.SetPtEtaPhiE(genMomentumBfield.Pt(),iEta,iPhi,genMomentumBfield.E());
                       elementCalo = (TRootCalo*) branchCalo->NewEntry();
@@ -441,8 +464,8 @@
 		      PhysicsTower Tower(LorentzVector(momentumCaloSegmentation.Px(),momentumCaloSegmentation.Py(),momentumCaloSegmentation.Pz(),momentumCaloSegmentation.E()));
                       towers.push_back(Tower);
-                  } // if iEta != -100
+                  } // if iEta != UNDEFINED
 		} // else : when particles reach the calos
 	      } // 2.1.2.3 calotowers
-	      
+
 
 	      // 2.1.2.4 ********************* central detector: tracks
@@ -451,5 +474,5 @@
 		 (genMomentumBfield.E()!=0) &&
 		 (fabs(genMomentumBfield.Eta()) < DET->CEN_max_tracker) && 
-		 (DET->FLAG_bfield || ( !DET->FLAG_bfield && genMomentumBfield.Pt() > DET->TRACK_ptmin )) &&     
+		 (DET->FLAG_bfield || ( !DET->FLAG_bfield && genMomentum.Pt() > DET->TRACK_ptmin )) &&     
 			// if bfield not simulated, pt should be high enough to be taken into account
 		 ((rand()%100) < DET->TRACK_eff)  &&
@@ -475,11 +498,199 @@
 	     VFD->RomanPots(treeWriter,branchRP220,branchFP420,particle);
 	  }
-	  
-       } // 2.1 while : loop on all particles of the event.
-     
+} // IF OLD ALGORITHM (= no energy flow )
+
+
+
+
+else // IF NEW ALGORITHM with energy flow              
+*/
+{
+      D_CaloTowerList list_of_active_towers; 
+      D_CaloTowerList list_of_towers_with_photon; // to speed up the code: will only look in interesting towers for gamma candidates
+
+      // 2.1a Loop over all particles in event, to fill the towers
+      itGen.Reset();
+      TRootGenParticle *particle; 
+      while( (particle = (TRootGenParticle*) itGen.Next()) ) {
+	  int pid = abs(particle->PID);
+	  particle->Charge=Charge(pid);
+	  particle->setFractions(); // init
+
+
+	  // 2.1a.1********************* preparation for the b-tagging
+	  //// This subarray is needed for the B-jet algorithm
+	  // optimization for speed : put first PID condition, then ETA condition, then either pt or status
+	  if( (pid <= pB || pid == pGLUON) &&// is it a light quark or a gluon, i.e. is it one of these : u,d,c,s,b,g ?
+              fabs(particle->Eta) < DET->CEN_max_tracker &&
+              particle->Status != 1 &&
+              particle->PT > DET->PT_QUARKS_MIN ) {
+            NFCentralQ.Add(particle);
+          }
+
+	  // 2.1a.2********************* visible particles only
+	  if( (particle->Status == 1) && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) ){
+
+            // 2.1a.2.1 Central solenoidal magnetic field 
+	    TRACP->bfield(particle); // fills in particle->EtaCalo et particle->PhiCalo
+	    particle->SetEtaPhi(particle->EtaCalo,particle->PhiCalo); // ???? check this line
+
+	    // 2.1a.2.2 Filling the calorimetric towers -- includes also forward detectors ?
+	    // first checks if the charged particles reach the calo!
+	    if( DET->FLAG_bfield || 
+                particle->getCharge()==0 ||
+                (!DET->FLAG_bfield && particle->getCharge()!=0 && particle->PT > DET->TRACK_ptmin)) 
+	    if(
+		(particle->EtaCalo > list_of_calorimeters.getEtamin() ) &&
+		(particle->EtaCalo < list_of_calorimeters.getEtamax() )
+              ) {
+		float iEta=UNDEFINED, iPhi=UNDEFINED;
+		DET->BinEtaPhi(particle->PhiCalo,particle->EtaCalo,iPhi,iEta); // fills in iPhi and iEta
+		if (iEta != UNDEFINED && iPhi != UNDEFINED) 
+                {
+		  D_CaloTower tower(iEta,iPhi);	                               // new tower
+		  tower.Set_Eem_Ehad_E_ET(particle->E*particle->getFem() , particle->E*particle->getFhad() );
+		
+		  list_of_active_towers.addTower(tower);
+		  // this list may contain several times the same calotower, as several particles 
+		  // may leave some energy in the same calotower
+		  // After the loop on particles, identical cells in the list should be merged 
+		} // iEta and iPhi must be defined
+	    }
+
+ 	    // 2.1a.2.3 charged particles in tracker: energy flow
+	    // if bfield not simulated, pt should be high enough to be taken into account
+	    // it is supposed here that DET->MAX_calo > DET->CEN_max_tracker > DET->CEN_max_mu > 0
+	    if( particle->getCharge() !=0 &&
+	        fabs(particle->EtaCalo)< DET->CEN_max_tracker && // stays in the tracker -> track available
+	       ( DET->FLAG_bfield || 
+                  (!DET->FLAG_bfield && particle->PT > DET->TRACK_ptmin)
+                )
+              ) {	
+		// 2.1a.2.3.1 Filling the particle properties + smearing
+		   // Hypothesis: the final eta/phi are the ones from the generator, thanks to the track reconstruction
+		   // This is the EnergyFlow hypothesis
+		   particle->SetEtaPhi(particle->Eta,particle->Phi);
+		   float sET=UNDEFINED; // smeared ET, computed from the smeared E -> needed for the tracks
+
+		   // 2.1a.2.3.2 Muons
+		   if (pid == pMU && fabs(particle->EtaCalo)< DET->CEN_max_mu && particle->PT > DET->PTCUT_muon ) {
+                        TLorentzVector p;
+                        float sPT = gRandom->Gaus(particle->PT, DET->MU_SmearPt*particle->PT );
+			if (sPT > 0 && sPT > DET->PTCUT_muon) {
+                           p.SetPtEtaPhiE(sPT,particle->Eta,particle->Phi,sPT*cosh(particle->Eta));
+                           muon.push_back(D_Particle(p,pMU,particle->EtaCalo,particle->PhiCalo));
+                        }
+			sET = (sPT >0)? sPT : 0;
+                   }  
+		   // 2.1a.2.3.3 Electrons
+		   else if (pid == pE) {
+		     // Finds in which calorimeter the particle has gone, to know its resolution
+
+		     D_CaloElement currentCalo = list_of_calorimeters.getElement(particle->EtaCalo);
+		     if(currentCalo.getName() == dummyCalo.getName()) {
+			cout << "** Warning: the calo coverage behind the tracker is not complete! **" << endl; }
+
+		     // final smeared EM energy   // electromagnetic fraction F_em =1 for electrons;
+	             float sE  = currentCalo.getElectromagneticResolution().Smear(particle->E);
+		     if (sE>0) {
+ 		       sET = sE/cosh(particle->Eta); 
+		       // NB: ET is found via the calorimetry and not via the track curvature
+
+		       TLorentzVector p;
+		       p.SetPtEtaPhiE(sET,particle->Eta,particle->Phi,sE);
+		       if (sET > DET->PTCUT_elec)
+		         electron.push_back(D_Particle(p,particle->PID,particle->EtaCalo,particle->PhiCalo));
+                     } else { sET=0;} // if negative smeared energy -- needed for the tracks
+		   }  
+		   // 2.1a.2.3.4 Other charged particles : smear them for the tracks!
+		   else {  //other particles
+			D_CaloElement currentCalo = list_of_calorimeters.getElement(particle->EtaCalo);
+			float sEem  = currentCalo.getElectromagneticResolution().Smear(particle->E * particle->getFem());
+			float sEhad = currentCalo.getHadronicResolution().Smear(particle->E * particle->getFhad());
+			float sE = ( (sEem>0)? sEem : 0 ) + ( (sEhad>0)? sEhad : 0 );
+			sET = sE/cosh(particle->EtaCalo);
+		   }
+
+		   // 2.1a.2.3.5 Tracks
+		   if( (rand()%100) < DET->TRACK_eff && sET!=0) {
+   		     elementTrack = (D_Track*) branchTrack->NewEntry();
+		     elementTrack->Set(particle->Eta, particle->Phi, particle->EtaCalo, particle->PhiCalo, sET);
+		     TrackCentral.push_back(elementTrack->GetFourVector()); // tracks at vertex!
+		     // TODO!!! associates the tracks to the calo where it points to
+		     // TODO!!! apply a smearing on the position of the origin of the track
+		     // TODO!!! elementTracks->SetPositionOut(Xout,Yout,Zout);
+		   }
+            } // 2.1a.2.3 : if tracker/energy-flow
+ 	    // 2.1a.2.4 Photons
+               // stays in the tracker -> track available -> gamma ID
+	    else if( (pid == pGAMMA) && fabs(particle->EtaCalo)< DET->CEN_max_tracker ) {
+		float iEta=UNDEFINED, iPhi=UNDEFINED;
+                DET->BinEtaPhi(particle->PhiCalo,particle->EtaCalo,iPhi,iEta); // fills in iPhi and iEta
+                D_CaloTower tower(iEta,iPhi); 
+		// stores the list of towers where to apply the photon ID algorithm. Just a trick for a faster search
+		list_of_towers_with_photon.addTower(tower); 
+	    }
+	    // 2.1a.2.5 : very forward detectors
+	    else if (DET->FLAG_vfd==1) {
+		// for the moment, only protons are transported
+		// BUT !!! could be a beam of other particles! (heavy ions?)
+		// BUT ALSO !!! if very forward muons, or others!
+              VFD->RomanPots(treeWriter,branchRP220,branchFP420,particle);
+              VFD->ZDC(treeWriter,branchZDC,particle);
+	    }
+	    // 2.1a.2.6: Zero degree calorimeter
+	    //else if(DET->FLAG_zdc==1) {
+              //VFD->ZDC(treeWriter,branchZDC,particle);
+	    //}
+
+          } // 2.1a.2 : if visible particle
+      } // loop on all particles 2.1a
+	
+
+      // 2.1b loop on all (activated) towers
+      // at this stage, list_of_active_towers may contain several times the same tower
+      // first step is to merge identical towers, by matching their (iEta,iPhi)  
+      list_of_active_towers.mergeDuplicates();
+      // Calotower smearing 
+      list_of_active_towers.smearTowers(list_of_calorimeters);
+
+      for(unsigned int i=0; i<list_of_active_towers.size(); i++) {
+	float iEta = list_of_active_towers[i].getEta();
+	float iPhi = list_of_active_towers[i].getPhi();
+	float e = list_of_active_towers[i].getE(); 
+	if(iEta != UNDEFINED && iPhi != UNDEFINED && e!=0) {
+           elementCalo = (TRootCalo*) branchCalo->NewEntry();
+	   elementCalo->set(list_of_active_towers[i]);
+           // not beautiful : should be improved!
+	   TLorentzVector p;
+	   p.SetPtEtaPhiE(list_of_active_towers[i].getET(), iEta, iPhi, e );
+           PhysicsTower Tower(LorentzVector(p.Px(),p.Py(),p.Pz(),p.E()));
+           towers.push_back(Tower);
+	}
+      } // loop on towers
+
+      // 2.1c photon ID
+      // list_of_towers_with_photon is the list of towers with photon candidates
+      // already smeared !
+      // sorts the vector and smears duplicates 
+      list_of_towers_with_photon.mergeDuplicates(); 
+
+      for(unsigned int i=0; i<list_of_towers_with_photon.size(); i++) {
+	  float eta = list_of_towers_with_photon[i].getEta(); 
+	  float phi = list_of_towers_with_photon[i].getPhi();
+	  D_CaloTower cal(list_of_active_towers.getElement(eta,phi));
+	  if(cal.getEta() != UNDEFINED && cal.getPhi() != UNDEFINED &&  cal.getE() > 0) {
+  	    TLorentzVector p;
+            p.SetPtEtaPhiE(cal.getET(), eta,phi,cal.getE() );
+	    if (cal.getET() > DET->PTCUT_gamma) { gamma.push_back(D_Particle(p,pGAMMA,p.Eta(),p.Phi())); }
+          }
+      } // for -- list of photons
+
+}  // IF NEW ALGORITHM with energy flow
+
+
 
 
       // 2.2 ********** Output preparation & complex objects  ***********
- 
       // 2.2.1 ********************* sorting collections by decreasing pt
       DET->SortedVector(electron); 
@@ -487,4 +698,6 @@
         elementElec = (TRootElectron*) branchElectron->NewEntry();
         elementElec->Set(electron[i].Px(),electron[i].Py(),electron[i].Pz(),electron[i].E());
+        elementElec->EtaCalo = electron[i].EtaCalo(); 
+        elementElec->PhiCalo = electron[i].PhiCalo();
         elementElec->Charge = sign(electron[i].PID());
         elementElec->IsolFlag = DET->Isolation(electron[i].Phi(),electron[i].Eta(),TrackCentral,2.0);//isolation based on tracks
@@ -495,4 +708,6 @@
         elementMu->Charge = sign(muon[i].PID());
         elementMu->Set(muon[i].Px(),muon[i].Py(),muon[i].Pz(),muon[i].E());
+        elementMu->EtaCalo = muon[i].EtaCalo(); 
+        elementMu->PhiCalo = muon[i].PhiCalo();
         elementMu->IsolFlag = DET->Isolation(muon[i].Phi(),muon[i].Eta(),TrackCentral,2.0);
       }
@@ -501,4 +716,6 @@
 	elementPhoton = (TRootPhoton*) branchPhoton->NewEntry();
 	elementPhoton->Set(gamma[i].Px(),gamma[i].Py(),gamma[i].Pz(),gamma[i].E());
+        elementPhoton->EtaCalo = gamma[i].EtaCalo();
+        elementPhoton->PhiCalo = gamma[i].PhiCalo();
       }
       
@@ -522,5 +739,5 @@
       elementEtmis->Py = (-PTmis).Py();
       
-      // 2.2.3 ************* B-tag, tau jets
+      // 2.2.3 ************* jets, B-tag, tau jets
       sorted_jets=JETRUN->RunJets(input_particles);
       JETRUN->RunJetBtagging(treeWriter, branchJet,sorted_jets,NFCentralQ);
@@ -529,5 +746,6 @@
       treeWriter->Fill();
     } // 2. Loop over all events ('for' loop)
- 
+
+
   cout <<"**                                                                 **"<< endl;
   cout <<"**                 Exiting detector simulation...                  **"<< endl;
@@ -649,5 +867,5 @@
   delete treeReader;
   delete DET;
-//  delete TRIGT;
+  delete TRIGT;
   delete TRACP;
   delete JETRUN;
