#include <hcalCalib.h>

Public Member Functions
virtual void	Begin (TTree *tree)
void	GetCoefFromMtrxInvOfAve ()
virtual Int_t	GetEntry (Long64_t entry, Int_t getall=0)
virtual TList *	GetOutputList () const
	hcalCalib (TTree *=0)
virtual void	Init (TTree *tree)
void	makeTextFile ()
virtual Bool_t	Notify ()
virtual Bool_t	Process (Long64_t entry)
Bool_t	ReadPhiSymCor ()
void	SetApplyPhiSymCorFlag (Bool_t b)
void	SetCalibAbsIEtaMax (Int_t i)
void	SetCalibAbsIEtaMin (Int_t i)
void	SetCalibMethod (TString s)
void	SetCalibType (TString s)
void	SetCaloGeometry (const CaloGeometry *g)
void	SetCombinePhiFlag (Bool_t b)
void	SetConeMaxDist (Float_t d)
void	SetHbClusterSize (Int_t i)
void	SetHeClusterSize (Int_t i)
void	SetHistoFileName (TString filename)
virtual void	SetInputList (TList *input)
void	SetMaxEOverP (Float_t e)
void	SetMaxEtThirdJet (Float_t et)
void	SetMaxProbeJetEmFrac (Float_t f)
void	SetMaxTagJetAbsEta (Float_t e)
void	SetMaxTagJetEmFrac (Float_t f)
void	SetMaxTargetE (Float_t e)
void	SetMaxTrkEmE (Float_t e)
void	SetMinCellE (Float_t e)
void	SetMinDPhiDiJets (Float_t dphi)
void	SetMinEOverP (Float_t e)
void	SetMinProbeJetAbsEta (Float_t e)
void	SetMinTagJetEt (Float_t e)
void	SetMinTargetE (Float_t e)
virtual void	SetObject (TObject *obj)
virtual void	SetOption (const char *option)
void	SetOutputCorCoefFileName (TString filename)
void	SetPhiSymCorFileName (TString filename)
void	SetSumDepthsFlag (Bool_t b)
void	SetSumSmallDepthsFlag (Bool_t b)
void	SetUseConeClustering (Bool_t b)
virtual void	Terminate ()
virtual Int_t	Version () const
virtual	~hcalCalib ()
Public Attributes
Bool_t	APPLY_PHI_SYM_COR_FLAG
TBranch *	b_cells
TBranch *	b_emEnergy
TBranch *	b_etVetoJet
TBranch *	b_eventNumber
TBranch *	b_iEtaHit
TBranch *	b_iPhiHit
TBranch *	b_probeJetEmFrac
TBranch *	b_probeJetP4
TBranch *	b_runNumber
TBranch *	b_tagJetEmFrac
TBranch *	b_tagJetP4
TBranch *	b_targetE
TBranch *	b_xTrkEcal
TBranch *	b_xTrkHcal
TBranch *	b_yTrkEcal
TBranch *	b_yTrkHcal
TBranch *	b_zTrkEcal
TBranch *	b_zTrkHcal
Int_t	CALIB_ABS_IETA_MAX
Int_t	CALIB_ABS_IETA_MIN
TString	CALIB_METHOD
TString	CALIB_TYPE
std::vector< std::vector < Float_t > >	cellEnergies
std::vector< std::vector < UInt_t > >	cellIds
TClonesArray *	cells
Bool_t	COMBINE_PHI
Float_t	emEnergy
Float_t	etVetoJet
UInt_t	eventNumber
	pointer to the analyzed TTree or TChain
TTree *	fChain
Int_t	HB_CLUSTER_SIZE
Int_t	HE_CLUSTER_SIZE
TString	HISTO_FILENAME
std::map< Int_t, Float_t >	iEtaCoefMap
Int_t	iEtaHit
UInt_t	iPhiHit
Float_t	MAX_CONE_DIST
Float_t	MAX_EOVERP
Float_t	MAX_ET_THIRD_JET
Float_t	MAX_PROBEJET_EMFRAC
Float_t	MAX_TAGJET_ABSETA
Float_t	MAX_TAGJET_EMFRAC
Float_t	MAX_TARGET_E
Float_t	MAX_TRK_EME
Float_t	MIN_CELL_E
Float_t	MIN_DPHI_DIJETS
Float_t	MIN_EOVERP
Float_t	MIN_PROBEJET_ABSETA
Float_t	MIN_TAGJET_ET
Float_t	MIN_TARGET_E
TString	OUTPUT_COR_COEF_FILENAME
TString	PHI_SYM_COR_FILENAME
std::map< UInt_t, Float_t >	phiSymCor
Float_t	probeJetEmFrac
TLorentzVector *	probeJetP4
std::vector< std::pair< Int_t, UInt_t > >	refIEtaIPhi
UInt_t	runNumber
std::map< UInt_t, Float_t >	solution
Bool_t	SUM_DEPTHS
Bool_t	SUM_SMALL_DEPTHS
Float_t	tagJetEmFrac
TLorentzVector *	tagJetP4
Float_t	targetE
std::vector< Float_t >	targetEnergies
const CaloGeometry *	theCaloGeometry
Bool_t	USE_CONE_CLUSTERING
Float_t	xTrkEcal
Float_t	xTrkHcal
Float_t	yTrkEcal
Float_t	yTrkHcal
Float_t	zTrkEcal
Float_t	zTrkHcal

Detailed Description

Definition at line 39 of file hcalCalib.h.

Constructor & Destructor Documentation

hcalCalib::hcalCalib ( TTree * = 0 ) [inline]

Definition at line 89 of file hcalCalib.h.

{ }

virtual hcalCalib::~hcalCalib ( ) [inline, virtual]

Definition at line 90 of file hcalCalib.h.

{ }

Member Function Documentation

void hcalCalib::Begin ( TTree * tree ) [virtual]

Definition at line 70 of file hcalCalib.cc.

References gather_cfg::cout, cmsRelvalreport::exit, h1_allTrkP, h1_corResp, h1_corRespBarrel, h1_corRespEndcap, h1_corRespIEta, h1_numEventsTwrIEta, h1_rawResp, h1_rawRespBarrel, h1_rawRespEndcap, h1_rawSumE, h1_selTrkP_iEta10, h1_trkP, h2_dHitRefBarrel, h2_dHitRefEndcap, histoFile, i, and nEvents.

                               {
  TString option = GetOption();

  nEvents = 0;

  if (APPLY_PHI_SYM_COR_FLAG && !ReadPhiSymCor()) {
    cout << "\nERROR: Failed to read the phi symmetry corrections." << endl;
    cout << "Check if the filename is correct. If the corrections are not needed, set the corresponding flag to \"false\"\n" << endl;

    cout << "\nThe program will be terminated\n" << endl;

    exit(1);

  }

    
  //  cellEnergies.reserve(1000000);
  //  cellIds.reserve(1000000);
  //  targetEnergies.reserve(1000000);
    
      histoFile = new TFile(HISTO_FILENAME.Data(), "RECREATE");


  h1_trkP    = new TH1F("h1_trkP", "Track momenta; p_{trk} (GeV); Number of tracks", 100, 0, 200);
  h1_allTrkP = new TH1F("h1_allTrkP", "Track momenta - all tracks; p_{trk} (GeV); Number of tracks", 100, 0, 200);

  h1_selTrkP_iEta10 = new TH1F("h1_selTrkP_iEta10", "Track momenta - tracks with |iEta|<10; p_{trk} (GeV); Number of tracks", 100, 0, 200); 



  if (CALIB_TYPE=="ISO_TRACK") 
    h1_rawSumE = new TH1F("h1_rawSumE", "Cluster Energy; E_{cl} (GeV); Number of tracks", 100, 0, 200);
  else 
    h1_rawSumE = new TH1F("h1_rawSumE", "Cluster Energy; E_{cl} (GeV); Number of tracks", 1000, 0, 2000);   



  h1_rawResp = new TH1F("h1_rawResp", "Uncorrected response: |iEta|<24;  E_{had}/p; Number of tracks", 300, 0, 3);
  h1_corResp = new TH1F("h1_corResp", "Corrected response: |iEta|<24; E_{had}/p; Number of tracks", 300, 0, 3);
    
  h1_rawRespBarrel = new TH1F("h1_rawRespBarrel", "Uncorrected response: |iEta|<15;  E_{had}/p; Number of tracks", 300, 0, 3);
  h1_corRespBarrel = new TH1F("h1_corRespBarrel", "Corrected response: |iEta|<15; E_{had}/p; Number of tracks", 300, 0, 3);
    
  h1_rawRespEndcap = new TH1F("h1_rawRespEndcap", "Uncorrected response:  17<|iEta|<24;  E_{had}/p; Number of tracks", 300, 0, 3);
  h1_corRespEndcap = new TH1F("h1_corRespEndcap", "Corrected response: 17<|iEta|<24; E_{had}/p; Number of tracks", 300, 0, 3);
  
  h1_numEventsTwrIEta = new TH1F("h1_numEventsTwrIEta", "h1_numEventsTwrIEta", 80, -40, 40);

  h2_dHitRefBarrel     = new TH2F("h2_dHitRefBarrel","{#Delta}i{#phi} vs {#Delta}i{#eta} of hit and most energetic tower(|i{#eta}|<16);{#Delta}i{#eta}; {#Delta}i{#phi}", 10, -5, 5, 10, -5, 5);
  h2_dHitRefEndcap     = new TH2F("h2_dHitRefEndcap","{#Delta}i{#phi} vs {#Delta}i{#eta} of hit and most energetic tower (16<|i{#eta}|<25) ;{#Delta}i{#eta}; {#Delta}i{#phi}", 10, -5, 5, 10, -5, 5);
      


  TString histoName = "isoTrack_";
 
  for (Int_t i=0; i<48; ++i) {
    Long_t iEta; 
    if   (i<24) iEta = i-24;
    else iEta = i-23;
    TString hn = histoName + iEta;
    h1_corRespIEta[i] = new TH1F(hn, hn, 300, 0, 3.0);
  }


}  // end of Begin()

void hcalCalib::GetCoefFromMtrxInvOfAve ( )

Definition at line 522 of file hcalCalib.cc.

References funct::A, abs, b, begin, end, spr::find(), i, HcalDetId::ieta(), j, lumiNorm::norm, findQualityFiles::size, and x.

                                        {

    // these maps are keyed by iEta        
    map<Int_t, Float_t> aveTargetE;     
    map<Int_t, Int_t>   nEntries;     // count hits
        
    //  iEtaRef  iEtaCell, energy 
    map<Int_t, map<Int_t, Float_t> > aveHitE; // add energies in the loop, normalize after that
 
    for (unsigned int i=0; i<cellEnergies.size(); ++i) {
      Int_t iEtaRef = refIEtaIPhi[i].first;
      aveTargetE[iEtaRef] += targetEnergies[i];
      nEntries[iEtaRef]++;

      // for hybrid method: matrix inv of averages preceeded by L3
      if (CALIB_METHOD=="L3_AND_MTRX_INV") {
        for (unsigned int j=0; j<(cellEnergies[i]).size(); ++j) {
          aveHitE[iEtaRef][HcalDetId(cellIds[i][j]).ieta()] += (solution[cellIds[i][j]] * cellEnergies[i][j]);
        }          
      }
      else if (CALIB_METHOD=="MATRIX_INV_OF_ETA_AVE") {
        for (unsigned int j=0; j<(cellEnergies[i]).size(); ++j) {
          aveHitE[iEtaRef][HcalDetId(cellIds[i][j]).ieta()] += cellEnergies[i][j];
        }          
      }

    } // end of loop of entries

           
    // scale by number of entries to get the averages
    Float_t norm = 1.0;
    for (map<Int_t, Float_t>::iterator m_it = aveTargetE.begin(); m_it!=aveTargetE.end(); ++m_it){            
      Int_t iEta = m_it->first;
      norm = (nEntries[iEta] > 0)? 1.0/(nEntries[iEta]) : 1.0;
      aveTargetE[iEta] *= norm;

      map<Int_t, Float_t>::iterator n_it = (aveHitE[iEta]).begin();

      Float_t sumRawE = 0;
      for (; n_it!=(aveHitE[iEta]).end(); ++n_it) {
        (n_it->second) *= norm;
        sumRawE += (n_it->second);
      }
        
    }  // end of scaling by number of entries
  
    Int_t ONE_SIDE_IETA_RANGE = CALIB_ABS_IETA_MAX - CALIB_ABS_IETA_MIN +1;

    // conversion from iEta to index for the linear system
    // contains elements only in the valid range for *matrix inversion*
    vector<Int_t> iEtaList;

    for (Int_t i=-CALIB_ABS_IETA_MAX; i<=CALIB_ABS_IETA_MAX; ++i) {
      if (abs(i)<CALIB_ABS_IETA_MIN) continue;
      iEtaList.push_back(i); 
    }   

    TMatrixD A(2*ONE_SIDE_IETA_RANGE, 2*ONE_SIDE_IETA_RANGE);
    TMatrixD b(2*ONE_SIDE_IETA_RANGE, 1);
    TMatrixD x(2*ONE_SIDE_IETA_RANGE, 1);
        
    for (Int_t i=0; i<2*ONE_SIDE_IETA_RANGE; ++i) {
      for (Int_t j=0; j<2*ONE_SIDE_IETA_RANGE; ++j) {
        A(i, j) = 0.0;
      }
    }

    for (UInt_t i=0; i<iEtaList.size(); ++i) {     
      b(i, 0) = aveTargetE[iEtaList[i]];        

      map<Int_t, Float_t>::iterator n_it = aveHitE[iEtaList[i]].begin();
      for (; n_it!=aveHitE[iEtaList[i]].end(); ++n_it){

        if (fabs(n_it->first)>CALIB_ABS_IETA_MAX || fabs(n_it->first)<CALIB_ABS_IETA_MIN) continue;
        Int_t j =  Int_t(find(iEtaList.begin(),iEtaList.end(), n_it->first) - iEtaList.begin());        
        A(i, j) = n_it->second;

      }

    }   

    TMatrixD coef = b;
    TDecompQRH qrh(A);
    Bool_t hasSolution = qrh.MultiSolve(coef);
   
    if (hasSolution && (CALIB_METHOD=="L3_AND_MTRX_INV" || CALIB_METHOD=="MATRIX_INV_OF_ETA_AVE")) {
      for (UInt_t i=0; i<iEtaList.size(); ++i) {
        iEtaCoefMap[iEtaList[i]] = coef(i, 0);
      }
    }
}

virtual Int_t hcalCalib::GetEntry	(	Long64_t	entry,
		Int_t	getall = `0`
	)		`[inline, virtual]`

Definition at line 97 of file hcalCalib.h.

References fChain.

{ return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; }

virtual TList* hcalCalib::GetOutputList ( ) const [inline, virtual]

Definition at line 101 of file hcalCalib.h.

{ return fOutput; }

virtual void hcalCalib::Init ( TTree * tree ) [virtual]

void hcalCalib::makeTextFile ( )

Definition at line 673 of file hcalCalib.cc.

References abs, eta(), HcalEndcap, HcalOuter, phi, and HcalDetId::validDetId().

                             {

  //{ HcalEmpty=0, HcalBarrel=1, HcalEndcap=2, HcalOuter=3, HcalForward=4, HcalTriggerTower=5, HcalOther=7 };
  
  TString sdName[8] = {"EMPTY", "HB", "HE", "HO", "HF", "TRITWR", "UNDEF", "OTHER"}; 

  FILE *constFile = fopen(OUTPUT_COR_COEF_FILENAME.Data(), "w+"); 

  // header of the constants file
  fprintf(constFile, "%1s%16s%16s%16s%16s%9s%11s\n", 
          "#", "eta", "phi", "depth", "det", "value", "DetId");

  // Order loops to produce sequence of constants as for phi symmetry 
  
  for(Int_t sd=1; sd<=4; sd++) {

    for(Int_t e=1; e<=50; e++) {
      Int_t eta;

      for(Int_t side=0; side<2; side++) {
        eta = (side==0)? -e : e; //ta *= (-1);

        for(Int_t phi=1; phi<=72; phi++) {

          for(Int_t d=1; d<5; d++) {

            if (!HcalDetId::validDetId(HcalSubdetector(sd), eta, phi, d)) continue;
            HcalDetId id(HcalSubdetector(sd), eta, phi, d);
            Float_t corrFactor = 1.0;
            

            if (abs(eta)>=CALIB_ABS_IETA_MIN && abs(eta)<=CALIB_ABS_IETA_MAX && HcalSubdetector(sd)!=HcalOuter) {
            //      if (abs(eta)>=CALIB_ABS_IETA_MIN && abs(eta)<=22 && HcalSubdetector(sd)!=HcalOuter) {


              // need some care when depths were summed for iEta=16 =>
              // the coeficients are saved in depth 1 of HB: affects
              Int_t subdetInd = sd;
              
              if (abs(eta)==16 && HcalSubdetector(sd)==HcalEndcap && SUM_DEPTHS) {
                subdetInd = 1;
              }

              if (CALIB_METHOD=="L3" || CALIB_METHOD=="L3_AND_MTRX_INV") {


                if (SUM_DEPTHS && COMBINE_PHI) corrFactor = solution[HcalDetId(HcalSubdetector(subdetInd), eta, 1, 1)];
                else if (SUM_DEPTHS) corrFactor = solution[HcalDetId(HcalSubdetector(subdetInd), eta, phi, 1)];
                else if (COMBINE_PHI) corrFactor = solution[HcalDetId(HcalSubdetector(sd), eta, 1, d)];
                else corrFactor = solution[HcalDetId(HcalSubdetector(sd), eta, phi, d)];


                // Remark: a new case was added (sumSmallDepths) where the first two depths in towers 15,16 
                // are summed and stored in  depth 1.
                // For now we create the correction coef for depth 2 (set equal to depth 1)
                // after the call to the L3 minimizer so that this case is also handled without modifying the
                // logic above. Probably it is better to move it here? 
                

              }// L3

              else if (CALIB_METHOD=="MATRIX_INV_OF_ETA_AVE") {
                corrFactor = iEtaCoefMap[eta];
              }

              else if (CALIB_METHOD=="ISO_TRK_PHI_SYM") {
                corrFactor = solution[HcalDetId(HcalSubdetector(sd), eta, phi, d)];
              }


            } // if within the calibration range
            
            fprintf(constFile, "%17i%16i%16i%16s%9.5f%11X\n", 
                    eta, phi, d, sdName[sd].Data(), corrFactor, id.rawId());

          }
        }
      }
    }
  }
   
  return;
}

virtual Bool_t hcalCalib::Notify ( ) [virtual]

Bool_t hcalCalib::Process ( Long64_t entry ) [virtual]

Definition at line 143 of file hcalCalib.cc.

References abs, diJetCalib::combinePhi, gather_cfg::cout, TCell::e(), filterCells3x3(), filterCells5x5(), filterCellsInCone(), getIEtaIPhiForHighestE(), h1_allTrkP, h1_rawSumE, h1_selTrkP_iEta10, h2_dHitRefBarrel, h2_dHitRefEndcap, HcalBarrel, HcalEndcap, HcalForward, HcalOuter, i, TCell::id(), M_PI, nEvents, edm::second(), TCell::SetE(), diJetCalib::sumDepths, and diJetCalib::sumSmallDepths.

                                        {
 
  //  fChain->GetTree()->GetEntry(entry);
  GetEntry(entry);


  set<UInt_t>  uniqueIds; // for testing: check if there are duplicate cells   (AA)

 
  Bool_t acceptEvent = kTRUE;
   
  ++nEvents;

  if (!(nEvents%100000)) cout << "event: " << nEvents << endl;


  h1_allTrkP->Fill(targetE);
   
  if (targetE < MIN_TARGET_E || targetE > MAX_TARGET_E) return kFALSE;;

 
  // make local copy as the cells may be modified  due to phi/depth sum, phi corrections etc  
  vector<TCell> selectCells;


  if (cells->GetSize()==0) return kFALSE;

  if (CALIB_TYPE=="DI_JET" && probeJetEmFrac > 0.999) return kTRUE;
  

  for (Int_t i=0; i<cells->GetSize(); ++i) {
    TCell* thisCell =  (TCell*) cells->At(i);

    if (HcalDetId(thisCell->id()).subdet()== HcalOuter)  continue;           // reject HO, make a switch!

    if (HcalDetId(thisCell->id()).subdet() != HcalBarrel &&
        HcalDetId(thisCell->id()).subdet() != HcalEndcap &&
        HcalDetId(thisCell->id()).subdet() != HcalForward) {
      
      cout << "Unknown or wrong hcal subdetector: " << HcalDetId(thisCell->id()).subdet() << endl;

    }


    // Apply phi symmetry corrections if the flag is set
    if (APPLY_PHI_SYM_COR_FLAG) thisCell->SetE(phiSymCor[thisCell->id()] * thisCell->e());

    if (thisCell->e() > MIN_CELL_E) selectCells.push_back(*thisCell);
  }


  if (selectCells.size()==0) {
    cout << "NO CELLS ABOVE THRESHOLD FOUND FOR TARGET!!!" << endl;
  }


  if (SUM_DEPTHS)  sumDepths(selectCells);
  else if (SUM_SMALL_DEPTHS) sumSmallDepths(selectCells); // depth 1,2 in twrs 15,16



  // most energetic tower (IsoTracks) or centroid of probe jet (DiJets)
  pair<Int_t, UInt_t> refPos;

  Int_t dEtaHitRef = 999;
  Int_t dPhiHitRef = 999;

  if (CALIB_TYPE=="ISO_TRACK") { 

    Int_t  iEtaMaxE;  // filled by reference in getIEtaIPhiForHighestE
    UInt_t iPhiMaxE;  // 
      
    getIEtaIPhiForHighestE(selectCells, iEtaMaxE, iPhiMaxE);   
 
    dEtaHitRef = iEtaMaxE - iEtaHit;
    dPhiHitRef = iPhiMaxE - iPhiHit;

    if (dPhiHitRef < -36) dPhiHitRef += 72;
    if (dPhiHitRef >  36) dPhiHitRef -= 72;

    if  (iEtaHit*iEtaMaxE < 0) {
      if (dEtaHitRef<0) dEtaHitRef += 1;
      if (dEtaHitRef>0) dEtaHitRef -= 1;     
    }

      
    if (abs(iEtaHit)<16) h2_dHitRefBarrel->Fill(dEtaHitRef, dPhiHitRef);
    if (abs(iEtaHit)>16 && abs(iEtaHit)<25) h2_dHitRefEndcap->Fill(dEtaHitRef, dPhiHitRef);

    // --------------------------------------------------
    // Choice of cluster definition 
    //
    // fixed size NxN clusters as specified in to config file 
    if (!USE_CONE_CLUSTERING) {
      if (abs(iEtaMaxE)<16 && HB_CLUSTER_SIZE==3) filterCells3x3(selectCells, iEtaMaxE, iPhiMaxE);
      if (abs(iEtaMaxE)>15 && HE_CLUSTER_SIZE==3) filterCells3x3(selectCells, iEtaMaxE, iPhiMaxE);
      
      if (abs(iEtaMaxE)<16 && HB_CLUSTER_SIZE==5) filterCells5x5(selectCells, iEtaMaxE, iPhiMaxE);
      if (abs(iEtaMaxE)>15 && HE_CLUSTER_SIZE==5) filterCells5x5(selectCells, iEtaMaxE, iPhiMaxE);    
    }      
    else {
      //  calculate distance at hcal surface
      const GlobalPoint hitPositionHcal(xTrkHcal, yTrkHcal, zTrkHcal); 
      filterCellsInCone(selectCells, hitPositionHcal, MAX_CONE_DIST, theCaloGeometry);
    }



    refPos.first  = iEtaMaxE;
    refPos.second = iPhiMaxE;

  }
  else if (CALIB_TYPE=="DI_JET") {          // Apply selection cuts on DiJet events here

    if (etVetoJet>MAX_ET_THIRD_JET) acceptEvent = kFALSE;
    
    Float_t jetsDPhi = probeJetP4->DeltaPhi(*tagJetP4); 
    if (fabs(jetsDPhi * 180.0/M_PI) < MIN_DPHI_DIJETS) acceptEvent = kFALSE;

    if (probeJetEmFrac>MAX_PROBEJET_EMFRAC)          acceptEvent = kFALSE;
    if (fabs(probeJetP4->Eta())<MIN_PROBEJET_ABSETA) acceptEvent = kFALSE;
    if (fabs(tagJetP4->Eta())>MAX_TAGJET_ABSETA)     acceptEvent = kFALSE;
    if (fabs(tagJetP4->Et())< MIN_TAGJET_ET)         acceptEvent = kFALSE;

    if (acceptEvent) {

      Int_t  iEtaMaxE;  // filled by reference in getIEtaIPhiForHighestE
      UInt_t iPhiMaxE;  // 
      
      getIEtaIPhiForHighestE(selectCells, iEtaMaxE, iPhiMaxE);   
      
      // The ref position for the jet is not used in the minimization at this time.
      // It will be needed if we attempt to do matrix inversion: then the question is
      // which value is better suited: the centroid of the jet or the hottest tower...

      //    refPos.first  = iEtaHit;
      //    refPos.second = iPhiHit;
      
      refPos.first  = iEtaMaxE;
      refPos.second = iPhiMaxE;  
      
      if (abs(iEtaMaxE)>40) acceptEvent = kFALSE;             // for testing :  set as parameter (AA)

    }

  }


  if (COMBINE_PHI) combinePhi(selectCells);
    
  // fill the containers for the minimization prcedures
  vector<Float_t> energies;
  vector<UInt_t>  ids;
    
  for (vector<TCell>::iterator i_it = selectCells.begin(); i_it!=selectCells.end(); ++i_it) {

    // for testing : fill only unique id's

    if (uniqueIds.insert(i_it->id()).second) {
      energies.push_back(i_it->e());
      ids.push_back(i_it->id());
    }

  }

  if (CALIB_TYPE=="ISO_TRACK") {
    if (accumulate(energies.begin(), energies.end(), 0.0) / targetE < MIN_EOVERP) acceptEvent=kFALSE;
    if (accumulate(energies.begin(), energies.end(), 0.0) / targetE > MAX_EOVERP) acceptEvent=kFALSE;

    if (emEnergy > MAX_TRK_EME) acceptEvent=kFALSE;

    if (abs(dEtaHitRef)>1 || abs(dPhiHitRef)>1) acceptEvent=kFALSE;

    // Have to check if for |iEta|>20 (and neighboring region) the dPhiHitRef
    // should be relaxed to 2. The neighboring towers have dPhi=2...


  }

    
  h1_rawSumE->Fill(accumulate(energies.begin(), energies.end(), 0.0));
    
     
  // here we fill the information for the minimization procedure
  if (acceptEvent) { 
    cellEnergies.push_back(energies);
    cellIds.push_back(ids);
    targetEnergies.push_back(targetE);
    refIEtaIPhi.push_back(refPos);

    if (abs(refPos.first)<=10) 
      h1_selTrkP_iEta10->Fill(targetE);


  }



    
  // Clean up
  energies.clear();
  ids.clear();
  selectCells.clear();
    
  return kTRUE;
}

Bool_t hcalCalib::ReadPhiSymCor ( )

Definition at line 615 of file hcalCalib.cc.

References gather_cfg::cout, HcalBarrel, HcalEndcap, HcalForward, HcalOuter, geometryCSVtoXML::line, HcalDetId::validDetId(), and relativeConstraints::value.

                                {

  ifstream phiSymFile(PHI_SYM_COR_FILENAME.Data());

  if (!phiSymFile) {
    cout << "\nERROR: Can not find file with phi symmetry constants \"" <<  PHI_SYM_COR_FILENAME.Data() << "\"" << endl;
    return kFALSE;
  }

  // assumes the format used in CSA08, first line is a comment

  Int_t   iEta;
  UInt_t  iPhi;
  UInt_t  depth;
  TString sdName;
  UInt_t  detId;

  Float_t value;
  HcalSubdetector sd;

  std::string line;

  while (getline(phiSymFile, line)) {

    if(!line.size() || line[0]=='#') continue;

    std::istringstream linestream(line);
    linestream >> iEta >> iPhi >> depth >> sdName >> value >> hex >> detId;

    if (sdName=="HB") sd = HcalBarrel;
    else if (sdName=="HE") sd = HcalEndcap;
    else if (sdName=="HO") sd = HcalOuter;
    else if (sdName=="HF") sd = HcalForward;
    else {
      cout << "\nInvalid detector name in phi symmetry constants file: " << sdName.Data() << endl;
      cout << "Check file and rerun!\n" << endl;
      return kFALSE;
    }

    // check if the data is consistent    
  
    if (HcalDetId(sd, iEta, iPhi, depth) != HcalDetId(detId)) {
      cout << "\nInconsistent info in phi symmetry file: subdet, iEta, iPhi, depth do not match rawId!\n" << endl;
      return kFALSE;    
    }
    if (!HcalDetId::validDetId(sd, iEta, iPhi, depth)) {
      cout << "\nInvalid DetId from: iEta=" << iEta << " iPhi=" << iPhi << " depth=" << depth 
           << " subdet=" << sdName.Data() << " detId=" << detId << endl << endl; 
      return kFALSE;    
    }

    phiSymCor[HcalDetId(sd, iEta, iPhi, depth)] = value;
  }

  return kTRUE;
}

void hcalCalib::SetApplyPhiSymCorFlag ( Bool_t b ) [inline]

Definition at line 171 of file hcalCalib.h.

References APPLY_PHI_SYM_COR_FLAG, and b.

Referenced by HcalCalibrator::endJob().

{ APPLY_PHI_SYM_COR_FLAG = b; }

void hcalCalib::SetCalibAbsIEtaMax ( Int_t i ) [inline]

Definition at line 167 of file hcalCalib.h.

References CALIB_ABS_IETA_MAX, and i.

Referenced by HcalCalibrator::endJob().

{ CALIB_ABS_IETA_MAX = i; }

void hcalCalib::SetCalibAbsIEtaMin ( Int_t i ) [inline]

Definition at line 168 of file hcalCalib.h.

References CALIB_ABS_IETA_MIN, and i.

Referenced by HcalCalibrator::endJob().

{ CALIB_ABS_IETA_MIN = i; }

void hcalCalib::SetCalibMethod ( TString s ) [inline]

Definition at line 160 of file hcalCalib.h.

References CALIB_METHOD, and asciidump::s.

Referenced by HcalCalibrator::endJob().

{ CALIB_METHOD = s; }

void hcalCalib::SetCalibType ( TString s ) [inline]

Definition at line 159 of file hcalCalib.h.

References CALIB_TYPE, and asciidump::s.

Referenced by HcalCalibrator::endJob().

{ CALIB_TYPE = s; }

void hcalCalib::SetCaloGeometry ( const CaloGeometry * g ) [inline]

Definition at line 182 of file hcalCalib.h.

References g, and theCaloGeometry.

Referenced by HcalCalibrator::endJob().

{ theCaloGeometry = g; }

void hcalCalib::SetCombinePhiFlag ( Bool_t b ) [inline]

Definition at line 154 of file hcalCalib.h.

References b, and COMBINE_PHI.

Referenced by HcalCalibrator::endJob().

{ COMBINE_PHI = b; }

void hcalCalib::SetConeMaxDist ( Float_t d ) [inline]

Definition at line 165 of file hcalCalib.h.

References MAX_CONE_DIST.

Referenced by HcalCalibrator::endJob().

{ MAX_CONE_DIST = d; }

void hcalCalib::SetHbClusterSize ( Int_t i ) [inline]

Definition at line 161 of file hcalCalib.h.

References HB_CLUSTER_SIZE, and i.

Referenced by HcalCalibrator::endJob().

{ HB_CLUSTER_SIZE = i; }

void hcalCalib::SetHeClusterSize ( Int_t i ) [inline]

Definition at line 162 of file hcalCalib.h.

References HE_CLUSTER_SIZE, and i.

Referenced by HcalCalibrator::endJob().

{ HE_CLUSTER_SIZE = i; }

void hcalCalib::SetHistoFileName ( TString filename ) [inline]

Definition at line 179 of file hcalCalib.h.

References lut2db_cfg::filename, and HISTO_FILENAME.

Referenced by HcalCalibrator::endJob().

{ HISTO_FILENAME = filename; }

virtual void hcalCalib::SetInputList ( TList * input ) [inline, virtual]

Definition at line 100 of file hcalCalib.h.

References collect_tpl::input.

{ fInput = input; }

void hcalCalib::SetMaxEOverP ( Float_t e ) [inline]

Definition at line 157 of file hcalCalib.h.

References MAX_EOVERP.

Referenced by HcalCalibrator::endJob().

{ MAX_EOVERP = e; }

void hcalCalib::SetMaxEtThirdJet ( Float_t et ) [inline]

Definition at line 169 of file hcalCalib.h.

References MAX_ET_THIRD_JET.

Referenced by HcalCalibrator::endJob().

{ MAX_ET_THIRD_JET = et; }

void hcalCalib::SetMaxProbeJetEmFrac ( Float_t f ) [inline]

Definition at line 173 of file hcalCalib.h.

References f, and MAX_PROBEJET_EMFRAC.

Referenced by HcalCalibrator::endJob().

{ MAX_PROBEJET_EMFRAC = f; }

void hcalCalib::SetMaxTagJetAbsEta ( Float_t e ) [inline]

Definition at line 175 of file hcalCalib.h.

References MAX_TAGJET_ABSETA.

Referenced by HcalCalibrator::endJob().

{ MAX_TAGJET_ABSETA = e; }

void hcalCalib::SetMaxTagJetEmFrac ( Float_t f ) [inline]

Definition at line 174 of file hcalCalib.h.

References f, and MAX_TAGJET_EMFRAC.

Referenced by HcalCalibrator::endJob().

{ MAX_TAGJET_EMFRAC = f; }

void hcalCalib::SetMaxTargetE ( Float_t e ) [inline]

Definition at line 151 of file hcalCalib.h.

References MAX_TARGET_E.

Referenced by HcalCalibrator::endJob().

{ MAX_TARGET_E = e; }

void hcalCalib::SetMaxTrkEmE ( Float_t e ) [inline]

Definition at line 158 of file hcalCalib.h.

References MAX_TRK_EME.

Referenced by HcalCalibrator::endJob().

{ MAX_TRK_EME = e; }

void hcalCalib::SetMinCellE ( Float_t e ) [inline]

Definition at line 155 of file hcalCalib.h.

References MIN_CELL_E.

Referenced by HcalCalibrator::endJob().

{ MIN_CELL_E = e; }

void hcalCalib::SetMinDPhiDiJets ( Float_t dphi ) [inline]

Definition at line 170 of file hcalCalib.h.

References MIN_DPHI_DIJETS.

Referenced by HcalCalibrator::endJob().

{ MIN_DPHI_DIJETS = dphi; }

void hcalCalib::SetMinEOverP ( Float_t e ) [inline]

Definition at line 156 of file hcalCalib.h.

References MIN_EOVERP.

Referenced by HcalCalibrator::endJob().

{ MIN_EOVERP = e; }

void hcalCalib::SetMinProbeJetAbsEta ( Float_t e ) [inline]

Definition at line 177 of file hcalCalib.h.

References MIN_PROBEJET_ABSETA.

Referenced by HcalCalibrator::endJob().

{ MIN_PROBEJET_ABSETA = e; }

void hcalCalib::SetMinTagJetEt ( Float_t e ) [inline]

Definition at line 176 of file hcalCalib.h.

References MIN_TAGJET_ET.

Referenced by HcalCalibrator::endJob().

{ MIN_TAGJET_ET = e; }

void hcalCalib::SetMinTargetE ( Float_t e ) [inline]

Definition at line 150 of file hcalCalib.h.

References MIN_TARGET_E.

Referenced by HcalCalibrator::endJob().

{ MIN_TARGET_E = e; }

virtual void hcalCalib::SetObject ( TObject * obj ) [inline, virtual]

Definition at line 99 of file hcalCalib.h.

References VarParsing::obj.

{ fObject = obj; }

virtual void hcalCalib::SetOption ( const char * option ) [inline, virtual]

Definition at line 98 of file hcalCalib.h.

{ fOption = option; }

void hcalCalib::SetOutputCorCoefFileName ( TString filename ) [inline]

Definition at line 178 of file hcalCalib.h.

References lut2db_cfg::filename, and OUTPUT_COR_COEF_FILENAME.

Referenced by HcalCalibrator::endJob().

{ OUTPUT_COR_COEF_FILENAME = filename; }

void hcalCalib::SetPhiSymCorFileName ( TString filename ) [inline]

Definition at line 172 of file hcalCalib.h.

References lut2db_cfg::filename, and PHI_SYM_COR_FILENAME.

Referenced by HcalCalibrator::endJob().

{ PHI_SYM_COR_FILENAME = filename; }

void hcalCalib::SetSumDepthsFlag ( Bool_t b ) [inline]

Definition at line 152 of file hcalCalib.h.

References b, and SUM_DEPTHS.

Referenced by HcalCalibrator::endJob().

{ SUM_DEPTHS = b; }

void hcalCalib::SetSumSmallDepthsFlag ( Bool_t b ) [inline]

Definition at line 153 of file hcalCalib.h.

References b, and SUM_SMALL_DEPTHS.

Referenced by HcalCalibrator::endJob().

{ SUM_SMALL_DEPTHS = b; }

void hcalCalib::SetUseConeClustering ( Bool_t b ) [inline]

Definition at line 164 of file hcalCalib.h.

References b, and USE_CONE_CLUSTERING.

Referenced by HcalCalibrator::endJob().

{ USE_CONE_CLUSTERING = b; }

void hcalCalib::Terminate ( ) [virtual]

Definition at line 352 of file hcalCalib.cc.

References abs, gather_cfg::cout, HcalDetId::depth(), h1_allTrkP, h1_corResp, h1_corRespBarrel, h1_corRespEndcap, h1_corRespIEta, h1_numEventsTwrIEta, h1_rawResp, h1_rawRespBarrel, h1_rawRespEndcap, h1_rawSumE, h1_selTrkP_iEta10, h1_trkP, h2_dHitRefBarrel, h2_dHitRefEndcap, HcalBarrel, histoFile, i, HcalDetId::ieta(), HcalDetId::ietaAbs(), HcalDetId::iphi(), MinL3AlgoUniv< IDdet >::iterate(), j, findQualityFiles::size, HcalDetId::subdet(), and CommonMethods::weight().

                          {

  cout << "\n\nFinished reading the events.\n";
  cout << "Number of input objects: " << cellIds.size()  << endl; 
  cout << "Performing minimization: depending on selected method can take some time...\n\n";
  
  for (vector<pair<Int_t, UInt_t> >::iterator it_rp=refIEtaIPhi.begin(); it_rp!=refIEtaIPhi.end(); ++it_rp ) {
    Float_t weight = (abs(it_rp->first)<21)? 1.0/72.0 : 1.0/36.0;
    h1_numEventsTwrIEta->Fill(it_rp->first, weight);
  }
 
  if (CALIB_METHOD=="MATRIX_INV_OF_ETA_AVE") {
    GetCoefFromMtrxInvOfAve();
   }    
  else if (CALIB_METHOD=="L3" || CALIB_METHOD=="L3_AND_MTRX_INV") {
    int eventWeight = 2;  // 2 is the default (try at some point 0,1,2,3)
    MinL3AlgoUniv<UInt_t>* thisL3Algo = new MinL3AlgoUniv<UInt_t>(eventWeight);
    int numIterations = 10;  // default 10
        
    solution = thisL3Algo->iterate(cellEnergies, cellIds, targetEnergies, numIterations);

    // in order to handle the case where sumDepths="false", but the flag to sum depths 1,2 in HB towers 15, 16
    // is set (sumSmallDepths) we create entries in "solution" to create equal correction here
    // for each encountered coef in depth one.

    if (!SUM_DEPTHS && SUM_SMALL_DEPTHS) {
      vector<UInt_t> idForSummedCells;

      for (map<UInt_t, Float_t>::iterator m_it = solution.begin(); m_it != solution.end(); ++m_it) {
        if (HcalDetId(m_it->first).ietaAbs()!=15 && HcalDetId(m_it->first).ietaAbs()!=16) continue;
        if (HcalDetId(m_it->first).subdet()!=HcalBarrel) continue;
        if (HcalDetId(m_it->first).depth()==1)
          idForSummedCells.push_back(HcalDetId(m_it->first));
      }

      for (vector<UInt_t>::iterator v_it=idForSummedCells.begin(); v_it!=idForSummedCells.end(); ++v_it) {
        UInt_t addCoefId = HcalDetId(HcalBarrel, HcalDetId(*v_it).ieta(), HcalDetId(*v_it).iphi(), 2);
        solution[addCoefId] = solution[*v_it];
      }
     
    }  // end of special treatment for "sumSmallDepths" mode


    if (CALIB_METHOD=="L3_AND_MTRX_INV") {
      GetCoefFromMtrxInvOfAve();      
      
      // loop over the solution from L3 and multiply by the additional factor from
      // the matrix inversion. Set coef outside of the valid calibration region =1.       
      for (map<UInt_t, Float_t>::iterator it_s=solution.begin(); it_s!=solution.end(); ++it_s) {
        Int_t iEtaSol = HcalDetId(it_s->first).ieta();
        if (abs(iEtaSol)<CALIB_ABS_IETA_MIN || abs(iEtaSol)> CALIB_ABS_IETA_MAX) it_s->second = 1.0;
        else  it_s->second *= iEtaCoefMap[iEtaSol]; 
      }    
      
    } // if CALIB_METHOD=="L3_AND_MTRX_INV"
    
  }  // end of L3 or L3 + mtrx inversion minimization

  // done getting the constants -> write the formatted file
  makeTextFile();
  
  // fill some histograms
  Float_t rawResp = 0;
  Float_t corResp = 0;
  Int_t maxIEta = 0;
  Int_t minIEta = 999;


  for (unsigned int i=0; i<cellEnergies.size(); ++i) {
    Int_t iEta;    
    for (unsigned int j=0; j<(cellEnergies[i]).size(); ++j) {
      iEta = HcalDetId(cellIds[i][j]).ieta();
      rawResp += (cellEnergies[i])[j];
      
      if (CALIB_METHOD=="L3_AND_MTRX_INV") {
        corResp += (solution[cellIds[i][j]] * cellEnergies[i][j]);
      }
      else if (CALIB_METHOD=="L3") {
        corResp += (solution[cellIds[i][j]] * (cellEnergies[i][j]));
      }
      else if (CALIB_METHOD=="MATRIX_INV_OF_ETA_AVE") {
        corResp += (iEtaCoefMap[iEta]* cellEnergies[i][j]);     
      }
      
      if (maxIEta < abs(iEta) ) maxIEta = abs(iEta);
      if (minIEta > abs(iEta) ) minIEta = abs(iEta); 
    }

    rawResp /= targetEnergies[i];
    corResp /= targetEnergies[i];


    // fill histograms based on iEta on ref point of the cluster (for now: hottest tower) 
    // expect range |iEta|<=24 (to do: add flexibility for arbitrary range)
 
    if (CALIB_TYPE=="ISO_TRACK") {
      Int_t ind = refIEtaIPhi[i].first;  
      (ind<0) ? (ind +=24) : (ind +=23);
      if (ind>=0 && ind<48) {
        h1_corRespIEta[ind]->Fill(corResp);  
      }
     
      // fill histograms for cases where all towers are in the barrel or endcap
      if (maxIEta<25) {
        h1_rawResp->Fill(rawResp);
        h1_corResp->Fill(corResp);
      }
      if (maxIEta<15) {
        h1_rawRespBarrel->Fill(rawResp);
        h1_corRespBarrel->Fill(corResp);
      }
      else if (maxIEta<25 && minIEta>16) {
        h1_rawRespEndcap->Fill(rawResp);
        h1_corRespEndcap->Fill(corResp);
      }
    }  // histograms for isotrack calibration


    else {
     
      // put jet plots here

    }



    rawResp = 0;
    corResp = 0;
    maxIEta = 0;
    minIEta = 999;
  }
   

  // save the histograms 
  h1_trkP->Write();
  h1_allTrkP->Write();

h1_selTrkP_iEta10->Write();

  h1_rawSumE->Write();
  h1_rawResp->Write();
  h1_corResp->Write();
  h1_rawRespBarrel->Write();
  h1_corRespBarrel->Write();
  h1_rawRespEndcap->Write();
  h1_corRespEndcap->Write();
  h1_numEventsTwrIEta->Write();
  h2_dHitRefBarrel->Write();
  h2_dHitRefEndcap->Write();
  for (Int_t i=0; i<48; ++i) {
    h1_corRespIEta[i]->Write();
  }





  histoFile->Write();
  histoFile->Close();

  cout << "\n Finished calibration.\n " << endl; 

    
}  // end of Terminate()