From 28f9dfe644439c35c228ddfa3bf9f33a875c377e Mon Sep 17 00:00:00 2001 From: John Chapman <jchapman@cern.ch> Date: Tue, 12 Sep 2017 15:01:24 +0200 Subject: [PATCH] Cleaned-up `CaloHitAna`. Reduced `CaloHitAna` memory usage. This is a manual sweep of the changes from !3978 into `master`. Added an option to enable/disable recording of g4hits. The reason for these changes is to fix the memory problems that were present in the previous version of the code. There were minor problems with the initialisation of new branches for the output tree, there were some vectors that weren't getting cleared after each loop, and there were a couple of objects that weren't being used anymore which took up too much memory. The initialisation of branches has been moved to the header file, lines have been added to the source file to clean the vectors after each loop and the unnecessary objects that accounted for most of the memory problems have been removed. --- .../tools/CaloHitAna.C | 1076 +++++++---------- .../tools/CaloHitAna.h | 232 ++-- 2 files changed, 571 insertions(+), 737 deletions(-) diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/CaloHitAna.C b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/CaloHitAna.C index f0928b8255f..061a6bbf6cf 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/CaloHitAna.C +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/CaloHitAna.C @@ -16,7 +16,6 @@ using namespace std; void CaloHitAna::Loop() { - //This Loop is supposed to read the output ntuple/tree //(vector of cells, vector of G4 hits, vector of FCS detailed hits, ...) //from the athena ISF_HitAnalysis algorithm (==ESD dumper) @@ -24,11 +23,11 @@ void CaloHitAna::Loop() //NB: E = Ehit * SamplingFraction for Tile or E = Ehit / SamplingFraction (LAr) conversion //is done here - if(fChain == 0) return; + if (fChain == 0) return; ProcInfo_t procinfo; Long64_t nentries = fChain->GetEntriesFast(); - if(m_max_nentries>=0 && m_max_nentries<nentries) nentries=m_max_nentries; + if (m_max_nentries >= 0 && m_max_nentries < nentries) nentries = m_max_nentries; std::map<Long64_t, FCS_cell> cells; //read all objects and collect them by identifier (Long64_t) std::map<Long64_t, std::vector<FCS_g4hit> > g4hits; @@ -43,618 +42,471 @@ void CaloHitAna::Loop() //From here: Loop over events: Long64_t nbytes = 0, nb = 0; std::cout << "before event loop" << std::endl; - for(Long64_t jentry=0; jentry<nentries;jentry++) + for (Long64_t jentry = 0; jentry < nentries; jentry++) + { + if (jentry % 100 == 0) + std::cout << "jentry " << jentry << endl; + + Long64_t ientry = LoadTree(jentry); + if (ientry < 0) break; + nb = fChain->GetEntry(jentry); nbytes += nb; + if (jentry % m_PrintOutFrequency == 0) + { + std::cout << "Processing entry " << jentry << std::endl; + gSystem->GetProcInfo(&procinfo); + std::cout << "Memory usage: " << procinfo.fMemResident << " " << procinfo.fMemVirtual << std::endl; + } + + //Copy truth + new_truthE->clear(); + new_truthPx->clear(); + new_truthPy->clear(); + new_truthPz->clear(); + new_truthPDG->clear(); + new_truthBarcode->clear(); + new_truthVtxBarcode->clear(); + //truthcollection->clear(); + m_newTTC_entrance_eta->clear(); + m_newTTC_entrance_phi->clear(); + m_newTTC_entrance_r->clear(); + m_newTTC_entrance_z->clear(); + m_newTTC_back_eta->clear(); + m_newTTC_back_phi->clear(); + m_newTTC_back_r->clear(); + m_newTTC_back_z->clear(); + m_newTTC_IDCaloBoundary_eta->clear(); + m_newTTC_IDCaloBoundary_phi->clear(); + m_newTTC_IDCaloBoundary_r->clear(); + m_newTTC_IDCaloBoundary_z->clear(); + m_newTTC_Angle3D->clear(); + m_newTTC_AngleEta->clear(); + + oneeventcells->m_vector.clear(); + + if (TruthE->size() != 1) + std::cout << "Strange! TruthE->size() is " << TruthE->size() << ", but it should be 1. Please check input." << std::endl; + + for (unsigned int truth_i = 0; truth_i < TruthE->size(); truth_i++) { - if(jentry%100==0) - std::cout<<"jentry "<<jentry<<endl; - - Long64_t ientry = LoadTree(jentry); - if (ientry < 0) break; - nb = fChain->GetEntry(jentry); nbytes += nb; - if (jentry % m_PrintOutFrequency == 0) - { - std::cout <<"Processing entry "<<jentry<<std::endl; - gSystem->GetProcInfo(&procinfo); - std::cout <<"Memory usage: "<<procinfo.fMemResident<<" "<<procinfo.fMemVirtual<<std::endl; - } - - //Copy truth - new_truthE->clear(); - new_truthPx->clear(); - new_truthPy->clear(); - new_truthPz->clear(); - new_truthPDG->clear(); - new_truthBarcode->clear(); - new_truthVtxBarcode->clear(); - //truthcollection->clear(); - m_newTTC_entrance_eta->clear(); - m_newTTC_entrance_phi->clear(); - m_newTTC_entrance_r->clear(); - m_newTTC_entrance_z->clear(); - m_newTTC_back_eta->clear(); - m_newTTC_back_phi->clear(); - m_newTTC_back_r->clear(); - m_newTTC_back_z->clear(); - m_newTTC_IDCaloBoundary_eta->clear(); - m_newTTC_IDCaloBoundary_phi->clear(); - m_newTTC_IDCaloBoundary_r->clear(); - m_newTTC_IDCaloBoundary_z->clear(); - m_newTTC_Angle3D->clear(); - m_newTTC_AngleEta->clear(); - - oneeventcells->m_vector.clear(); - - if(TruthE->size()!=1) - std::cout<<"Strange! TruthE->size() is "<<TruthE->size()<<", but it should be 1. Please check input."<<std::endl; - - for(unsigned int truth_i = 0; truth_i < TruthE->size(); truth_i++) - { - //std::cout <<truth_i<<" "<<TruthE->size()<<" "<<TruthPDG->at(truth_i)<<std::endl; - //std::cout <<"TTC2: "<<(*TTC_entrance_eta)[truth_i].size()<<std::endl;//this one has 24 layers ->ok - - new_truthE->push_back(TruthE->at(truth_i)); - new_truthPx->push_back(TruthPx->at(truth_i)); - new_truthPy->push_back(TruthPy->at(truth_i)); - new_truthPz->push_back(TruthPz->at(truth_i)); - new_truthPDG->push_back(TruthPDG->at(truth_i)); - new_truthBarcode->push_back(TruthBarcode->at(truth_i)); - //new_truthVtxBarcode->push_back(TruthVtxBarcode->at(truth_i)); - - m_newTTC_IDCaloBoundary_eta->push_back(newTTC_IDCaloBoundary_eta->at(truth_i)); - m_newTTC_IDCaloBoundary_phi->push_back(newTTC_IDCaloBoundary_phi->at(truth_i)); - m_newTTC_IDCaloBoundary_r ->push_back(newTTC_IDCaloBoundary_r->at(truth_i)); - m_newTTC_IDCaloBoundary_z ->push_back(newTTC_IDCaloBoundary_z->at(truth_i)); - m_newTTC_Angle3D ->push_back(newTTC_Angle3D->at(truth_i)); - m_newTTC_AngleEta ->push_back(newTTC_AngleEta->at(truth_i)); - - //create temporary vectors - vector<float> new_entrance_eta; - vector<float> new_entrance_phi; - vector<float> new_entrance_r; - vector<float> new_entrance_z; - vector<float> new_back_eta; - vector<float> new_back_phi; - vector<float> new_back_r; - vector<float> new_back_z; - - for(unsigned int s=0;s<24;s++) - { - new_entrance_eta.push_back((newTTC_entrance_eta->at(truth_i))[s]); - new_entrance_phi.push_back((newTTC_entrance_phi->at(truth_i))[s]); - new_entrance_r.push_back((newTTC_entrance_r->at(truth_i))[s]); - new_entrance_z.push_back((newTTC_entrance_z->at(truth_i))[s]); - new_back_eta.push_back((newTTC_back_eta->at(truth_i))[s]); - new_back_phi.push_back((newTTC_back_phi->at(truth_i))[s]); - new_back_r.push_back((newTTC_back_r->at(truth_i))[s]); - new_back_z.push_back((newTTC_back_z->at(truth_i))[s]); - } - - //push back temporary vectors - m_newTTC_entrance_eta->push_back(new_entrance_eta); - m_newTTC_entrance_phi->push_back(new_entrance_phi); - m_newTTC_entrance_r->push_back(new_entrance_r); - m_newTTC_entrance_z->push_back(new_entrance_z); - m_newTTC_back_eta->push_back(new_back_eta); - m_newTTC_back_phi->push_back(new_back_phi); - m_newTTC_back_r->push_back(new_back_r); - m_newTTC_back_z->push_back(new_back_z); - - /* - one_truth.SetPxPyPzE((*TruthPx)[truth_i], (*TruthPy)[truth_i], (*TruthPz)[truth_i], (*TruthE)[truth_i]); - one_truth.pdgid = TruthPDG->at(truth_i); - one_truth.barcode = TruthBarcode->at(truth_i); - one_truth.vtxbarcode = TruthVtxBarcode->at(truth_i); - truthcollection->push_back(one_truth); - */ - - } - - //std::cout<<"check. after truth block"<<std::endl; - - //Now matching between cells, G4hits and G4detailed hits - //sort cells by identifier: - //clear first the containers for this event: - cells.clear(); - g4hits.clear(); - hits.clear(); - //std::cout <<"Check Original size: Cells: "<<CellIdentifier->size()<<" G4Hits: "<<G4HitCellIdentifier->size()<<" FCSHits: "<<HitCellIdentifier->size()<<std::endl; - if(m_Debug > 1) std::cout <<"Reading cells..."; - - for (unsigned int cell_i = 0; cell_i <CellIdentifier->size(); cell_i++) - { - if (cells.find((*CellIdentifier)[cell_i]) == cells.end()) //doesn't exist - { - one_cell.cell_identifier = (*CellIdentifier)[cell_i]; - one_cell.sampling = (*CellSampling)[cell_i]; - one_cell.energy = (*CellE)[cell_i]; - one_cell.center_x = 0.0; //for now - one_cell.center_y = 0.0; - one_cell.center_z = 0.0; - cells.insert(std::pair<Long64_t, FCS_cell>(one_cell.cell_identifier, one_cell)); - } - else - { - //there shouldn't be a cell with the same identifier in this event - std::cout <<"ISF_HitAnalysis: Same cell???? ERROR"<<std::endl; - } - } - - if(m_Debug > 1) std::cout <<" Done"<<std::endl; - - if (m_Debug > 1) std::cout <<"Reading G4hits"; - - if (m_do_g4_hits){ - for (unsigned int g4hit_i = 0; g4hit_i <G4HitIdentifier->size(); g4hit_i++) - { - if ((*G4HitSampling)[g4hit_i] >=0 && (*G4HitSampling)[g4hit_i]<=25 && (*G4HitT)[g4hit_i]>m_TimingCut) - { - if (m_Debug > 1) std::cout <<"Ignoring G4hit, time too large: "<<g4hit_i<<" time: "<<(*G4HitT)[g4hit_i]<<std::endl; - continue; - } - - if (g4hits.find((*G4HitCellIdentifier)[g4hit_i]) == g4hits.end()) - { - //this G4 hit doesn't exist yet - one_g4hit.identifier = (*G4HitIdentifier)[g4hit_i]; - one_g4hit.cell_identifier = (*G4HitCellIdentifier)[g4hit_i]; - one_g4hit.sampling = (*G4HitSampling)[g4hit_i]; - one_g4hit.hit_time = (*G4HitT)[g4hit_i]; - //one_g4hit.hit_energy = (*G4HitE)[g4hit_i]; - //scale the hit energy with the sampling fraction - if (one_g4hit.sampling >=12 && one_g4hit.sampling <=20) - {//tile - //std::cout <<"Tile: "<<(*G4HitE)[g4hit_i]<<" "<<(*G4HitSamplingFraction)[g4hit_i]<<std::endl; - if((*G4HitSamplingFraction)[g4hit_i]) - { - one_g4hit.hit_energy = (*G4HitE)[g4hit_i]*(*G4HitSamplingFraction)[g4hit_i]; - //std::cout <<"TileE: "<<one_g4hit.hit_energy<<std::endl; - } - else one_g4hit.hit_energy = 0.; - } - else - { - //std::cout <<"LAr: "<<(*G4HitE)[g4hit_i]<<" "<<(*G4HitSamplingFraction)[g4hit_i]<<std::endl; - one_g4hit.hit_energy = (*G4HitE)[g4hit_i]/(*G4HitSamplingFraction)[g4hit_i]; - } - //one_g4hit.hit_sampfrac = (*G4HitSamplingFraction)[g4hit_i]; - //new g4hit -> insert vector with 1 element - g4hits.insert(std::pair<Long64_t, std::vector<FCS_g4hit> >(one_g4hit.cell_identifier, std::vector<FCS_g4hit>(1,one_g4hit))); - } - else - { - //G4 hit exists in this identifier -> push_back new to the vector //FCS_g4hit one_g4hit; - one_g4hit.identifier = (*G4HitIdentifier)[g4hit_i]; - one_g4hit.cell_identifier = (*G4HitCellIdentifier)[g4hit_i]; - one_g4hit.sampling = (*G4HitSampling)[g4hit_i]; - one_g4hit.hit_time = (*G4HitT)[g4hit_i]; - if (one_g4hit.sampling >=12 && one_g4hit.sampling <=20) - {//tile - //std::cout <<"Tile2: "<<(*G4HitE)[g4hit_i]<<" "<<(*G4HitSamplingFraction)[g4hit_i]<<std::endl; - if((*G4HitSamplingFraction)[g4hit_i]) - { - one_g4hit.hit_energy = (*G4HitE)[g4hit_i]*(*G4HitSamplingFraction)[g4hit_i]; - //std::cout <<"TileE2: "<<one_g4hit.hit_energy<<std::endl; - } - else one_g4hit.hit_energy = 0.; - } - else - { - //std::cout <<"LAr2: "<<(*G4HitE)[g4hit_i]<<" "<<(*G4HitSamplingFraction)[g4hit_i]<<std::endl; - one_g4hit.hit_energy = (*G4HitE)[g4hit_i]/(*G4HitSamplingFraction)[g4hit_i]; - } - //one_g4hit.hit_sampfrac = (*G4HitSamplingFraction)[g4hit_i]; - g4hits[(*G4HitCellIdentifier)[g4hit_i]].push_back(one_g4hit); - } - } + //std::cout <<truth_i<<" "<<TruthE->size()<<" "<<TruthPDG->at(truth_i)<<std::endl; + //std::cout <<"TTC2: "<<(*TTC_entrance_eta)[truth_i].size()<<std::endl;//this one has 24 layers ->ok + + new_truthE->push_back(TruthE->at(truth_i)); + new_truthPx->push_back(TruthPx->at(truth_i)); + new_truthPy->push_back(TruthPy->at(truth_i)); + new_truthPz->push_back(TruthPz->at(truth_i)); + new_truthPDG->push_back(TruthPDG->at(truth_i)); + new_truthBarcode->push_back(TruthBarcode->at(truth_i)); + //new_truthVtxBarcode->push_back(TruthVtxBarcode->at(truth_i)); + + m_newTTC_IDCaloBoundary_eta->push_back(newTTC_IDCaloBoundary_eta->at(truth_i)); + m_newTTC_IDCaloBoundary_phi->push_back(newTTC_IDCaloBoundary_phi->at(truth_i)); + m_newTTC_IDCaloBoundary_r ->push_back(newTTC_IDCaloBoundary_r->at(truth_i)); + m_newTTC_IDCaloBoundary_z ->push_back(newTTC_IDCaloBoundary_z->at(truth_i)); + m_newTTC_Angle3D ->push_back(newTTC_Angle3D->at(truth_i)); + m_newTTC_AngleEta ->push_back(newTTC_AngleEta->at(truth_i)); + + //create temporary vectors + vector<float> new_entrance_eta; + vector<float> new_entrance_phi; + vector<float> new_entrance_r; + vector<float> new_entrance_z; + vector<float> new_back_eta; + vector<float> new_back_phi; + vector<float> new_back_r; + vector<float> new_back_z; + + for (unsigned int s = 0; s < 24; s++) + { + new_entrance_eta.push_back((newTTC_entrance_eta->at(truth_i))[s]); + new_entrance_phi.push_back((newTTC_entrance_phi->at(truth_i))[s]); + new_entrance_r.push_back((newTTC_entrance_r->at(truth_i))[s]); + new_entrance_z.push_back((newTTC_entrance_z->at(truth_i))[s]); + new_back_eta.push_back((newTTC_back_eta->at(truth_i))[s]); + new_back_phi.push_back((newTTC_back_phi->at(truth_i))[s]); + new_back_r.push_back((newTTC_back_r->at(truth_i))[s]); + new_back_z.push_back((newTTC_back_z->at(truth_i))[s]); } - if (m_Debug > 1) std::cout <<" Done"<<std::endl; - - if (m_Debug > 1) std::cout <<"Reading detailed FCS hits "<< HitIdentifier->size()<<std::endl; - - for (unsigned int hit_i = 0; hit_i < HitIdentifier->size(); hit_i++) - { - if((*HitSampling)[hit_i] >=0 && (*HitSampling)[hit_i]<=25 && (*HitT)[hit_i]>m_TimingCut) - { - if (m_Debug > 1) - std::cout <<"Ignoring FCS hit, time too large: "<<hit_i<<" time: "<<(*HitT)[hit_i]<<std::endl; - continue; - } - if(hits.find((*HitCellIdentifier)[hit_i]) == hits.end()) - { - //Detailed hit doesn't exist yet - one_hit.identifier = (*HitIdentifier)[hit_i]; - one_hit.cell_identifier = (*HitCellIdentifier)[hit_i]; - one_hit.sampling = (*HitSampling)[hit_i]; - - if (one_hit.sampling >=12 && one_hit.sampling <=20) - {//tile - if((*HitSamplingFraction)[hit_i]) - { - one_hit.hit_energy = (*HitE)[hit_i]*(*HitSamplingFraction)[hit_i]; - } - else one_hit.hit_energy = 0.; - } - else - { - one_hit.hit_energy = (*HitE)[hit_i]/(*HitSamplingFraction)[hit_i]; - } - //one_hit.hit_sampfrac = (*HitSamplingFraction)[hit_i]; - one_hit.hit_time = (*HitT)[hit_i]; - one_hit.hit_x = (*HitX)[hit_i]; - one_hit.hit_y = (*HitY)[hit_i]; - one_hit.hit_z = (*HitZ)[hit_i]; - hits.insert(std::pair<Long64_t, std::vector<FCS_hit> >(one_hit.cell_identifier, std::vector<FCS_hit>(1,one_hit))); - } - else - { - //Detailed hit exists in this identifier -> push_back new to the vector - one_hit.identifier = (*HitIdentifier)[hit_i]; - one_hit.cell_identifier = (*HitCellIdentifier)[hit_i]; - one_hit.sampling = (*HitSampling)[hit_i]; - //one_hit.hit_energy = (*HitE)[hit_i]; - if (one_hit.sampling >=12 && one_hit.sampling <=20) - {//tile - if ((*HitSamplingFraction)[hit_i]) - { - one_hit.hit_energy = (*HitE)[hit_i]*(*HitSamplingFraction)[hit_i]; - } - else one_hit.hit_energy = 0.; - } - else - { - one_hit.hit_energy = (*HitE)[hit_i]/(*HitSamplingFraction)[hit_i]; - } - //one_hit.hit_sampfrac = (*HitSamplingFraction)[hit_i]; - one_hit.hit_time = (*HitT)[hit_i]; - one_hit.hit_x = (*HitX)[hit_i]; - one_hit.hit_y = (*HitY)[hit_i]; - one_hit.hit_z = (*HitZ)[hit_i]; - hits[(*HitCellIdentifier)[hit_i]].push_back(one_hit); - } - } - - if (m_Debug > 1) std::cout <<" Done"<<std::endl; - - if (m_Debug > 1) std::cout <<"ISF_HitAnalysis Check: Cells: "<<cells.size()<<" G4hits: "<<g4hits.size()<<" FCS detailed hits: "<<hits.size()<<std::endl; - - //Start matching: - Int_t iindex = 0; - for (std::map<Long64_t, FCS_cell>::iterator it = cells.begin(); it!= cells.end(); ) - { - iindex++; - // std::cout <<iindex<<std::endl; - one_matchedcell.clear(); //maybe not completely necessery, as we're not pushing_back into vectors - //set the cell part - one_matchedcell.cell = it->second; - //now look for FCS detailed hits in this cell - std::map<Long64_t, std::vector<FCS_hit> >::iterator it2 = hits.find(it->first); - if (it2!=hits.end()) - { - //std::cout <<"FCS hits found in this cell"<<std::endl; - one_matchedcell.hit = it2->second; - hits.erase(it2); //remove it - } - else - { - //no hit found for this cell - one_matchedcell.hit.clear(); //important! - } - //now look for G4hits in this cell - std::map<Long64_t, std::vector<FCS_g4hit> >::iterator it3 = g4hits.find(it->first); - if (it3 != g4hits.end()) - { - //std::cout <<"G4 hits found in this cell"<<std::endl; - one_matchedcell.g4hit = it3->second; - g4hits.erase(it3); - } - else - { - //no g4hit found for this cell - one_matchedcell.g4hit.clear();//important! - } - //std::cout <<"Erase cell"<<std::endl; - cells.erase(it++); - //std::cout <<"Insert matched object"<<std::endl; - //push_back matched cell for event jentry - oneeventcells->push_back(one_matchedcell); - //std::cout <<"Go next"<<std::endl; - - } - - //ok, cells should be empty, what about hits and g4hits? - //There could be G4hits/FCS hits for which we don't have a cell ->create a dummy empty cell with 0 energy, take the cell identifier from the hit - if (m_Debug > 1) std::cout <<"ISF_HitAnalysis Check after cells: "<<cells.size()<<" "<<g4hits.size()<<" "<<hits.size()<<std::endl; - - for (std::map<Long64_t, std::vector<FCS_hit> >::iterator it = hits.begin(); it!= hits.end();) - { - one_matchedcell.clear(); - one_matchedcell.cell.cell_identifier = it->first; - //std::cout <<"This hit didn't exist in cell: "<<it->first<<std::endl; - if (it->second.size()) - { - one_matchedcell.cell.sampling = (it->second)[0].sampling; - } - else - { - one_matchedcell.cell.sampling = -1; // - //ok, but you really shouldn't be here - std::cout <<"ERROR: You shouldn't really be here"<<std::endl; - } - one_matchedcell.cell.energy = 0.; - one_matchedcell.cell.center_x = 0.0; - one_matchedcell.cell.center_y = 0.0; - one_matchedcell.cell.center_z = 0.0; - one_matchedcell.hit = it->second; - std::map<Long64_t, std::vector<FCS_g4hit> >::iterator it3 = g4hits.find(it->first); - if (it3 != g4hits.end()) - { - one_matchedcell.g4hit = it3->second; - g4hits.erase(it3); - } - else - { - //no g4hit found for this cell - one_matchedcell.g4hit.clear(); //important! - } - hits.erase(it++); - oneeventcells->push_back(one_matchedcell); - - } - - //ok, hits should be empty, what about g4hits? - if (m_Debug > 1 )std::cout <<"ISF_HitAnalysis Check after hits: "<<cells.size()<<" "<<g4hits.size()<<" "<<hits.size()<<std::endl; - for (std::map<Long64_t, std::vector<FCS_g4hit> >::iterator it = g4hits.begin(); it!= g4hits.end();) - { - one_matchedcell.clear(); //maybe not so important - one_matchedcell.cell.cell_identifier = it->first; - if (it->second.size()) - { - one_matchedcell.cell.sampling = (it->second)[0].sampling; - } - else - { - one_matchedcell.cell.sampling = -1; // - //not really - std::cout <<"ERROR: You shouldn't really be here"<<std::endl; - } - one_matchedcell.cell.energy = 0.; - one_matchedcell.cell.center_x = 0.0; - one_matchedcell.cell.center_y = 0.0; - one_matchedcell.cell.center_z = 0.0; - one_matchedcell.g4hit = it->second; - one_matchedcell.hit.clear(); //important!! - g4hits.erase(it++); - oneeventcells->push_back(one_matchedcell); - } - - if (m_Debug > 1) std::cout <<"ISF_HitAnalysis Check after g4hits: "<<cells.size()<<" "<<g4hits.size()<<" "<<hits.size()<<std::endl; - //Final size for this event - if (m_Debug > 1) std::cout <<"ISF_HitAnalysis Matched cells size: "<<oneeventcells->size()<<std::endl; - - //Can fill the output tree already here: - total_cell_e = 0; - total_hit_e = 0; - total_g4hit_e = 0; - - for (int j=0; j<MAX_LAYER-1; j++) - { - layercells[j]->m_vector = oneeventcells->GetLayer(j); - } - - //this is for invalid cells - layercells[MAX_LAYER-1]->m_vector = oneeventcells->GetLayer(-1); - - for (int i=0; i<MAX_LAYER; i++) - { - cell_energy->at(i) = 0.0; //zero for each event! - hit_energy->at(i) = 0.0; - g4hit_energy->at(i) = 0.0; - - for (unsigned int cellindex = 0; cellindex < layercells[i]->size(); cellindex++) - { - if (i!=MAX_LAYER-1) - { - cell_energy->at(i)+=layercells[i]->m_vector.at(cellindex).cell.energy; - total_cell_e+=layercells[i]->m_vector.at(cellindex).cell.energy; - } - else - { - //don't add the energy in the invalid layer to the total energy (if there is any (shouldn't) - cell_energy->at(i)+=layercells[i]->m_vector.at(cellindex).cell.energy; //this should be here anyway - } - - //sum energy of all FCS detailed hits in this layer/cell - for (unsigned int j=0; j<layercells[i]->m_vector.at(cellindex).hit.size(); j++) - { - if (i!=MAX_LAYER-1) - { - total_hit_e+=layercells[i]->m_vector.at(cellindex).hit[j].hit_energy; - hit_energy->at(i)+=layercells[i]->m_vector.at(cellindex).hit[j].hit_energy; - } - else - { - //again, don't add invalid layer energy to the sum - hit_energy->at(i)+=layercells[i]->m_vector.at(cellindex).hit[j].hit_energy; - } - } - - //sum energy of all G4 hits in this layer/cell - for (unsigned int j=0; j<layercells[i]->m_vector.at(cellindex).g4hit.size(); j++) - { - if (i!=MAX_LAYER-1) - { - total_g4hit_e+=layercells[i]->m_vector.at(cellindex).g4hit[j].hit_energy; - g4hit_energy->at(i)+=layercells[i]->m_vector.at(cellindex).g4hit[j].hit_energy; - } - else - { - //don't add invalied layer energy to the sum - g4hit_energy->at(i)+=layercells[i]->m_vector.at(cellindex).g4hit[j].hit_energy; - } - } - } - } - - cell_energy->at(MAX_LAYER) = total_cell_e; - hit_energy->at(MAX_LAYER) = total_hit_e; - g4hit_energy->at(MAX_LAYER) = total_g4hit_e; - - //Uncomment this to get energy cross check - if (jentry % m_PrintOutFrequency ==0) std::cout <<"Energy check: event: "<<jentry<<" cell: "<<total_cell_e<<" g4hits: "<<total_g4hit_e<<" hits: "<<total_hit_e<<std::endl; - //tree gets filled - m_OutputTree->Fill(); - - }//loop over entries + //push back temporary vectors + m_newTTC_entrance_eta->push_back(new_entrance_eta); + m_newTTC_entrance_phi->push_back(new_entrance_phi); + m_newTTC_entrance_r->push_back(new_entrance_r); + m_newTTC_entrance_z->push_back(new_entrance_z); + m_newTTC_back_eta->push_back(new_back_eta); + m_newTTC_back_phi->push_back(new_back_phi); + m_newTTC_back_r->push_back(new_back_r); + m_newTTC_back_z->push_back(new_back_z); + + /* + one_truth.SetPxPyPzE((*TruthPx)[truth_i], (*TruthPy)[truth_i], (*TruthPz)[truth_i], (*TruthE)[truth_i]); + one_truth.pdgid = TruthPDG->at(truth_i); + one_truth.barcode = TruthBarcode->at(truth_i); + one_truth.vtxbarcode = TruthVtxBarcode->at(truth_i); + truthcollection->push_back(one_truth); + */ + + } + + //std::cout<<"check. after truth block"<<std::endl; + + //Now matching between cells, G4hits and G4detailed hits + //sort cells by identifier: + //clear first the containers for this event: + cells.clear(); + g4hits.clear(); + hits.clear(); + //std::cout <<"Check Original size: Cells: "<<CellIdentifier->size()<<" G4Hits: "<<G4HitCellIdentifier->size()<<" FCSHits: "<<HitCellIdentifier->size()<<std::endl; + if (m_Debug > 1) std::cout << "Reading cells..."; + + for (unsigned int cell_i = 0; cell_i < CellIdentifier->size(); cell_i++) + { + if (cells.find((*CellIdentifier)[cell_i]) == cells.end()) //doesn't exist + { + one_cell.cell_identifier = (*CellIdentifier)[cell_i]; + one_cell.sampling = (*CellSampling)[cell_i]; + one_cell.energy = (*CellE)[cell_i]; + one_cell.center_x = 0.0; //for now + one_cell.center_y = 0.0; + one_cell.center_z = 0.0; + cells.insert(std::pair<Long64_t, FCS_cell>(one_cell.cell_identifier, one_cell)); + } + else + { + //there shouldn't be a cell with the same identifier in this event + std::cout << "ISF_HitAnalysis: Same cell???? ERROR" << std::endl; + } + } - m_OutputTree->Write(); - m_Output->Close(); + if (m_Debug > 1) std::cout << " Done" << std::endl; -}; + if (m_Debug > 1) std::cout << "Reading G4hits"; -void CaloHitAna::Finish(std::vector<Int_t> settings, TString outputname) -{ -} -// //Not necessary to call now (it is done in Loop directly) - -// //This writes out the tree out of m_all_cells vector -// //settings vector specifies what branches are written out -// //if settings[0]==1 -> write all cells for all events -// //if settings[1]==1 -> write all cells for all layers(samplings) separately -// //if settings[2]==1 -> write fractions of energies in all layers - -// TFile *fout = TFile::Open(outputname,"RECREATE"); -// TTree *tout = new TTree("FCS_ParametrizationInput_TEST","Output_Matched_cell_Tree"); - -// FCS_matchedcellvector oneeventcells; //these are all matched cells in a single event -// FCS_matchedcellvector layercells[MAX_LAYER]; //these are all matched cells in a given layer in a given event - -// Float_t total_cell_e = 0.0; -// Float_t total_hit_e = 0.0; -// Float_t total_g4hit_e = 0.0; - -// std::vector<Float_t> cell_energy; -// std::vector<Float_t> hit_energy; -// std::vector<Float_t> g4hit_energy; - -// //Store energies (not fractions in each layer 0-23 (BPS-FCAL2), 24 = INVALID, 25 = TOTAL -// cell_energy.resize(MAX_LAYER+1); -// hit_energy.resize(MAX_LAYER+1); -// g4hit_energy.resize(MAX_LAYER+1); - -// for(int i =0 ; i<MAX_LAYER+1; i++) -// { -// cell_energy[i] = 0.0; -// hit_energy[i] = 0.0; -// g4hit_energy[i] = 0.0; -// } - -// //create branches in the output tree according to the settings vector -// if(settings[0]==1 || ! settings.size()) -// { -// //Write all FCS_matchedcells -// tout->Branch("AllCells", &oneeventcells); -// } -// if (settings.size()>=2 && settings[1]==1) -// { -// //write cells per layer -// for (Int_t i=0; i<MAX_LAYER; i++) -// { -// TString branchname="Sampling_"; -// branchname+=i; -// tout->Branch(branchname, &layercells[i]); -// } -// } - -// if(settings.size()>=3 && settings[2]==1) -// { -// //write also energies per layer: -// tout->Branch("cell_energy",&cell_energy); -// tout->Branch("hit_energy",&hit_energy); -// tout->Branch("g4hit_energy",&g4hit_energy); - -// //This is a duplicate of cell_energy[25] -// tout->Branch("total_cell_energy", &total_cell_e); -// tout->Branch("total_hit_energy", &total_hit_e); -// tout->Branch("total_g4hit_energy", &total_g4hit_e); -// } - -// //Loop over all events -// for(unsigned int event=0; event<m_all_cells.size(); event++) -// { - -// if (event % m_PrintOutFrequency == 0) -// std::cout <<"Reprocessing: "<<event<<std::endl; - -// total_cell_e = 0.0; -// total_hit_e = 0.0; -// total_g4hit_e = 0.0; -// oneeventcells = m_all_cells[event]; //This gets cells for particular event i (m_all_cells holds now all cells from all events (huge!) -// for (int j=0; j<MAX_LAYER-1; j++) -// layercells[j].m_vector = oneeventcells.GetLayer(j); - -// //this is for invalid cells -// layercells[MAX_LAYER-1].m_vector = oneeventcells.GetLayer(-1); - -// for (int i=0; i<MAX_LAYER; i++) -// { -// cell_energy[i] = 0.0; //zero for each event! -// hit_energy[i] = 0.0; -// g4hit_energy[i] = 0.0; -// for (unsigned int cellindex = 0; cellindex < layercells[i].size(); cellindex++) -// { -// if (i!=MAX_LAYER-1) -// { -// cell_energy[i]+=layercells[i][cellindex].cell.energy; -// total_cell_e+=layercells[i][cellindex].cell.energy; -// } -// else -// { -// //don't add the energy in the invalid layer to the total energy (if there is any (shouldn't) -// cell_energy[i]+=layercells[i][cellindex].cell.energy; //this should be here anyway -// } - -// //sum energy of all FCS detailed hits in this layer/cell -// for (unsigned int j=0; j<layercells[i][cellindex].hit.size(); j++) -// { -// if (i!=MAX_LAYER-1) -// { -// total_hit_e+=layercells[i][cellindex].hit[j].hit_energy; -// hit_energy[i]+=layercells[i][cellindex].hit[j].hit_energy; -// } -// else -// { -// //again, don't add invalid layer energy to the sum -// hit_energy[i]+=layercells[i][cellindex].hit[j].hit_energy; -// } -// } - -// //sum energy of all G4 hits in this layer/cell -// for (unsigned int j=0; j<layercells[i][cellindex].g4hit.size(); j++) -// { -// if (i!=MAX_LAYER-1) -// { -// total_g4hit_e+=layercells[i][cellindex].g4hit[j].hit_energy; -// g4hit_energy[i]+=layercells[i][cellindex].g4hit[j].hit_energy; -// } -// else -// { -// //don't add invalied layer energy to the sum -// g4hit_energy[i]+=layercells[i][cellindex].g4hit[j].hit_energy; -// } -// } -// } -// } - -// //Uncomment this to get energy cross check -// //std::cout <<"Energy check: event: "<<event<<" cell: "<<total_cell_e<<" g4hits: "<<total_g4hit_e<<" hits: "<<total_hit_e<<std::endl; -// //tree gets filled -// tout->Fill(); - -// } - -// //and written out -// tout->Write(); -// fout->Close(); - -// } + if (m_do_g4_hits) + { + for (unsigned int g4hit_i = 0; g4hit_i < G4HitIdentifier->size(); g4hit_i++) + { + if ((*G4HitSampling)[g4hit_i] >= 0 && (*G4HitSampling)[g4hit_i] <= 25 && (*G4HitT)[g4hit_i] > m_TimingCut) + { + if (m_Debug > 1) std::cout << "Ignoring G4hit, time too large: " << g4hit_i << " time: " << (*G4HitT)[g4hit_i] << std::endl; + continue; + } + + if (g4hits.find((*G4HitCellIdentifier)[g4hit_i]) == g4hits.end()) + { + //this G4 hit doesn't exist yet + one_g4hit.identifier = (*G4HitIdentifier)[g4hit_i]; + one_g4hit.cell_identifier = (*G4HitCellIdentifier)[g4hit_i]; + one_g4hit.sampling = (*G4HitSampling)[g4hit_i]; + one_g4hit.hit_time = (*G4HitT)[g4hit_i]; + //one_g4hit.hit_energy = (*G4HitE)[g4hit_i]; + //scale the hit energy with the sampling fraction + if (one_g4hit.sampling >= 12 && one_g4hit.sampling <= 20) + { //tile + //std::cout <<"Tile: "<<(*G4HitE)[g4hit_i]<<" "<<(*G4HitSamplingFraction)[g4hit_i]<<std::endl; + if ((*G4HitSamplingFraction)[g4hit_i]) + { + one_g4hit.hit_energy = (*G4HitE)[g4hit_i] * (*G4HitSamplingFraction)[g4hit_i]; + //std::cout <<"TileE: "<<one_g4hit.hit_energy<<std::endl; + } + else one_g4hit.hit_energy = 0.; + } + else + { + //std::cout <<"LAr: "<<(*G4HitE)[g4hit_i]<<" "<<(*G4HitSamplingFraction)[g4hit_i]<<std::endl; + one_g4hit.hit_energy = (*G4HitE)[g4hit_i] / (*G4HitSamplingFraction)[g4hit_i]; + } + //one_g4hit.hit_sampfrac = (*G4HitSamplingFraction)[g4hit_i]; + //new g4hit -> insert vector with 1 element + g4hits.insert(std::pair<Long64_t, std::vector<FCS_g4hit> >(one_g4hit.cell_identifier, std::vector<FCS_g4hit>(1, one_g4hit))); + } + else + { + //G4 hit exists in this identifier -> push_back new to the vector //FCS_g4hit one_g4hit; + one_g4hit.identifier = (*G4HitIdentifier)[g4hit_i]; + one_g4hit.cell_identifier = (*G4HitCellIdentifier)[g4hit_i]; + one_g4hit.sampling = (*G4HitSampling)[g4hit_i]; + one_g4hit.hit_time = (*G4HitT)[g4hit_i]; + if (one_g4hit.sampling >= 12 && one_g4hit.sampling <= 20) + { //tile + //std::cout <<"Tile2: "<<(*G4HitE)[g4hit_i]<<" "<<(*G4HitSamplingFraction)[g4hit_i]<<std::endl; + if ((*G4HitSamplingFraction)[g4hit_i]) + { + one_g4hit.hit_energy = (*G4HitE)[g4hit_i] * (*G4HitSamplingFraction)[g4hit_i]; + //std::cout <<"TileE2: "<<one_g4hit.hit_energy<<std::endl; + } + else one_g4hit.hit_energy = 0.; + } + else + { + //std::cout <<"LAr2: "<<(*G4HitE)[g4hit_i]<<" "<<(*G4HitSamplingFraction)[g4hit_i]<<std::endl; + one_g4hit.hit_energy = (*G4HitE)[g4hit_i] / (*G4HitSamplingFraction)[g4hit_i]; + } + //one_g4hit.hit_sampfrac = (*G4HitSamplingFraction)[g4hit_i]; + g4hits[(*G4HitCellIdentifier)[g4hit_i]].push_back(one_g4hit); + } + } + } + + if (m_Debug > 1) std::cout << " Done" << std::endl; + + if (m_Debug > 1) std::cout << "Reading detailed FCS hits " << HitIdentifier->size() << std::endl; + + for (unsigned int hit_i = 0; hit_i < HitIdentifier->size(); hit_i++) + { + if ((*HitSampling)[hit_i] >= 0 && (*HitSampling)[hit_i] <= 25 && (*HitT)[hit_i] > m_TimingCut) + { + if (m_Debug > 1) + std::cout << "Ignoring FCS hit, time too large: " << hit_i << " time: " << (*HitT)[hit_i] << std::endl; + continue; + } + if (hits.find((*HitCellIdentifier)[hit_i]) == hits.end()) + { + //Detailed hit doesn't exist yet + one_hit.identifier = (*HitIdentifier)[hit_i]; + one_hit.cell_identifier = (*HitCellIdentifier)[hit_i]; + one_hit.sampling = (*HitSampling)[hit_i]; + + if (one_hit.sampling >= 12 && one_hit.sampling <= 20) + { //tile + if ((*HitSamplingFraction)[hit_i]) + { + one_hit.hit_energy = (*HitE)[hit_i] * (*HitSamplingFraction)[hit_i]; + } + else one_hit.hit_energy = 0.; + } + else + { + one_hit.hit_energy = (*HitE)[hit_i] / (*HitSamplingFraction)[hit_i]; + } + //one_hit.hit_sampfrac = (*HitSamplingFraction)[hit_i]; + one_hit.hit_time = (*HitT)[hit_i]; + one_hit.hit_x = (*HitX)[hit_i]; + one_hit.hit_y = (*HitY)[hit_i]; + one_hit.hit_z = (*HitZ)[hit_i]; + hits.insert(std::pair<Long64_t, std::vector<FCS_hit> >(one_hit.cell_identifier, std::vector<FCS_hit>(1, one_hit))); + } + else + { + //Detailed hit exists in this identifier -> push_back new to the vector + one_hit.identifier = (*HitIdentifier)[hit_i]; + one_hit.cell_identifier = (*HitCellIdentifier)[hit_i]; + one_hit.sampling = (*HitSampling)[hit_i]; + //one_hit.hit_energy = (*HitE)[hit_i]; + if (one_hit.sampling >= 12 && one_hit.sampling <= 20) + { //tile + if ((*HitSamplingFraction)[hit_i]) + { + one_hit.hit_energy = (*HitE)[hit_i] * (*HitSamplingFraction)[hit_i]; + } + else one_hit.hit_energy = 0.; + } + else + { + one_hit.hit_energy = (*HitE)[hit_i] / (*HitSamplingFraction)[hit_i]; + } + //one_hit.hit_sampfrac = (*HitSamplingFraction)[hit_i]; + one_hit.hit_time = (*HitT)[hit_i]; + one_hit.hit_x = (*HitX)[hit_i]; + one_hit.hit_y = (*HitY)[hit_i]; + one_hit.hit_z = (*HitZ)[hit_i]; + hits[(*HitCellIdentifier)[hit_i]].push_back(one_hit); + } + } + + if (m_Debug > 1) std::cout << " Done" << std::endl; + + if (m_Debug > 1) std::cout << "ISF_HitAnalysis Check: Cells: " << cells.size() << " G4hits: " << g4hits.size() << " FCS detailed hits: " << hits.size() << std::endl; + + //Start matching: + Int_t iindex = 0; + for (std::map<Long64_t, FCS_cell>::iterator it = cells.begin(); it != cells.end(); ) + { + iindex++; + // std::cout <<iindex<<std::endl; + one_matchedcell.clear(); //maybe not completely necessery, as we're not pushing_back into vectors + //set the cell part + one_matchedcell.cell = it->second; + //now look for FCS detailed hits in this cell + std::map<Long64_t, std::vector<FCS_hit> >::iterator it2 = hits.find(it->first); + if (it2 != hits.end()) + { + //std::cout <<"FCS hits found in this cell"<<std::endl; + one_matchedcell.hit = it2->second; + hits.erase(it2); //remove it + } + else + { + //no hit found for this cell + one_matchedcell.hit.clear(); //important! + } + //now look for G4hits in this cell + std::map<Long64_t, std::vector<FCS_g4hit> >::iterator it3 = g4hits.find(it->first); + if (it3 != g4hits.end()) + { + //std::cout <<"G4 hits found in this cell"<<std::endl; + one_matchedcell.g4hit = it3->second; + g4hits.erase(it3); + } + else + { + //no g4hit found for this cell + one_matchedcell.g4hit.clear();//important! + } + //std::cout <<"Erase cell"<<std::endl; + cells.erase(it++); + //std::cout <<"Insert matched object"<<std::endl; + //push_back matched cell for event jentry + oneeventcells->push_back(one_matchedcell); + //std::cout <<"Go next"<<std::endl; + + } + + //ok, cells should be empty, what about hits and g4hits? + //There could be G4hits/FCS hits for which we don't have a cell ->create a dummy empty cell with 0 energy, take the cell identifier from the hit + if (m_Debug > 1) std::cout << "ISF_HitAnalysis Check after cells: " << cells.size() << " " << g4hits.size() << " " << hits.size() << std::endl; + + for (std::map<Long64_t, std::vector<FCS_hit> >::iterator it = hits.begin(); it != hits.end();) + { + one_matchedcell.clear(); + one_matchedcell.cell.cell_identifier = it->first; + //std::cout <<"This hit didn't exist in cell: "<<it->first<<std::endl; + if (it->second.size()) + { + one_matchedcell.cell.sampling = (it->second)[0].sampling; + } + else + { + one_matchedcell.cell.sampling = -1; // + //ok, but you really shouldn't be here + std::cout << "ERROR: You shouldn't really be here" << std::endl; + } + one_matchedcell.cell.energy = 0.; + one_matchedcell.cell.center_x = 0.0; + one_matchedcell.cell.center_y = 0.0; + one_matchedcell.cell.center_z = 0.0; + one_matchedcell.hit = it->second; + std::map<Long64_t, std::vector<FCS_g4hit> >::iterator it3 = g4hits.find(it->first); + if (it3 != g4hits.end()) + { + one_matchedcell.g4hit = it3->second; + g4hits.erase(it3); + } + else + { + //no g4hit found for this cell + one_matchedcell.g4hit.clear(); //important! + } + hits.erase(it++); + oneeventcells->push_back(one_matchedcell); + + } + + //ok, hits should be empty, what about g4hits? + if (m_Debug > 1 )std::cout << "ISF_HitAnalysis Check after hits: " << cells.size() << " " << g4hits.size() << " " << hits.size() << std::endl; + for (std::map<Long64_t, std::vector<FCS_g4hit> >::iterator it = g4hits.begin(); it != g4hits.end();) + { + one_matchedcell.clear(); //maybe not so important + one_matchedcell.cell.cell_identifier = it->first; + if (it->second.size()) + { + one_matchedcell.cell.sampling = (it->second)[0].sampling; + } + else + { + one_matchedcell.cell.sampling = -1; // + //not really + std::cout << "ERROR: You shouldn't really be here" << std::endl; + } + one_matchedcell.cell.energy = 0.; + one_matchedcell.cell.center_x = 0.0; + one_matchedcell.cell.center_y = 0.0; + one_matchedcell.cell.center_z = 0.0; + one_matchedcell.g4hit = it->second; + one_matchedcell.hit.clear(); //important!! + g4hits.erase(it++); + oneeventcells->push_back(one_matchedcell); + } + + if (m_Debug > 1) std::cout << "ISF_HitAnalysis Check after g4hits: " << cells.size() << " " << g4hits.size() << " " << hits.size() << std::endl; + //Final size for this event + if (m_Debug > 1) std::cout << "ISF_HitAnalysis Matched cells size: " << oneeventcells->size() << std::endl; + + //Can fill the output tree already here: + total_cell_e = 0; + total_hit_e = 0; + total_g4hit_e = 0; + + for (int j = 0; j < MAX_LAYER - 1; j++) + { + layercells[j]->m_vector = oneeventcells->GetLayer(j); + } + + //this is for invalid cells + layercells[MAX_LAYER - 1]->m_vector = oneeventcells->GetLayer(-1); + + for (int i = 0; i < MAX_LAYER; i++) + { + cell_energy->at(i) = 0.0; //zero for each event! + hit_energy->at(i) = 0.0; + g4hit_energy->at(i) = 0.0; + + for (unsigned int cellindex = 0; cellindex < layercells[i]->size(); cellindex++) + { + if (i != MAX_LAYER - 1) + { + cell_energy->at(i) += layercells[i]->m_vector.at(cellindex).cell.energy; + total_cell_e += layercells[i]->m_vector.at(cellindex).cell.energy; + } + else + { + //don't add the energy in the invalid layer to the total energy (if there is any (shouldn't) + cell_energy->at(i) += layercells[i]->m_vector.at(cellindex).cell.energy; //this should be here anyway + } + + //sum energy of all FCS detailed hits in this layer/cell + for (unsigned int j = 0; j < layercells[i]->m_vector.at(cellindex).hit.size(); j++) + { + if (i != MAX_LAYER - 1) + { + total_hit_e += layercells[i]->m_vector.at(cellindex).hit[j].hit_energy; + hit_energy->at(i) += layercells[i]->m_vector.at(cellindex).hit[j].hit_energy; + } + else + { + //again, don't add invalid layer energy to the sum + hit_energy->at(i) += layercells[i]->m_vector.at(cellindex).hit[j].hit_energy; + } + } + + //sum energy of all G4 hits in this layer/cell + for (unsigned int j = 0; j < layercells[i]->m_vector.at(cellindex).g4hit.size(); j++) + { + if (i != MAX_LAYER - 1) + { + total_g4hit_e += layercells[i]->m_vector.at(cellindex).g4hit[j].hit_energy; + g4hit_energy->at(i) += layercells[i]->m_vector.at(cellindex).g4hit[j].hit_energy; + } + else + { + //don't add invalied layer energy to the sum + g4hit_energy->at(i) += layercells[i]->m_vector.at(cellindex).g4hit[j].hit_energy; + } + } + } + } + + cell_energy->at(MAX_LAYER) = total_cell_e; + hit_energy->at(MAX_LAYER) = total_hit_e; + g4hit_energy->at(MAX_LAYER) = total_g4hit_e; + + //Uncomment this to get energy cross check + if (jentry % m_PrintOutFrequency == 0) std::cout << "Energy check: event: " << jentry << " cell: " << total_cell_e << " g4hits: " << total_g4hit_e << " hits: " << total_hit_e << std::endl; + //tree gets filled + m_OutputTree->Fill(); + + }//loop over entries + + m_OutputTree->Write(); + m_Output->Close +}; diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/CaloHitAna.h b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/CaloHitAna.h index be84063994c..302b44f9638 100644 --- a/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/CaloHitAna.h +++ b/Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/CaloHitAna.h @@ -21,12 +21,10 @@ #include <FCS_Cell.h> #include <iostream> -// Fixed size dimensions of array or collections stored in the TTree if any. -// const int CaloHitAna::MAX_LAYER =25; //number of calorimeter layers/samplings, use last (MAX_LAYER-1) for invalid cells/hits - -class CaloHitAna { - public : +class CaloHitAna +{ +public : TTree *fChain; //!pointer to the analyzed TTree or TChain Int_t fCurrent; //!current Tree number in a TChain TString fFilename; @@ -113,18 +111,13 @@ class CaloHitAna { std::vector<FCS_truth>* truthcollection = new std::vector<FCS_truth>; - Float_t total_cell_e = 0; Float_t total_hit_e = 0; Float_t total_g4hit_e = 0; - // Float_t* p_total_cell_e = &total_cell_e; - // Float_t* p_total_hit_e = &total_hit_e; - // Float_t* p_total_g4hit_e = &total_g4hit_e; - - std::vector<Float_t>* cell_energy = new std::vector<Float_t>(MAX_LAYER+1); - std::vector<Float_t>* hit_energy = new std::vector<Float_t>(MAX_LAYER+1); - std::vector<Float_t>* g4hit_energy = new std::vector<Float_t>(MAX_LAYER+1); + std::vector<Float_t>* cell_energy = new std::vector<Float_t>(MAX_LAYER + 1); + std::vector<Float_t>* hit_energy = new std::vector<Float_t>(MAX_LAYER + 1); + std::vector<Float_t>* g4hit_energy = new std::vector<Float_t>(MAX_LAYER + 1); std::vector<Float_t>* new_truthE = new std::vector<Float_t>; std::vector<Float_t>* new_truthPx = new std::vector<Float_t>; @@ -134,9 +127,6 @@ class CaloHitAna { std::vector<int>* new_truthPDG = new std::vector<int>; std::vector<int>* new_truthVtxBarcode = new std::vector<int>; - //If this won't work, we will have to change it... (memory??) - // std::vector<FCS_matchedcellvector*> m_all_cells; //hm, make it a vector of (vector of FCS_matchedcell) and have it for all 1000 events at once in memory?? - // List of branches TBranch *b_HitX; //! TBranch *b_HitY; //! @@ -169,20 +159,20 @@ class CaloHitAna { TBranch *b_G4HitSamplingFraction; //! TBranch *b_G4HitSampling; //! /* - TBranch *b_TTC_back_eta; //! - TBranch *b_TTC_back_phi; //! - TBranch *b_TTC_back_r; //! - TBranch *b_TTC_back_z; //! - TBranch *b_TTC_entrance_eta; //! - TBranch *b_TTC_entrance_phi; //! - TBranch *b_TTC_entrance_r; //! - TBranch *b_TTC_entrance_z; //! - TBranch *b_TTC_IDCaloBoundary_eta; //! - TBranch *b_TTC_IDCaloBoundary_phi; //! - TBranch *b_TTC_IDCaloBoundary_r; //! - TBranch *b_TTC_IDCaloBoundary_z; //! - TBranch *b_TTC_Angle3D; //! - TBranch *b_TTC_AngleEta; //! + TBranch *b_TTC_back_eta; //! + TBranch *b_TTC_back_phi; //! + TBranch *b_TTC_back_r; //! + TBranch *b_TTC_back_z; //! + TBranch *b_TTC_entrance_eta; //! + TBranch *b_TTC_entrance_phi; //! + TBranch *b_TTC_entrance_r; //! + TBranch *b_TTC_entrance_z; //! + TBranch *b_TTC_IDCaloBoundary_eta; //! + TBranch *b_TTC_IDCaloBoundary_phi; //! + TBranch *b_TTC_IDCaloBoundary_r; //! + TBranch *b_TTC_IDCaloBoundary_z; //! + TBranch *b_TTC_Angle3D; //! + TBranch *b_TTC_AngleEta; //! */ TBranch *b_newTTC_back_eta; //! @@ -201,7 +191,7 @@ class CaloHitAna { TBranch *b_newTTC_AngleEta; //! - CaloHitAna(TString filename="ISF_HitAnalysispion_eta1.root", TString outputname="output1.root", std::vector<Int_t> settings=std::vector<Int_t>(), Float_t timingcut=999999., Int_t debug=0, TTree *tree=0); + CaloHitAna(TString filename = "ISF_HitAnalysispion_eta1.root", TString outputname = "output1.root", std::vector<Int_t> settings = std::vector<Int_t>(), Float_t timingcut = 999999., Int_t debug = 0, TTree *tree = 0); virtual ~CaloHitAna(); virtual Int_t Cut(Long64_t entry); virtual Int_t GetEntry(Long64_t entry); @@ -209,15 +199,14 @@ class CaloHitAna { virtual void Init(TTree *tree); virtual void InitOutTree(); virtual void Loop(); - virtual void Finish(std::vector<Int_t> settings, TString outputname="output_cells.root"); virtual Bool_t Notify(); virtual void Show(Long64_t entry = -1); - void SetDebug(Int_t debug=0) {m_Debug=debug;}; - void SetTimingCut(Int_t timingcut=999999) {m_TimingCut=timingcut;}; - void SetDoAllCells(Int_t doit=0) {if(m_Settings.size()<1) m_Settings.resize(1,0);m_Settings[0]=doit;}; - void SetDoLayers(Int_t doit=0) {if(m_Settings.size()<2) m_Settings.resize(2,0);m_Settings[1]=doit;}; - void SetDoLayerSums(Int_t doit=0) {if(m_Settings.size()<3) m_Settings.resize(3,0);m_Settings[2]=doit;}; + void SetDebug(Int_t debug = 0) {m_Debug = debug;}; + void SetTimingCut(Int_t timingcut = 999999) {m_TimingCut = timingcut;}; + void SetDoAllCells(Int_t doit = 0) {if (m_Settings.size() < 1) m_Settings.resize(1, 0); m_Settings[0] = doit;}; + void SetDoLayers(Int_t doit = 0) {if (m_Settings.size() < 2) m_Settings.resize(2, 0); m_Settings[1] = doit;}; + void SetDoLayerSums(Int_t doit = 0) {if (m_Settings.size() < 3) m_Settings.resize(3, 0); m_Settings[2] = doit;}; }; #endif @@ -225,29 +214,26 @@ class CaloHitAna { #ifdef CaloHitAna_cxx CaloHitAna::CaloHitAna(TString filename, TString outputname, std::vector<Int_t> settings, Float_t timingcut, Int_t debug, TTree *tree) : fChain(0) { - fFilename= filename; + fFilename = filename; // if parameter tree is not specified (or zero), connect the file // used to generate this class and read the Tree. if (tree == 0) { - TFile *f = TFile::Open(filename,"READ"); - TString dirname=filename; - dirname+=":/ISF_HitAnalysis"; + TFile *f = TFile::Open(filename, "READ"); + TString dirname = filename; + dirname += ":/ISF_HitAnalysis"; TDirectory * dir = (TDirectory*)f->Get(dirname); - dir->GetObject("CaloHitAna",tree); - + dir->GetObject("CaloHitAna", tree); } - //tree->Print(); Init(tree); m_Settings = settings; m_Debug = debug; m_TimingCut = timingcut; m_OutputName = outputname; m_PrintOutFrequency = 100; - m_max_nentries=-1; + m_max_nentries = -1; m_Output = new TFile(m_OutputName, "RECREATE"); - m_OutputTree = new TTree("FCS_ParametrizationInput","Output_Matched_cell_Tree"); + m_OutputTree = new TTree("FCS_ParametrizationInput", "Output_Matched_cell_Tree"); InitOutTree(); - //std::cout <<"Input: "<<fFilename<<" output: "<<m_OutputName<<" debug: "<<m_Debug<<" TC: "<<m_TimingCut<<std::endl; } CaloHitAna::~CaloHitAna() @@ -279,68 +265,64 @@ Long64_t CaloHitAna::LoadTree(Long64_t entry) void CaloHitAna::InitOutTree() { - - m_OutputTree->Branch("TruthE",&new_truthE); - m_OutputTree->Branch("TruthPx",&new_truthPx); - m_OutputTree->Branch("TruthPy",&new_truthPy); - m_OutputTree->Branch("TruthPz",&new_truthPz); - m_OutputTree->Branch("TruthPDG",&new_truthPDG); - m_OutputTree->Branch("TruthBarcode",&new_truthBarcode); - m_OutputTree->Branch("TruthVtxBarcode",&new_truthVtxBarcode); //this is duplicate of what is in the truth collection, will be good to remove/hide at some point - - m_OutputTree->Branch("newTTC_back_eta",&m_newTTC_back_eta); - m_OutputTree->Branch("newTTC_back_phi",&m_newTTC_back_phi); + m_OutputTree->Branch("TruthE", &new_truthE); + m_OutputTree->Branch("TruthPx", &new_truthPx); + m_OutputTree->Branch("TruthPy", &new_truthPy); + m_OutputTree->Branch("TruthPz", &new_truthPz); + m_OutputTree->Branch("TruthPDG", &new_truthPDG); + m_OutputTree->Branch("TruthBarcode", &new_truthBarcode); + m_OutputTree->Branch("TruthVtxBarcode", &new_truthVtxBarcode); //this is duplicate of what is in the truth collection, will be good to remove/hide at some point + + m_OutputTree->Branch("newTTC_back_eta", &m_newTTC_back_eta); + m_OutputTree->Branch("newTTC_back_phi", &m_newTTC_back_phi); m_OutputTree->Branch("newTTC_back_r", &m_newTTC_back_r); m_OutputTree->Branch("newTTC_back_z", &m_newTTC_back_z); - m_OutputTree->Branch("newTTC_entrance_eta",&m_newTTC_entrance_eta); - m_OutputTree->Branch("newTTC_entrance_phi",&m_newTTC_entrance_phi); + m_OutputTree->Branch("newTTC_entrance_eta", &m_newTTC_entrance_eta); + m_OutputTree->Branch("newTTC_entrance_phi", &m_newTTC_entrance_phi); m_OutputTree->Branch("newTTC_entrance_r", &m_newTTC_entrance_r); m_OutputTree->Branch("newTTC_entrance_z", &m_newTTC_entrance_z); - m_OutputTree->Branch("newTTC_IDCaloBoundary_eta",&m_newTTC_IDCaloBoundary_eta); - m_OutputTree->Branch("newTTC_IDCaloBoundary_phi",&m_newTTC_IDCaloBoundary_phi); - m_OutputTree->Branch("newTTC_IDCaloBoundary_r",&m_newTTC_IDCaloBoundary_r); - m_OutputTree->Branch("newTTC_IDCaloBoundary_z",&m_newTTC_IDCaloBoundary_z); + m_OutputTree->Branch("newTTC_IDCaloBoundary_eta", &m_newTTC_IDCaloBoundary_eta); + m_OutputTree->Branch("newTTC_IDCaloBoundary_phi", &m_newTTC_IDCaloBoundary_phi); + m_OutputTree->Branch("newTTC_IDCaloBoundary_r", &m_newTTC_IDCaloBoundary_r); + m_OutputTree->Branch("newTTC_IDCaloBoundary_z", &m_newTTC_IDCaloBoundary_z); m_OutputTree->Branch("newTTC_Angle3D", &m_newTTC_Angle3D); - m_OutputTree->Branch("newTTC_AngleEta",&m_newTTC_AngleEta); + m_OutputTree->Branch("newTTC_AngleEta", &m_newTTC_AngleEta); //create branches in the output tree according to the settings vector - if(! m_Settings.size() || m_Settings[0]==1) - { - //Write all FCS_matchedcells - m_OutputTree->Branch("AllCells", &oneeventcells); - } - if (m_Settings.size()>=2 && m_Settings[1]==1) - { - //write cells per layer - for (Int_t i=0; i<MAX_LAYER; i++) - { - TString branchname="Sampling_"; - branchname+=i; - // std::cout<< "fail?" << std::endl; - layercells[i] = new FCS_matchedcellvector; - m_OutputTree->Branch(branchname, &layercells[i]); - } - } - if (m_Settings.size()>=3 && m_Settings[2]==1) + if (! m_Settings.size() || m_Settings[0] == 1) + { + //Write all FCS_matchedcells + m_OutputTree->Branch("AllCells", &oneeventcells); + } + if (m_Settings.size() >= 2 && m_Settings[1] == 1) + { + //write cells per layer + for (Int_t i = 0; i < MAX_LAYER; i++) { - //write also energies per layer: - m_OutputTree->Branch("cell_energy", &cell_energy); - m_OutputTree->Branch("hit_energy", &hit_energy); - m_OutputTree->Branch("g4hit_energy",&g4hit_energy); - - //This is a duplicate of cell_energy[25] - m_OutputTree->Branch("total_cell_energy", &total_cell_e); - m_OutputTree->Branch("total_hit_energy", &total_hit_e); - m_OutputTree->Branch("total_g4hit_energy", &total_g4hit_e); + TString branchname = "Sampling_"; + branchname += i; + // std::cout<< "fail?" << std::endl; + layercells[i] = new FCS_matchedcellvector; + m_OutputTree->Branch(branchname, &layercells[i]); } + } + if (m_Settings.size() >= 3 && m_Settings[2] == 1) + { + //write also energies per layer: + m_OutputTree->Branch("cell_energy", &cell_energy); + m_OutputTree->Branch("hit_energy", &hit_energy); + m_OutputTree->Branch("g4hit_energy", &g4hit_energy); + + //This is a duplicate of cell_energy[25] + m_OutputTree->Branch("total_cell_energy", &total_cell_e); + m_OutputTree->Branch("total_hit_energy", &total_hit_e); + m_OutputTree->Branch("total_g4hit_energy", &total_g4hit_e); + } // Enable/Disable recording of g4hits - if (m_Settings.size()>=4 && m_Settings[3]==1){ + if (m_Settings.size() >= 4 && m_Settings[3] == 1) { std::cout << "do g4hits!" << std::endl; m_do_g4_hits = true; } - - - } void CaloHitAna::Init(TTree *tree) @@ -385,20 +367,20 @@ void CaloHitAna::Init(TTree *tree) G4HitSamplingFraction = 0; G4HitSampling = 0; /* - TTC_back_eta = 0; - TTC_back_phi = 0; - TTC_back_r = 0; - TTC_back_z = 0; - TTC_entrance_eta = 0; - TTC_entrance_phi = 0; - TTC_entrance_r = 0; - TTC_entrance_z = 0; - TTC_IDCaloBoundary_eta = 0; - TTC_IDCaloBoundary_phi = 0; - TTC_IDCaloBoundary_r = 0; - TTC_IDCaloBoundary_z = 0; - TTC_Angle3D = 0; - TTC_AngleEta = 0; + TTC_back_eta = 0; + TTC_back_phi = 0; + TTC_back_r = 0; + TTC_back_z = 0; + TTC_entrance_eta = 0; + TTC_entrance_phi = 0; + TTC_entrance_r = 0; + TTC_entrance_z = 0; + TTC_IDCaloBoundary_eta = 0; + TTC_IDCaloBoundary_phi = 0; + TTC_IDCaloBoundary_r = 0; + TTC_IDCaloBoundary_z = 0; + TTC_Angle3D = 0; + TTC_AngleEta = 0; */ newTTC_back_eta = 0; newTTC_back_phi = 0; @@ -453,20 +435,20 @@ void CaloHitAna::Init(TTree *tree) fChain->SetBranchAddress("G4HitSamplingFraction", &G4HitSamplingFraction, &b_G4HitSamplingFraction); fChain->SetBranchAddress("G4HitSampling", &G4HitSampling, &b_G4HitSampling); /* - fChain->SetBranchAddress("TTC_back_eta", &TTC_back_eta, &b_TTC_back_eta); - fChain->SetBranchAddress("TTC_back_phi", &TTC_back_phi, &b_TTC_back_phi); - fChain->SetBranchAddress("TTC_back_r", &TTC_back_r, &b_TTC_back_r); - fChain->SetBranchAddress("TTC_back_z", &TTC_back_z, &b_TTC_back_z); - fChain->SetBranchAddress("TTC_entrance_eta", &TTC_entrance_eta, &b_TTC_entrance_eta); - fChain->SetBranchAddress("TTC_entrance_phi", &TTC_entrance_phi, &b_TTC_entrance_phi); - fChain->SetBranchAddress("TTC_entrance_r", &TTC_entrance_r, &b_TTC_entrance_r); - fChain->SetBranchAddress("TTC_entrance_z", &TTC_entrance_z, &b_TTC_entrance_z); - fChain->SetBranchAddress("TTC_IDCaloBoundary_eta", &TTC_IDCaloBoundary_eta, &b_TTC_IDCaloBoundary_eta); - fChain->SetBranchAddress("TTC_IDCaloBoundary_phi", &TTC_IDCaloBoundary_phi, &b_TTC_IDCaloBoundary_phi); - fChain->SetBranchAddress("TTC_IDCaloBoundary_r", &TTC_IDCaloBoundary_r, &b_TTC_IDCaloBoundary_r); - fChain->SetBranchAddress("TTC_IDCaloBoundary_z", &TTC_IDCaloBoundary_z, &b_TTC_IDCaloBoundary_z); - fChain->SetBranchAddress("TTC_Angle3D", &TTC_Angle3D, &b_TTC_Angle3D); - fChain->SetBranchAddress("TTC_AngleEta", &TTC_AngleEta, &b_TTC_AngleEta); + fChain->SetBranchAddress("TTC_back_eta", &TTC_back_eta, &b_TTC_back_eta); + fChain->SetBranchAddress("TTC_back_phi", &TTC_back_phi, &b_TTC_back_phi); + fChain->SetBranchAddress("TTC_back_r", &TTC_back_r, &b_TTC_back_r); + fChain->SetBranchAddress("TTC_back_z", &TTC_back_z, &b_TTC_back_z); + fChain->SetBranchAddress("TTC_entrance_eta", &TTC_entrance_eta, &b_TTC_entrance_eta); + fChain->SetBranchAddress("TTC_entrance_phi", &TTC_entrance_phi, &b_TTC_entrance_phi); + fChain->SetBranchAddress("TTC_entrance_r", &TTC_entrance_r, &b_TTC_entrance_r); + fChain->SetBranchAddress("TTC_entrance_z", &TTC_entrance_z, &b_TTC_entrance_z); + fChain->SetBranchAddress("TTC_IDCaloBoundary_eta", &TTC_IDCaloBoundary_eta, &b_TTC_IDCaloBoundary_eta); + fChain->SetBranchAddress("TTC_IDCaloBoundary_phi", &TTC_IDCaloBoundary_phi, &b_TTC_IDCaloBoundary_phi); + fChain->SetBranchAddress("TTC_IDCaloBoundary_r", &TTC_IDCaloBoundary_r, &b_TTC_IDCaloBoundary_r); + fChain->SetBranchAddress("TTC_IDCaloBoundary_z", &TTC_IDCaloBoundary_z, &b_TTC_IDCaloBoundary_z); + fChain->SetBranchAddress("TTC_Angle3D", &TTC_Angle3D, &b_TTC_Angle3D); + fChain->SetBranchAddress("TTC_AngleEta", &TTC_AngleEta, &b_TTC_AngleEta); */ fChain->SetBranchAddress("newTTC_back_eta", &newTTC_back_eta, &b_newTTC_back_eta); fChain->SetBranchAddress("newTTC_back_phi", &newTTC_back_phi, &b_newTTC_back_phi); @@ -509,7 +491,7 @@ Int_t CaloHitAna::Cut(Long64_t entry) // This function may be called from Loop. // returns 1 if entry is accepted. // returns -1 otherwise. - std::cout <<entry<<std::endl; + std::cout << entry << std::endl; return 1; } #endif // #ifdef CaloHitAna_cxx -- GitLab