Fork me on GitHub

Changeset 82575a3 in git for modules/Calorimeter.cc


Ignore:
Timestamp:
Apr 16, 2014, 5:17:35 PM (11 years ago)
Author:
mselvaggi <mselvaggi@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
ab16005
Parents:
9687203
Message:

added simplecalo

File:
1 edited

Legend:

Unmodified
Added
Removed
  • modules/Calorimeter.cc

    r9687203 r82575a3  
    4141
    4242Calorimeter::Calorimeter() :
    43   fResolutionFormula(0),
     43  fECalResolutionFormula(0), fHCalResolutionFormula(0),
    4444  fItParticleInputArray(0), fItTrackInputArray(0),
    4545  fTowerTrackArray(0), fItTowerTrackArray(0)
    4646{
    47   fResolutionFormula = new DelphesFormula;
    48  
     47  fECalResolutionFormula = new DelphesFormula;
     48  fHCalResolutionFormula = new DelphesFormula;
     49
    4950  fTowerTrackArray = new TObjArray;
    5051  fItTowerTrackArray = fTowerTrackArray->MakeIterator();
     
    5556Calorimeter::~Calorimeter()
    5657{
    57   if(fResolutionFormula) delete fResolutionFormula;
    58  
     58  if(fECalResolutionFormula) delete fECalResolutionFormula;
     59  if(fHCalResolutionFormula) delete fHCalResolutionFormula;
     60
    5961  if(fTowerTrackArray) delete fTowerTrackArray;
    6062  if(fItTowerTrackArray) delete fItTowerTrackArray;
     
    6668{
    6769  ExRootConfParam param, paramEtaBins, paramPhiBins, paramFractions;
    68   Long_t i, j, k, size, sizeEtaBins, sizePhiBins;
    69   Double_t fraction;
     70  Long_t i, j, k, size, sizeEtaBins, sizePhiBins, sizeFractions;
     71  Double_t ecalFraction, hcalFraction;
    7072  TBinMap::iterator itEtaBin;
    7173  set< Double_t >::iterator itPhiBin;
     
    114116  // set default energy fractions values
    115117  fFractionMap.clear();
    116   fFractionMap[0] = 1.0;
     118  fFractionMap[0] = make_pair(0.0, 1.0);
    117119
    118120  for(i = 0; i < size/2; ++i)
    119121  {
    120122    paramFractions = param[i*2 + 1];
    121     fraction = paramFractions[0].GetDouble();
    122     fFractionMap[param[i*2].GetInt()] = fraction;
     123    sizeFractions = paramFractions.GetSize();
     124
     125    ecalFraction = paramFractions[0].GetDouble();
     126    hcalFraction = paramFractions[1].GetDouble();
     127
     128    fFractionMap[param[i*2].GetInt()] = make_pair(ecalFraction, hcalFraction);
    123129  }
    124130/*
     
    130136*/
    131137  // read resolution formulas
    132   fResolutionFormula->Compile(GetString("ResolutionFormula", "0"));
    133  
     138  fECalResolutionFormula->Compile(GetString("ECalResolutionFormula", "0"));
     139  fHCalResolutionFormula->Compile(GetString("HCalResolutionFormula", "0"));
     140
    134141  // import array with output from other modules
    135142  fParticleInputArray = ImportArray(GetString("ParticleInputArray", "ParticlePropagator/particles"));
     
    141148  // create output arrays
    142149  fTowerOutputArray = ExportArray(GetString("TowerOutputArray", "towers"));
    143   fEFlowTowerOutputArray = ExportArray(GetString("EFlowTowerOutputArray", "eflowTowers"));
    144  
     150  fPhotonOutputArray = ExportArray(GetString("PhotonOutputArray", "photons"));
     151 
     152  fEFlowTrackOutputArray = ExportArray(GetString("EFlowTrackOutputArray", "eflowTracks"));
     153  fEFlowPhotonOutputArray = ExportArray(GetString("EFlowPhotonOutputArray", "eflowPhotons"));
     154  fEFlowNeutralHadronOutputArray = ExportArray(GetString("EFlowNeutralHadronOutputArray", "eflowNeutralHadrons"));
     155
     156
    145157}
    146158
     
    167179  Int_t number;
    168180  Long64_t towerHit, towerEtaPhi, hitEtaPhi;
    169   Double_t fraction;
    170   Double_t energy;
     181  Double_t ecalFraction, hcalFraction;
     182  Double_t ecalEnergy, hcalEnergy;
    171183  Int_t pdgCode;
    172184
     
    181193  DelphesFactory *factory = GetFactory();
    182194  fTowerHits.clear();
    183   fTowerFractions.clear();
    184   fTrackFractions.clear();
    185  
     195  fTowerECalFractions.clear();
     196  fTowerHCalFractions.clear();
     197  fTrackECalFractions.clear();
     198  fTrackHCalFractions.clear();
     199
    186200  // loop over all particles
    187201  fItParticleInputArray->Reset();
     
    200214    }
    201215
    202     fraction = itFractionMap->second;
    203     fTowerFractions.push_back(fraction);
    204    
    205     if(fraction < 1.0E-9) continue;
     216    ecalFraction = itFractionMap->second.first;
     217    hcalFraction = itFractionMap->second.second;
     218
     219    fTowerECalFractions.push_back(ecalFraction);
     220    fTowerHCalFractions.push_back(hcalFraction);
     221
     222    if(ecalFraction < 1.0E-9 && hcalFraction < 1.0E-9) continue;
    206223
    207224    // find eta bin [1, fEtaBins.size - 1]
     
    243260    }
    244261
    245     fraction = itFractionMap->second;
    246  
    247     fTrackFractions.push_back(fraction);
    248  
     262    ecalFraction = itFractionMap->second.first;
     263    hcalFraction = itFractionMap->second.second;
     264
     265    fTrackECalFractions.push_back(ecalFraction);
     266    fTrackHCalFractions.push_back(hcalFraction);
     267
    249268    // find eta bin [1, fEtaBins.size - 1]
    250269    itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta());
     
    308327      fTowerEdges[3] = (*phiBins)[phiBin];
    309328
    310       fTowerEnergy = 0.0;
    311       fTrackEnergy = 0.0;
    312      
    313       fTowerTime = 0.0;
    314       fTrackTime = 0.0;
     329      fTowerECalEnergy = 0.0;
     330      fTowerHCalEnergy = 0.0;
     331
     332      fTrackECalEnergy = 0.0;
     333      fTrackHCalEnergy = 0.0;
     334
     335      fTowerECalTime = 0.0;
     336      fTowerHCalTime = 0.0;
     337
     338      fTrackECalTime = 0.0;
     339      fTrackHCalTime = 0.0;
     340
     341      fTowerECalWeightTime = 0.0;
     342      fTowerHCalWeightTime = 0.0;
    315343     
    316       fTowerWeightTime = 0.0;
    317      
    318344      fTowerTrackHits = 0;
    319345      fTowerPhotonHits = 0;
     
    330356      momentum = track->Momentum;
    331357      position = track->Position;
    332  
    333       energy = momentum.E() * fTrackFractions[number];
     358
    334359     
    335       fTrackEnergy += energy;
     360      ecalEnergy = momentum.E() * fTrackECalFractions[number];
     361      hcalEnergy = momentum.E() * fTrackHCalFractions[number];
     362
     363      fTrackECalEnergy += ecalEnergy;
     364      fTrackHCalEnergy += hcalEnergy;
    336365     
    337       fTrackTime += TMath::Sqrt(energy)*position.T();
    338       fTrackWeightTime += TMath::Sqrt(energy);
    339    
     366      fTrackECalTime += TMath::Sqrt(ecalEnergy)*position.T();
     367      fTrackHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
     368       
     369      fTrackECalWeightTime += TMath::Sqrt(ecalEnergy);
     370      fTrackHCalWeightTime += TMath::Sqrt(hcalEnergy);
     371
    340372      fTowerTrackArray->Add(track);
    341373
     
    351383
    352384    // fill current tower
    353     energy = momentum.E() * fTowerFractions[number];
    354    
    355     fTowerEnergy += energy;
     385    ecalEnergy = momentum.E() * fTowerECalFractions[number];
     386    hcalEnergy = momentum.E() * fTowerHCalFractions[number];
     387
     388    fTowerECalEnergy += ecalEnergy;
     389    fTowerHCalEnergy += hcalEnergy;
     390
     391    fTowerECalTime += TMath::Sqrt(ecalEnergy)*position.T();
     392    fTowerHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
     393
     394    fTowerECalWeightTime += TMath::Sqrt(ecalEnergy);
     395    fTowerHCalWeightTime += TMath::Sqrt(hcalEnergy);
    356396   
    357     fTowerTime += TMath::Sqrt(energy)*position.T();
    358     fTowerWeightTime += TMath::Sqrt(energy);
    359    
     397
    360398    fTower->AddCandidate(particle);
    361399  }
     
    369407void Calorimeter::FinalizeTower()
    370408{
    371   Candidate *tower;
     409  Candidate *track, *tower;
    372410  Double_t energy, pt, eta, phi;
    373   Double_t sigma;
    374   Double_t time;
     411  Double_t ecalEnergy, hcalEnergy;
     412  Double_t ecalSigma, hcalSigma;
     413  Double_t ecalTime, hcalTime, time;
    375414
    376415  if(!fTower) return;
    377416
    378   sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerEnergy);
    379 
    380 //  energy = gRandom->Gaus(fTowerEnergy, sigma);
    381 //  if(energy < 0.0) energy = 0.0;
    382 
    383   energy = LogNormal(fTowerEnergy, sigma);
    384   time = (fTowerWeightTime < 1.0E-09 ) ? 0 : fTowerTime/fTowerWeightTime;
     417  ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerECalEnergy);
     418
     419//  ecalEnergy = gRandom->Gaus(fTowerECalEnergy, ecalSigma);
     420//  if(ecalEnergy < 0.0) ecalEnergy = 0.0;
     421
     422  ecalEnergy = LogNormal(fTowerECalEnergy, ecalSigma);
     423  ecalTime = (fTowerECalWeightTime < 1.0E-09 ) ? 0 : fTowerECalTime/fTowerECalWeightTime;
     424
     425  hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy);
     426
     427//  hcalEnergy = gRandom->Gaus(fTowerHCalEnergy, hcalSigma);
     428//  if(hcalEnergy < 0.0) hcalEnergy = 0.0;
     429
     430  hcalEnergy = LogNormal(fTowerHCalEnergy, hcalSigma);
     431  hcalTime = (fTowerHCalWeightTime < 1.0E-09 ) ? 0 : fTowerHCalTime/fTowerHCalWeightTime;
     432
     433  energy = ecalEnergy + hcalEnergy;
     434  time = (TMath::Sqrt(ecalEnergy)*ecalTime + TMath::Sqrt(hcalEnergy)*hcalTime)/(TMath::Sqrt(ecalEnergy) + TMath::Sqrt(hcalEnergy));
     435
     436//  eta = fTowerEta;
     437//  phi = fTowerPhi;
    385438
    386439  eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
     
    392445  fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time);
    393446  fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
    394  
     447  fTower->Eem = ecalEnergy;
     448  fTower->Ehad = hcalEnergy;
     449
    395450  fTower->Edges[0] = fTowerEdges[0];
    396451  fTower->Edges[1] = fTowerEdges[1];
     
    400455
    401456  // fill calorimeter towers
    402   if(energy > 0.0) fTowerOutputArray->Add(fTower);
     457  if(energy > 0.0)
     458  {
     459    if(fTowerPhotonHits > 0 && fTowerTrackHits == 0)
     460    {
     461      fPhotonOutputArray->Add(fTower);
     462    }
     463   
     464    fTowerOutputArray->Add(fTower);
     465  }
     466
     467  // fill energy flow candidates
     468
     469  // save all the tracks as energy flow tracks
     470  fItTowerTrackArray->Reset();
     471  while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
     472  {
     473    fEFlowTrackOutputArray->Add(track);
     474  }
     475
     476  ecalEnergy -= fTrackECalEnergy;
     477  if(ecalEnergy < 0.0) ecalEnergy = 0.0;
     478
     479  hcalEnergy -= fTrackHCalEnergy;
     480  if(hcalEnergy < 0.0) hcalEnergy = 0.0;
     481
     482  energy = ecalEnergy + hcalEnergy;
    403483
    404484 
    405   // fill energy flow candidates
    406   energy -= fTrackEnergy;
    407   if(energy < 0.0) energy = 0.0;
    408    
    409   // save energy excess as an energy flow tower
    410   if(energy > 0.0)
     485  // save ECAL and/or HCAL energy excess as an energy flow tower
     486  if(ecalEnergy > 0.0)
    411487  {
    412488    // create new photon tower
    413489    tower = static_cast<Candidate*>(fTower->Clone());
    414     pt = energy / TMath::CosH(eta);
    415 
    416     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
    417     fEFlowTowerOutputArray->Add(tower);
    418   }
     490
     491    pt = ecalEnergy / TMath::CosH(eta);
     492
     493    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy);
     494    tower->Eem = ecalEnergy;
     495    tower->Ehad = 0;
     496
     497    fEFlowPhotonOutputArray->Add(tower);
     498  }
     499
     500  if(hcalEnergy > 0.0)
     501  {
     502    // create new neutral hadron tower
     503    tower = static_cast<Candidate*>(fTower->Clone());
     504
     505    pt = hcalEnergy / TMath::CosH(eta);
     506
     507    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy);
     508    tower->Eem = 0;
     509    tower->Ehad = hcalEnergy;
     510
     511    fEFlowNeutralHadronOutputArray->Add(tower);
     512  }
     513
     514
     515
    419516
    420517}
Note: See TracChangeset for help on using the changeset viewer.