Fork me on GitHub

Changeset df636c74 in git


Ignore:
Timestamp:
Oct 8, 2015, 4:39:04 PM (9 years ago)
Author:
Pavel Demin <pavel.demin@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
d67b86b
Parents:
9da65a5
Message:

update eflow in SimpleCalorimeter

Location:
modules
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • modules/SimpleCalorimeter.cc

    r9da65a5 rdf636c74  
    5656SimpleCalorimeter::SimpleCalorimeter() :
    5757  fResolutionFormula(0),
    58   fItParticleInputArray(0), fItTrackInputArray(0),
    59   fTowerTrackArray(0), fItTowerTrackArray(0)
    60 {
     58  fItParticleInputArray(0), fItTrackInputArray(0)
     59{
     60  Int_t i;
     61
    6162  fResolutionFormula = new DelphesFormula;
    6263
    63   fTowerTrackArray = new TObjArray;
    64   fItTowerTrackArray = fTowerTrackArray->MakeIterator();
     64  for(i = 0; i < 2; ++i)
     65  {
     66    fTowerTrackArray[i] = new TObjArray;
     67    fItTowerTrackArray[i] = fTowerTrackArray[i]->MakeIterator();
     68  }
    6569}
    6670
     
    6973SimpleCalorimeter::~SimpleCalorimeter()
    7074{
     75  Int_t i;
     76
    7177  if(fResolutionFormula) delete fResolutionFormula;
    7278
    73   if(fTowerTrackArray) delete fTowerTrackArray;
    74   if(fItTowerTrackArray) delete fItTowerTrackArray;
     79  for(i = 0; i < 2; ++i)
     80  {
     81    if(fTowerTrackArray[i]) delete fTowerTrackArray[i];
     82    if(fItTowerTrackArray[i]) delete fItTowerTrackArray[i];
     83  }
    7584}
    7685
     
    137146  }
    138147
    139 /*
    140   TFractionMap::iterator itFractionMap;
    141   for(itFractionMap = fFractionMap.begin(); itFractionMap != fFractionMap.end(); ++itFractionMap)
    142   {
    143     cout << itFractionMap->first << "   " << itFractionMap->second.first  << "   " << itFractionMap->second.second << endl;
    144   }
    145 */
    146 
    147148  // read min E value for towers to be saved
    148149  fEnergyMin = GetDouble("EnergyMin", 0.0);
     
    168169  // create output arrays
    169170  fTowerOutputArray = ExportArray(GetString("TowerOutputArray", "towers"));
    170  
     171
    171172  fEFlowTrackOutputArray = ExportArray(GetString("EFlowTrackOutputArray", "eflowTracks"));
    172173  fEFlowTowerOutputArray = ExportArray(GetString("EFlowTowerOutputArray", "eflowTowers"));
     
    197198  Double_t fraction;
    198199  Double_t energy;
     200  Double_t sigma;
    199201  Int_t pdgCode;
    200202
     
    337339
    338340      fTowerEnergy = 0.0;
    339       fTrackEnergy = 0.0;
     341
     342      fTrackEnergy[0] = 0.0;
     343      fTrackEnergy[1] = 0.0;
    340344
    341345      fTowerTime = 0.0;
     
    347351      fTowerPhotonHits = 0;
    348352
    349       fTowerTrackArray->Clear();
     353      fTowerTrackArray[0]->Clear();
     354      fTowerTrackArray[1]->Clear();
    350355    }
    351356
     
    361366      energy = momentum.E() * fTrackFractions[number];
    362367
    363       fTrackEnergy += energy;
    364 
    365368      fTrackTime += TMath::Sqrt(energy)*position.T();
    366369      fTrackTimeWeight += TMath::Sqrt(energy);
    367370
    368       fTowerTrackArray->Add(track);
     371      if(fTrackFractions[number] > 1.0E-9)
     372      {
     373        sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
     374        if(sigma/momentum.E() < track->TrackResolution)
     375        {
     376          fTrackEnergy[0] += energy;
     377          fTowerTrackArray[0]->Add(track);
     378        }
     379        else
     380        {
     381          fTrackEnergy[1] += energy;
     382          fTowerTrackArray[1]->Add(track);
     383        }
     384      }
     385      else if(fTrackFractions[number] < 1.0E-9)
     386      {
     387        fEFlowTrackOutputArray->Add(track);
     388      }
    369389
    370390      continue;
     
    397417void SimpleCalorimeter::FinalizeTower()
    398418{
    399   Candidate *tower, *track;
     419  Candidate *tower, *track, *mother;
    400420  Double_t energy, pt, eta, phi;
    401421  Double_t sigma;
    402422  Double_t time;
    403423
    404   Double_t trkSigma, fraction;
    405  
    406   Int_t pdgCode;
    407424  TLorentzVector momentum;
    408425  TFractionMap::iterator itFractionMap;
    409  
     426
    410427  if(!fTower) return;
    411428
     
    437454
    438455  fTower->Eem = (!fIsEcal) ? 0 : energy;
    439   fTower->Ehad = (fIsEcal) ? 0 : energy; 
     456  fTower->Ehad = (fIsEcal) ? 0 : energy;
    440457
    441458  fTower->Edges[0] = fTowerEdges[0];
     
    447464  if(energy > 0.0) fTowerOutputArray->Add(fTower);
    448465
    449 
    450 
    451   // fill e-flow candidates
    452   fItTowerTrackArray->Reset();
    453   while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
    454   {
    455      momentum = track->Momentum;
    456    
    457      pdgCode = TMath::Abs(track->PID);
    458 
    459      itFractionMap = fFractionMap.find(pdgCode);
    460      if(itFractionMap == fFractionMap.end())
    461      {
    462        itFractionMap = fFractionMap.find(0);
    463      }
    464 
    465      fraction = itFractionMap->second;
    466    
    467      // charged particle has to deposit either in ECAL or HCAL
    468      if(fraction > 1.0E-9)
    469      {
    470        trkSigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());       
    471        if(track->TrackResolution < trkSigma/momentum.E())
    472        {
    473          energy -= momentum.E();
    474          fEFlowTrackOutputArray->Add(track);
    475        }
    476      }
    477      //forward all tracks from ECAL to HCAL
    478      else if(fIsEcal)
    479      {
    480        fEFlowTrackOutputArray->Add(track);
    481      }
    482      //store muons from HCAL
    483      else if(pdgCode == 13)
    484      {
    485        fEFlowTrackOutputArray->Add(track);
    486      }
    487   }
    488 
     466  // fill e-flow candidates
     467
     468  energy -= fTrackEnergy[1];
     469
     470  fItTowerTrackArray[0]->Reset();
     471  while((track = static_cast<Candidate*>(fItTowerTrackArray[0]->Next())))
     472  {
     473    mother = track;
     474    track = static_cast<Candidate*>(track->Clone());
     475    track->AddCandidate(mother);
     476
     477    track->Momentum *= energy/fTrackEnergy[0];
     478
     479    fEFlowTrackOutputArray->Add(track);
     480  }
     481
     482  fItTowerTrackArray[1]->Reset();
     483  while((track = static_cast<Candidate*>(fItTowerTrackArray[1]->Next())))
     484  {
     485    mother = track;
     486    track = static_cast<Candidate*>(track->Clone());
     487    track->AddCandidate(mother);
     488
     489    fEFlowTrackOutputArray->Add(track);
     490  }
     491
     492  if(fTowerTrackArray[0]->GetEntriesFast() > 0) energy = 0.0;
    489493
    490494  sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
     
    499503
    500504    tower->Eem = (!fIsEcal) ? 0 : energy;
    501     tower->Ehad = (fIsEcal) ? 0 : energy; 
     505    tower->Ehad = (fIsEcal) ? 0 : energy;
    502506
    503507    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
  • modules/SimpleCalorimeter.h

    r9da65a5 rdf636c74  
    5959  Double_t fTowerEta, fTowerPhi, fTowerEdges[4];
    6060  Double_t fTowerEnergy;
    61   Double_t fTrackEnergy;
     61  Double_t fTrackEnergy[2];
    6262
    6363  Double_t fTowerTime;
     
    7575  Bool_t fSmearTowerCenter;
    7676
    77   Bool_t fIsEcal; //! 
     77  Bool_t fIsEcal; //!
    7878
    7979  TFractionMap fFractionMap; //!
     
    101101  TObjArray *fEFlowTrackOutputArray; //!
    102102  TObjArray *fEFlowTowerOutputArray; //!
    103  
    104   TObjArray *fTowerTrackArray; //!
    105   TIterator *fItTowerTrackArray; //!
     103
     104  TObjArray *fTowerTrackArray[2]; //!
     105  TIterator *fItTowerTrackArray[2]; //!
    106106
    107107  void FinalizeTower();
Note: See TracChangeset for help on using the changeset viewer.