//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// #include "MyReconstruction.h" #include "TFile.h" using namespace std; //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Constructor //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// MyReconstruction::MyReconstruction() { std::cout << "Entering constructor." << std::endl; char outputFileName[10000]; sprintf(outputFileName, "/fermion/users/krishna/ParticleDensities_%d.root", sDTH.iDate); outputFile = new TFile(outputFileName, "RECREATE"); if (outputFile->IsOpen()) { std::cout << "Output file opened successfully: " << outputFileName << std::endl; } else { std::cerr << "Failed to open output file: " << outputFileName << std::endl; } for (int det = 0; det < Det; ++det) { TString histName = Form("hist_%d_%d", det, sDTH.iDate); TH1F* hist = new TH1F(histName, "Particle Density", 1010, 0.0, 1000.0); hParticleDensity[det] = hist; if (hist) { std::cout << "Created histogram for detector " << det << std::endl; hist->Write(); // Write the histogram to the output file } else { std::cerr << "Failed to create histogram for detector " << det << std::endl; } } std::cout << "Finished histogram creation." << std::endl; std::cout << "Exiting constructor." << std::endl; /*std::cout << "Entering constructor." << std::endl; char outputFileName[10000]; sprintf(outputFileName, "/fermion/users/krishna/ParticleDensities_%d.root", sDTH.iDate); outputFile = new TFile(outputFileName, "RECREATE"); if (outputFile->IsOpen()) { std::cout << "Output file opened successfully: " << outputFileName << std::endl; } else { std::cerr << "Failed to open output file: " << outputFileName << std::endl; } for (int det = 0; det < Det; ++det) { TString histName = Form("hist_%d_%d", det, sDTH.iDate); TH1F* hist = new TH1F(histName, "Particle Density", 1010, 0.0, 1000.0); hParticleDensity[det] = hist; if (hist) { std::cout << "Created histogram for detector " << det << std::endl; } else { std::cerr << "Failed to create histogram for detector " << det << std::endl; } } std::cout << "Finished histogram creation." << std::endl; std::cout << "Exiting constructor." << std::endl;*/ } // ::MyReconstruction() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Destructor //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// MyReconstruction::~MyReconstruction() { std::cout << "Entering destructor." << std::endl; if (outputFile && outputFile->IsOpen() && outputFile->IsWritable()) { std::cout << "Output file is open and writable." << std::endl; outputFile->cd(); for (auto& entry : hParticleDensity) { if (entry.second != nullptr) { entry.second->Write(); std::cout << "Wrote histogram: " << entry.second->GetName() << std::endl; } } outputFile->Close(); std::cout << "Output file closed." << std::endl; } else { std::cerr << "Output file is not open or not writable." << std::endl; } std::cout << "Exiting destructor." << std::endl; } // ::~MyReconstruction() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Setting inputs //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::SetInputs(struct ControlFile control) { inputs = control; Initialise(); } // ::SetInputs() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Initialise variables here //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::Initialise() { strcpy(inputs.Path[0], "/fermion/users/krishna/RecoIn/CALDB/STAGE1"); strcpy(inputs.Path[1], "/fermion/GRAPES-3/ROOTOUT/SC"); strcpy(inputs.Path[2], "/fermion/users/krishna/RecoIn/SCHEADER"); inputs.StartDate = 20140101; // Start date inputs.StopDate = 20140101; sSC.ndet = 10000; for (int i = 0; i < sSC.ndet; ++i) { sSC.detno[i] = i + 1; } CurvatureCorrection = inputs.dictionary.RECONSTRUCTION_CURVATURE_CORRECTION; //ReadDetectorCoordinates(); RefSHIndex = 0; } // ::Initialise() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Core processing //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::Reconstruct() { if (sSC.ndet <= 0) { std::cerr << "Error: Number of detectors (sSC.ndet) is not properly initialized." << std::endl; return; // Exit the method or handle the error as appropriate } std::cout << "Number of detectors: " << sSC.ndet << std::endl; for(int DayNo(inputs.StartDayNo); DayNo <= inputs.StopDayNo; DayNo++) { GetDate(DayNo); if(0 == OpenScRootFile()) { continue; } if(0 == OpenScHeaderFile()) { continue; } if(0 == OpenDBRootFile()) { continue; } // If curvature corrected angle is needed if(inputs.dictionary.RECONSTRUCTION_CURVATURE_CORRECTION) { if(0 == OpenReconstructedFile()) { continue; } } if(0 == ValidateEntries()) { continue; } SetOutputs(); ProcessEventTree(); FinaliseDay(); } // for DayNo <= inputs.StopDayNo // Delete the histograms for (auto& entry : hParticleDensity) { if (entry.second != nullptr) { delete entry.second; entry.second = nullptr; // Set the pointer to nullptr after deleting } } hParticleDensity.clear(); // Clear the map } // ::Reconstruct() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Get date for day number //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::GetDate(int dayno) { cDTH.SetCurrentDayNo(inputs.Year, dayno); cDTH.GetAllParameters(sDTH); } // ::GetDate() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Open gr3dbYYYYMMDD.root files /// RETURN STATUS /// 0 : FAILURE /// 1 : SUCCESS //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// bool MyReconstruction::OpenDBRootFile() { FileFoundFlag = 0; // getting the input file from the path for(unsigned short PathNo(0); PathNo < inputs.NPath; PathNo++) { sprintf(ReadName, "%s/caldb%d.root", inputs.Path[PathNo], sDTH.iDate); Check = fopen(ReadName, "r"); // Check if this file is opened if(Check == NULL) { FileFoundFlag = 0; } // if SC2ROOT else { FileFoundFlag = 1; break; } } // for PathNo 4 // checking the status of input file if(FileFoundFlag == 0) { printf(COLOR_ERROR "\n IN %s::%s(%d)" COLOR_RESET, __FILE__, __func__, __LINE__); printf(COLOR_ERROR "\n **** ERROR : OPENING FILE 'caldb%d.root' ****\n\n" COLOR_RESET, sDTH.iDate); return 0; } // if FileFoundFlag else { fclose(Check); //printf("\n\t%s", ReadName); DB = new TFile(ReadName); } // process this file if(DB->IsOpen()) { scdbtree = (TTree *)DB->Get("caldb"); DBRunEntries = scdbtree->GetEntries(); //GGetDBaseBranches(scdbtree, sDBRP, sDBDI, sDBPED, sDBGRHP, sDBGRLP, sDBLGR, sDBMC, sDBMCL, sDBCG, sDBTC, sDBTZ, sDBDQF, sDBWP); cGDB = new GetDBaseBranches(scdbtree); printf(COLOR_INPUT "\n %4d %10.3f MBytes %s [Input]" COLOR_RESET, DBRunEntries, GetFileSize(ReadName) / 1.E6, ReadName); ReadDBData(); return 1; } else { return 0; } } // ::OpenDBRootFile() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Read only run index DB data //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::ReadDBData() { ResetDBData(); for(int iEntry(0); iEntry < DBRunEntries; iEntry++) { cGDB->ResetData(); scdbtree->GetEntry(iEntry); DBIndex[cGDB->sDBRP.runno] = iEntry; //printf("\n iEntry = %4d, cGDB->sDBRP.runno = %6d", iEntry, cGDB->sDBRP.runno); } // for iEntry < DBEntries } // ::ReadDBData() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Resetting DB data //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::ResetDBData() { PreviousRun = -1; CurrentRun = -1; for(int i(0); i < MAXRUN; i++) { DBIndex[i] = -1; } // for i < MAXRUN } // ::ResetDBData() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Get DB data for the required run number //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::GetDBData(int runno) { cGDB->ResetData(); scdbtree->GetEntry(DBIndex[runno]); cGDB->GetDetectorwise(); //DisplayDBData(); } // ::GetDBData() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Display DB data //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::DisplayDBData() { printf("\n DBIndex[%d] = %3d", CurrentRun, DBIndex[CurrentRun]); printf("\n cGDB->dDBRP.runno = %6d", cGDB->dDBRP.runno); printf("\n cGDB->dDBRP.firsteventdate = %8d", cGDB->dDBRP.firsteventdate); printf("\n cGDB->dDBRP.firsteventtime1 = %6d", cGDB->dDBRP.firsteventdate); printf("\n cGDB->dDBRP.runtime = %f", cGDB->dDBRP.runtime); printf("\n cGDB->dDBDI.ndet = %d", cGDB->dDBDI.ndet); for(int iDet(1); iDet < MAXD; iDet++) { /* printf("\n %3d %2d %10.3f %10.3f %10.3f", iDet, cGDB->dDBDI.status[iDet], cGDB->dDBDI.detX[iDet], cGDB->dDBDI.detY[iDet], cGDB->dDBDI.detZ[iDet]); printf("\n %3d %10.3f %10.3f %10.3f %10.3f %2d", iDet, cGDB->dDBPED.pedhhmean[iDet], cGDB->dDBPED.pedhhrms[iDet], cGDB->dDBPED.pedhhpeak[iDet], cGDB->dDBPED.pedhhsigma[iDet], cGDB->dDBPED.pedhhfitflag[iDet]); printf("\n %10.3f %10.3f %10.3f %10.3f %2d", cGDB->dDBPED.pedhlmean[iDet], cGDB->dDBPED.pedhlrms[iDet], cGDB->dDBPED.pedhlpeak[iDet], cGDB->dDBPED.pedhlsigma[iDet], cGDB->dDBPED.pedhlfitflag[iDet]); printf("\n %10.3f %10.3f %10.3f %10.3f %2d", cGDB->dDBPED.pedlhmean[iDet], cGDB->dDBPED.pedlhrms[iDet], cGDB->dDBPED.pedlhpeak[iDet], cGDB->dDBPED.pedlhsigma[iDet], cGDB->dDBPED.pedlhfitflag[iDet]); printf("\n %10.3f %10.3f %10.3f %10.3f %2d", cGDB->dDBPED.pedllmean[iDet], cGDB->dDBPED.pedllrms[iDet], cGDB->dDBPED.pedllpeak[iDet], cGDB->dDBPED.pedllsigma[iDet], cGDB->dDBPED.pedllfitflag[iDet]); printf("\n %3d %10.3f %2d", iDet, cGDB->dDBGRHP.adcgr[iDet], cGDB->dDBGRHP.adcgrDQF[iDet]); printf("\n %3d %10.3f %10.3f %10.3f %10.3f", iDet, cGDB->dDBMC.muadcpeak[iDet], cGDB->dDBMC.muadcmedian[iDet], cGDB->dDBMC.mutdcpeak[iDet], cGDB->dDBMC.mutdcmedian[iDet]); printf("\n %3d %10.3f %10.3f %10.3f %10.3f", iDet, cGDB->dDBMCL.muadcpeak[iDet], cGDB->dDBMCL.muadcmedian[iDet], cGDB->dDBMCL.mutdcpeak[iDet], cGDB->dDBMCL.mutdcmedian[iDet]); printf("\n %3d %10.3f %2d", iDet, cGDB->dDBCG.computedadc[iDet], cGDB->dDBCG.adcDQF[iDet]); printf("\n %3d %10.3f %10.3f %10.3f", iDet, cGDB->dDBTC.a0[iDet], cGDB->dDBTC.a1[iDet], cGDB->dDBTC.a2[iDet]); printf("\n %3d %10.3f %10.3f %10.3f", iDet, cGDB->dDBTZ.tdc0[iDet], cGDB->dDBTZ.tdc0rms[iDet], cGDB->dDBTZ.tdc0err[iDet]); printf("\n %3d %2d %2d", iDet, cGDB->dDBDQF.ADCFlag[iDet], cGDB->dDBDQF.TDCFlag[iDet]); */ } // for iDet < MAXD printf("\n %10.3f %10.3f %10.3f", cGDB->dDBWP.tempin, cGDB->dDBWP.tempout, cGDB->dDBWP.phPa); getchar(); } // ::DisplayDBData() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Open scYYYYMMDD.root file /// RETURN STATUS /// 0 : FAILURE /// 1 : SUCCESS //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// bool MyReconstruction::OpenScRootFile() { FileFoundFlag = 0; // getting the input file from the path for(unsigned short PathNo(0); PathNo < inputs.NPath; PathNo++) { sprintf(ReadName, "%s/sc%d.root", inputs.Path[PathNo], sDTH.iDate); Check = fopen(ReadName, "r"); // Check if this file is opened if(Check == NULL) { FileFoundFlag = 0; } // if SC2ROOT else { FileFoundFlag = 1; break; } } // for PathNo 4 // checking the status of input file if(FileFoundFlag == 0) { printf(COLOR_ERROR "\n IN %s::%s(%d)" COLOR_RESET, __FILE__, __func__, __LINE__); printf(COLOR_ERROR "\n **** ERROR : OPENING FILE 'sc%d.root' ****\n\n" COLOR_RESET, sDTH.iDate); return 0; } // if FileFoundFlag else { fclose(Check); //printf("\n\t%s", ReadName); SC2ROOT = new TFile(ReadName); } // process this file if(SC2ROOT->IsOpen()) { GetTrees(); printf(COLOR_INPUT "\n %4d %10d %10.3f MBytes %s [Input]" COLOR_RESET, RunEntries, EventEntries, GetFileSize(ReadName) / 1.E6, ReadName); ReadRunData(); return 1; } else { return 0; } } // ::OpenScRootFile() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Open scYYYYMMDD.root file /// RETURN STATUS /// 0 : FAILURE /// 1 : SUCCESS //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// bool MyReconstruction::OpenReconstructedFile() { FileFoundFlag = 0; // getting the input file from the path for(unsigned short PathNo(0); PathNo < inputs.NPath; PathNo++) { sprintf(ReadName, "%s/sh%d.root", inputs.Path[PathNo], sDTH.iDate); Check = fopen(ReadName, "r"); // Check if this file is opened if(Check == NULL) { FileFoundFlag = 0; } // if SC2ROOT else { FileFoundFlag = 1; break; } } // for PathNo 4 // checking the status of input file if(FileFoundFlag == 0) { printf(COLOR_ERROR "\n IN %s::%s(%d)" COLOR_RESET, __FILE__, __func__, __LINE__); printf(COLOR_ERROR "\n **** ERROR : OPENING FILE 'sh%d.root' ****\n\n" COLOR_RESET, sDTH.iDate); return 0; } // if FileFoundFlag else { fclose(Check); //printf("\n\t%s", ReadName); SH2ROOT = new TFile(ReadName); } // process this file if(SH2ROOT->IsOpen()) { //GetTrees(); shtree = (TTree *) SH2ROOT->Get("sh"); if(!shtree) { printf(COLOR_ERROR "\n **** ERROR : NO TREE FOUND IN %s " COLOR_RESET , ReadName); return 0; } cShP = new SHPAR3(); cShP->GetShParData(shtree); printf(COLOR_INPUT "\n %s [input]" COLOR_RESET, ReadName); return 1; } else { return 0; } } // ::OpenReconstructedFile() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Open ScHeaderYYYYMMDD.root file /// RETURN STATUS /// 0 : FAILURE /// 1 : SUCCESS //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// bool MyReconstruction::OpenScHeaderFile() { FileFoundFlag = 0; // getting the input file from the path for(unsigned short PathNo(0); PathNo < inputs.NPath; PathNo++) { sprintf(ReadName, "%s/ScHeader%d.root", inputs.Path[PathNo], sDTH.iDate); Check = fopen(ReadName, "r"); // Check if this file is opened if(Check == NULL) { FileFoundFlag = 0; } // if SC2ROOT else { FileFoundFlag = 1; break; } } // for PathNo 4 // checking the status of input file if(FileFoundFlag == 0) { printf(COLOR_ERROR "\n IN %s::%s(%d)" COLOR_RESET, __FILE__, __func__, __LINE__); printf(COLOR_ERROR "\n **** ERROR : OPENING FILE 'ScHeader%d.root' ****\n\n" COLOR_RESET, sDTH.iDate); return 0; } // if FileFoundFlag else { fclose(Check); //printf("\n\t%s", ReadName); SCH2ROOT = new TFile(ReadName); } // process this file if(SCH2ROOT->IsOpen()) { //GetTrees(); schtree = (TTree *) SCH2ROOT->Get("ScHeader"); if(!schtree) { printf(COLOR_ERROR "\n **** ERROR : NO TREE FOUND IN %s " COLOR_RESET , ReadName); return 0; } ScHEntries = schtree->GetEntries(); //cShP = new SHPAR3(); //cShP->GetShParData(shtree); GetScHeaderBranches(); printf(COLOR_INPUT "\n %s [input]" COLOR_RESET, ReadName); return 1; } else { return 0; } } // ::OpenScHeaderFile() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Read only run index DB data //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::ReadRunData() { //printf("\n %s:%s()", __FILE__, __func__); ResetRunData(); for(int iEntry(0); iEntry < RunEntries; iEntry++) { scruntree->GetEntry(iEntry); RunIndex[sSRP.runno] = iEntry; //printf("\n iEntry = %4d, Runno = %6d", iEntry, sSRP.runno); //getchar(); } // for iEntry < DBEntries } // ::ReadRunData() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Resetting DB data //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::ResetRunData() { for(int i(0); i < MAXRUN; i++) { RunIndex[i] = -1; } // for i < MAXRUN } // ::ResetRunData() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Reset DB data for one entry //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::ResetRunData(int runno) { runno--; runno++; } // ::ResetRunData() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Get DB data for the required run number //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::GetRunData(int runno) { ResetRunData(runno); scruntree->GetEntry(RunIndex[runno]); //printf("\n RunIndex[%6d] = %3d", runno, RunIndex[runno]); //getchar(); } // ::GetRunData() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Get both TTree //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::GetTrees() { // getting the run tree scruntree = (TTree *)SC2ROOT->Get("scruntree"); // linking variables GGetScRunTrBranches(scruntree, sSRP, sPED, sSMP, sSCMI, sSDI); RunEntries = scruntree->GetEntries(); // getting the event tree scevtree = (TTree *)SC2ROOT->Get("scevtree"); // linking variables GGetScEvTrBranches(scevtree, sEVTH, sSC, sTMH, sMonTDC); // getting total entries EventEntries = scevtree->GetEntries(); } // ::GetTrees() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Finalise day //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::FinaliseDay() { scevtree->Delete(); SC2ROOT->Close(); scdbtree->Delete(); DB->Close(); delete cGDB; if(inputs.dictionary.RECONSTRUCTION_RECONSTRUCT || inputs.dictionary.RECONSTRUCTION_CURVATURE_CORRECTION) { shRootFile->cd(); sh->Write(); sh->Delete(); shRootFile->Close(); } if(inputs.dictionary.RECONSTRUCTION_GENERATE_PARTICLE_TOF) { ptRootFile->cd(); pt->Write(); pt->Delete(); ptRootFile->Close(); } } // ::FinaliseDay() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Process event tree //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::ProcessEventTree() { bool ENTRYFOUND(0); printf("\n"); ProgressBar bar; //EventEntries = 1000; if(inputs.SplitIndexFlag && (inputs.dictionary.RECONSTRUCTION_RECONSTRUCT || inputs.dictionary.RECONSTRUCTION_CURVATURE_CORRECTION)) { DivisionEntry = EventEntries / inputs.SplitIndexMax; RemainderEntry = EventEntries % inputs.SplitIndexMax; StartEntry = (inputs.SplitIndexMin - 1) * DivisionEntry; StopEntry = inputs.SplitIndexMin * DivisionEntry; if(inputs.SplitIndexMin == inputs.SplitIndexMax) { StopEntry += RemainderEntry; } /* printf("\n YES "); printf("\n SplitIndexMin %d, SplitIndexMax %d", inputs.SplitIndexMin, inputs.SplitIndexMax); printf("\n EventEntries %d", EventEntries); printf("\n DivisionEntry %d", DivisionEntry); printf("\n RemainderEntry %d", RemainderEntry); printf("\n StartEntry %d", StartEntry); printf("\n StopEntry %d", StopEntry); getchar(); */ printf(COLOR_WARNING "\n **** WARNING : PROCESSING EVENTS FROM %d to %d ****\n\n" COLOR_RESET, StartEntry, StopEntry); } else { StartEntry = 0; StopEntry = EventEntries; } //StopEntry = 1000; RefSHIndex = 0; for(int iEntry(StartEntry); iEntry < StopEntry; iEntry++) { bar.Send(iEntry - StartEntry, (StopEntry - StartEntry)); scevtree->GetEntry(iEntry); schtree->GetEntry(iEntry); //if(!EventsMatch()) continue; // If Sc and ScHeader events mis-match CurrentEntry = iEntry; // Take only air showers if(sEVTH.trigger != 2) { continue; } // if(sEVTH.trigger != 2) // If split index is given during curvature correction if(inputs.SplitIndexFlag && inputs.dictionary.RECONSTRUCTION_CURVATURE_CORRECTION && ENTRYFOUND == 0) { ENTRYFOUND = FindShEntry(); if(ENTRYFOUND==0) continue; } ProcessEvent(); RefSHIndex ++; } // for iEntry < EventEntries } // ::ProcessEventTree() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Sc events and ScHeader events matchin //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// bool MyReconstruction::EventsMatch() { if((sEVTH.runno == p_ScRunNo) && (sEVTH.eventno == p_ScEvNo)) { AssignCorrectTime(); return 1; } return 0; } // ::EventsMatch() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Assign correct time stamp from sc header //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::AssignCorrectTime() { sEVTH.evdate = p_ScEvDate; sEVTH.evtime1 = p_ScEvTime1; sEVTH.evtime2 = p_ScEvTime2; } // ::AssignCorrectTime() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Process one event here //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::ProcessEvent() { // Access the sh tree for NKG parameters if(inputs.dictionary.RECONSTRUCTION_CURVATURE_CORRECTION) { cShP->ResetData(); shtree->GetEntry(RefSHIndex); //if(cShP->Data.NKGFitFlag <1) return; //GetNKGInfo(); // Validate the mapping between the SC and SH events if(sEVTH.runno != cShP->Data.RunNo || sEVTH.eventno != cShP->Data.EventNo) { printf(COLOR_ERROR "\n **** ERROR : SC - SH mapping mis-match **** \n"); printf(COLOR_ERROR "\n SC RUN NO : %4d EventNo = %5d", sEVTH.runno, sEVTH.eventno); printf(COLOR_ERROR "\n SH RUN NO : %4d EventNo = %5d \n", cShP->Data.RunNo, cShP->Data.EventNo); return; } } CurrentRun = sEVTH.runno; // If run number changes then only get the DB entry // This method saves time and resources if(CurrentRun > PreviousRun) { GetDBData(CurrentRun); GetRunData(CurrentRun); ValidateRunno(CurrentRun); CalculateArrayCenter(); } // if(CurrentRun > PreviousRun) PreviousRun = CurrentRun; FillEventData(); } // ::ProcessEvent() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Process event data //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::FillEventData() { if(inputs.dictionary.RECONSTRUCTION_RECONSTRUCT) { ReconstructionCore::ResetShowerData(&sSD); UpdateWeights(&sSD, 0); //Reset PSum values ResetPSum(); if(inputs.dictionary.RECONSTRUCTION_RECONSTRUCT_EVENODD) { ReconstructionCore::ResetShowerData(&sSDOdd); ReconstructionCore::ResetShowerData(&sSDEven); UpdateWeights(&sSDOdd, 0); UpdateWeights(&sSDEven, 0); } // EVENODD if(inputs.dictionary.RECONSTRUCTION_RECONSTRUCT_LEFTRIGHT) { ReconstructionCore::ResetShowerData(&sSDLeft); ReconstructionCore::ResetShowerData(&sSDRight); UpdateWeights(&sSDLeft, 0); UpdateWeights(&sSDRight, 0); } // LEFTRIGHT } // RECONSTRUCT if(inputs.dictionary.RECONSTRUCTION_CURVATURE_CORRECTION) { ReconstructionAngle::ResetShowerData(&sSD); UpdateWeights(&sSD, 0); if(inputs.dictionary.RECONSTRUCTION_RECONSTRUCT_EVENODD) { ReconstructionAngle::ResetShowerData(&sSDOdd); ReconstructionAngle::ResetShowerData(&sSDEven); UpdateWeights(&sSDOdd, 0); UpdateWeights(&sSDEven, 0); } // EVENODD if(inputs.dictionary.RECONSTRUCTION_RECONSTRUCT_LEFTRIGHT) { ReconstructionAngle::ResetShowerData(&sSDLeft); ReconstructionAngle::ResetShowerData(&sSDRight); UpdateWeights(&sSDLeft, 0); UpdateWeights(&sSDRight, 0); } // LEFTRIGHT } // CURVATURE_CORRECTION if(inputs.dictionary.RECONSTRUCTION_GENERATE_PARTICLE_TOF) { ResetPTData(&sPT); } //printf("\n\n\n"); //printf("\n runno = %6d, eventno = %4d, ndet = %3d", sEVTH.runno, sEVTH.eventno, sSC.ndet); //getchar(); ProcessDetectors(1); // For curvature correction if(inputs.dictionary.RECONSTRUCTION_CURVATURE_CORRECTION) { //printf("\n runno = %6d, eventno = %4d, ndet = %3d", sEVTH.runno, sEVTH.eventno, sSC.ndet); //getchar(); cAng = new ReconstructionAngle(); cAng->ArrivalTimeOffset(&sSD); cAng->SetInputs(&sSD, 0, 1, 0); // Without NKG cAng->Reconstruct(&sSD, 1); // 1 for PMT uncorrected data if(inputs.dictionary.RECONSTRUCTION_RECONSTRUCT_EVENODD) { // Even-Odd method without NKG cAng->ArrivalTimeOffset(&sSDEven); cAng->SetInputs(&sSDEven, 0, 1, 0); cAng->Reconstruct(&sSDEven, 1); cAng->ArrivalTimeOffset(&sSDOdd); cAng->SetInputs(&sSDOdd, 0, 1, 0); cAng->Reconstruct(&sSDOdd, 1); } // Even-Odd if(inputs.dictionary.RECONSTRUCTION_RECONSTRUCT_LEFTRIGHT) { ProcessDetectors(2); // Left-Right method without NKG cAng->ArrivalTimeOffset(&sSDLeft); cAng->SetInputs(&sSDLeft, 0, 1, 0); cAng->Reconstruct(&sSDLeft, 1); cAng->ArrivalTimeOffset(&sSDRight); cAng->SetInputs(&sSDRight, 0, 1, 0); cAng->Reconstruct(&sSDRight, 1); } // Left-Right FillShTree(); delete cAng; } // Curvature Correction // Reconstrut here if(inputs.dictionary.RECONSTRUCTION_RECONSTRUCT) { // All detectors method //PrintShowerData(); //getchar(); cR = new ReconstructionCore(); cR->ArrivalTimeOffset(&sSD); cR->SetInputs(&sSD, 0, 1, inputs.dictionary.RECONSTRUCTION_RECONSTRUCT_NKGFIT); cR->Reconstruct(&sSD, 1); // 1 for PMT uncorrected data // Even-Odd and Left-Right methods are only for all high-gain PMTs if(inputs.dictionary.RECONSTRUCTION_RECONSTRUCT_PMT == 1) { if(inputs.dictionary.RECONSTRUCTION_RECONSTRUCT_EVENODD) { // Even method without NKG fit cR->ArrivalTimeOffset(&sSDEven); cR->SetInputs(&sSDEven, 0, 1, 0); cR->Reconstruct(&sSDEven, 1); // 1 for PMT uncorrected data // Odd method without NKG fit cR->ArrivalTimeOffset(&sSDOdd); cR->SetInputs(&sSDOdd, 0, 1, 0); cR->Reconstruct(&sSDOdd, 1); // 1 for PMT uncorrected data } // EVENODD if(inputs.dictionary.RECONSTRUCTION_RECONSTRUCT_LEFTRIGHT) { ProcessDetectors(2); // For Left-Right arrray /* printf("\n COGX = %10.3f, COGY = %10.3f", sSD.COGX, sSD.COGY); for(int iDet(1); iDet < MAXD; iDet++) { if(sSDLeft.HitFlag[iDet] == 1) { printf("\n iDet Left %4d", iDet); } } printf("\n-----------------------------------------------------"); for(int iDet(1); iDet < MAXD; iDet++) { if(sSDRight.HitFlag[iDet] == 1) { printf("\n iDet Right %4d", iDet); } } printf("\n====================================================="); getchar(); */ // Left array fit without NKG fit cR->ArrivalTimeOffset(&sSDLeft); cR->SetInputs(&sSDLeft, 0, 1, 0); cR->Reconstruct(&sSDLeft, 1); // For PMT uncorrected data // Right array fit without NKG fit cR->ArrivalTimeOffset(&sSDRight); cR->SetInputs(&sSDRight, 0, 1, 0); cR->Reconstruct(&sSDRight, 1); // For PMT uncorrected data } // LEFTRIGHT } // RECONSTRUCTION_RECONSTRUCT_PMT //PrintReconstructedData(); FillShTree(); delete cR; } // RECONSTRUCTION_RECONSTRUCT if(inputs.dictionary.RECONSTRUCTION_GENERATE_PARTICLE_TOF) { FillPTTree(); } // RECONSTRUCTION_GENERATE_PARTICLE_TOF } // ::FillEventData() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Entry found //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// bool MyReconstruction::FindShEntry() { for(int ish(0); ish < shtree->GetEntries(); ish++ ) { shtree->GetEntry(ish); if(sEVTH.runno == cShP->Data.RunNo && sEVTH.eventno == cShP->Data.EventNo) { RefSHIndex = ish; return 1; cout << sEVTH.runno <<" " << sEVTH.eventno << endl; cout << cShP->Data.RunNo <<" " << cShP->Data.EventNo << endl; getchar(); } } return 0; } // ::EntryFound() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Get NKG information for Curvature correction //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /*void MyReconstruction::GetNKGInfo() { float age(0),size(0); int fitflag(-10); shtree->GetEntry(RefSHIndex); //age=shtree->GetLeaf("Age")->GetValue(); //size=shtree->GetLeaf("NKGSize")->GetValue(); //fitflag=shtree->GetLeaf("NKGFitFlag")->GetValue(); age = cShP->Data.Age; size = cShP->Data.NKGSize; fitflag = cShP->Data.NKGFitFlag; sSD.NKGFitFlag=fitflag; sSD.Age_previous = age; sSDEven.Age_previous = age; sSDOdd.Age_previous = age; sSDLeft.Age_previous = age; sSDRight.Age_previous = age; sSD.Size_previous = size; sSDOdd.Size_previous = size; sSDEven.Size_previous = size; sSDLeft.Size_previous = size; sSDRight.Size_previous = size; } */ // ::GetNKGInfo() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Create output shYYYYMMDD.root files and corresponding trees //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::SetOutputs() { printf("\n"); if(inputs.dictionary.RECONSTRUCTION_RECONSTRUCT || inputs.dictionary.RECONSTRUCTION_CURVATURE_CORRECTION) { if(inputs.SplitIndexFlag) { sprintf(ReadName, "%s/sh%d_%06d.root", inputs.Output, sDTH.iDate, inputs.SplitIndexMin); } else { sprintf(ReadName, "%s/sh%d.root", inputs.Output, sDTH.iDate); } shRootFile = new TFile(ReadName, "RECREATE", "Reconstructed Shower Parameters"); if(shRootFile->IsOpen()) { printf(COLOR_OUTPUT "\n %s [Output]" COLOR_RESET, ReadName); CreateShParBranches(); } // if(shRootFile->IsOpen()) else { printf(COLOR_ERROR "\n IN %s::%s(%d)" COLOR_RESET, __FILE__, __func__, __LINE__); printf(COLOR_ERROR "\n **** ERROR : OPENING FILE %s ****\n\n" COLOR_RESET, ReadName); } } // RECONSTRUCTION_RECONSTRUCT if(inputs.dictionary.RECONSTRUCTION_GENERATE_PARTICLE_TOF) { sprintf(ReadName, "%s/pt%d.root", inputs.Output, sDTH.iDate); ptRootFile = new TFile(ReadName, "RECREATE", "Density & TOF Parameters"); if(ptRootFile->IsOpen()) { printf(COLOR_OUTPUT "\n %s [Output]" COLOR_RESET, ReadName); CreatePTBranches(); } else { printf(COLOR_ERROR "\n IN %s::%s(%d)" COLOR_RESET, __FILE__, __func__, __LINE__); printf(COLOR_ERROR "\n **** ERROR : OPENING FILE %s ****\n\n" COLOR_RESET, ReadName); } } printf("\n"); } // ::SetOutputs() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Create Shower Parameters branches //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::CreateShParBranches() { sh = new TTree("sh", "Shower Parameters Tree"); sh->Branch("RunNo", &sSP.RunNo, "RunNo/I"); sh->Branch("EventNo", &sSP.EventNo, "EventNo/I"); sh->Branch("EvDate", &sSP.EvDate, "EvDate/I"); sh->Branch("EvTime1", &sSP.EvTime1, "EvTime1/I"); sh->Branch("EvTime2", &sSP.EvTime2, "EvTime2/I"); sh->Branch("EvTimeFlag", &sSP.EvTimeFlag, "EvTimeFlag/I"); sh->Branch("nHitDetectors", &sSP.nHitDetectors, "nHitDetectors/I"); sh->Branch("PSum", &sSP.PSum, "PSum/D"); sh->Branch("PSumIn", &sSP.PSumIn, "PSumIn/D"); sh->Branch("PSumOut", &sSP.PSumOut, "PSumOut/D"); sh->Branch("PSumRatio", &sSP.PSumRatio, "PSumRatio/D"); sh->Branch("COGX", &sSP.COGX, "COGX/D"); sh->Branch("COGY", &sSP.COGY, "COGY/D"); sh->Branch("NKGX", &sSP.NKGX, "NKGX/D"); sh->Branch("NKGY", &sSP.NKGY, "NKGY/D"); sh->Branch("NKGZ", &sSP.NKGZ, "NKGZ/D"); sh->Branch("NKGXerr", &sSP.NKGXerr, "NKGXerr/D"); sh->Branch("NKGYerr", &sSP.NKGYerr, "NKGYerr/D"); sh->Branch("NKGZerr", &sSP.NKGZerr, "NKGZerr/D"); sh->Branch("LnCerr", &sSP.LnCerr, "LnCerr/D"); sh->Branch("Ageerr", &sSP.Ageerr, "Ageerr/D"); sh->Branch("NKGFitFlag", &sSP.NKGFitFlag, "NKGFitFlag/I"); sh->Branch("NKGStatus", &sSP.NKGStatus, "NKGStatus/I"); sh->Branch("NKGLoop", &sSP.NKGLoop, "NKGLoop/I"); sh->Branch("NKGSize", &sSP.NKGSize, "NKGSize/D"); sh->Branch("Age", &sSP.Age, "Age/D"); sh->Branch("Size", &sSP.Size, "Size/D"); sh->Branch("Energy", &sSP.Energy, "Energy/D"); sh->Branch("LnNKGP", &sSP.LnNKGP, "LnNKGP/D"); sh->Branch("EDM", &sSP.EDM, "EDM/D"); sh->Branch("Theta1", &sSP.Theta1, "Theta1/D"); sh->Branch("Phi1", &sSP.Phi1, "Phi1/D"); sh->Branch("ChiSq1", &sSP.ChiSq1, "ChiSq1/D"); sh->Branch("TZero1", &sSP.TZero1, "TZero1/D"); //sh->Branch("Theta2", &sSP.Theta2, "Theta2/D"); //sh->Branch("Phi2", &sSP.Phi2, "Phi2/D"); sh->Branch("Theta3", &sSP.Theta3, "Theta3/D"); sh->Branch("Phi3", &sSP.Phi3, "Phi3/D"); sh->Branch("Theta4", &sSP.Theta4, "Theta4/D"); sh->Branch("Phi4", &sSP.Phi4, "Phi4/D"); // Odd method sh->Branch("nHitDetectorsOdd", &sSP.nHitDetectorsOdd, "nHitDetectorsOdd/I"); sh->Branch("ThetaOdd1", &sSP.ThetaOdd1, "ThetaOdd1/D"); sh->Branch("PhiOdd1", &sSP.PhiOdd1, "PhiOdd1/D"); sh->Branch("ChiSqOdd1", &sSP.ChiSqOdd1, "ChiSqOdd1/D"); sh->Branch("TZeroOdd1", &sSP.TZeroOdd1, "TZeroOdd1/D"); //sh->Branch("ThetaOdd2", &sSP.ThetaOdd2, "ThetaOdd2/D"); //sh->Branch("PhiOdd2", &sSP.PhiOdd2, "PhiOdd2/D"); sh->Branch("ThetaOdd3", &sSP.ThetaOdd3, "ThetaOdd3/D"); sh->Branch("PhiOdd3", &sSP.PhiOdd3, "PhiOdd3/D"); sh->Branch("ThetaOdd4", &sSP.ThetaOdd4, "ThetaOdd4/D"); sh->Branch("PhiOdd4", &sSP.PhiOdd4, "PhiOdd4/D"); // Even method sh->Branch("nHitDetectorsEven", &sSP.nHitDetectorsEven, "nHitDetectorsEven/I"); sh->Branch("ThetaEven1", &sSP.ThetaEven1, "ThetaEven1/D"); sh->Branch("PhiEven1", &sSP.PhiEven1, "PhiEven1/D"); sh->Branch("ChiSqEven1", &sSP.ChiSqEven1, "ChiSqEven1/D"); sh->Branch("TZeroEven1", &sSP.TZeroEven1, "TZeroEven1/D"); //sh->Branch("ThetaEven2", &sSP.ThetaEven2, "ThetaEven2/D"); //sh->Branch("PhiEven2", &sSP.PhiEven2, "PhiEven2/D"); sh->Branch("ThetaEven3", &sSP.ThetaEven3, "ThetaEven3/D"); sh->Branch("PhiEven3", &sSP.PhiEven3, "PhiEven3/D"); sh->Branch("ThetaEven4", &sSP.ThetaEven4, "ThetaEven4/D"); sh->Branch("PhiEven4", &sSP.PhiEven4, "PhiEven4/D"); // Left method sh->Branch("nHitDetectorsLeft", &sSP.nHitDetectorsLeft, "nHitDetectorsLeft/I"); sh->Branch("ThetaLeft1", &sSP.ThetaLeft1, "ThetaLeft1/D"); sh->Branch("PhiLeft1", &sSP.PhiLeft1, "PhiLeft1/D"); sh->Branch("ChiSqLeft1", &sSP.ChiSqLeft1, "ChiSqLeft1/D"); sh->Branch("TZeroLeft1", &sSP.TZeroLeft1, "TZeroLeft1/D"); //sh->Branch("ThetaLeft2", &sSP.ThetaLeft2, "ThetaLeft2/D"); //sh->Branch("PhiLeft2", &sSP.PhiLeft2, "PhiLeft2/D"); sh->Branch("ThetaLeft3", &sSP.ThetaLeft3, "ThetaLeft3/D"); sh->Branch("PhiLeft3", &sSP.PhiLeft3, "PhiLeft3/D"); sh->Branch("ThetaLeft4", &sSP.ThetaLeft4, "ThetaLeft4/D"); sh->Branch("PhiLeft4", &sSP.PhiLeft4, "PhiLeft4/D"); // Right method sh->Branch("nHitDetectorsRight", &sSP.nHitDetectorsRight, "nHitDetectorsRight/I"); sh->Branch("ThetaRight1", &sSP.ThetaRight1, "ThetaRight1/D"); sh->Branch("PhiRight1", &sSP.PhiRight1, "PhiRight1/D"); sh->Branch("ChiSqRight1", &sSP.ChiSqRight1, "ChiSqRight1/D"); sh->Branch("TZeroRight1", &sSP.TZeroRight1, "TZeroRight1/D"); //sh->Branch("ThetaRight2", &sSP.ThetaRight2, "ThetaRight2/D"); //sh->Branch("PhiRight2", &sSP.PhiRight2, "PhiRight2/D"); sh->Branch("ThetaRight3", &sSP.ThetaRight3, "ThetaRight3/D"); sh->Branch("PhiRight3", &sSP.PhiRight3, "PhiRight3/D"); sh->Branch("ThetaRight4", &sSP.ThetaRight4, "ThetaRight4/D"); sh->Branch("PhiRight4", &sSP.PhiRight4, "PhiRight4/D"); if(inputs.dictionary.RECONSTRUCTION_CURVATURE_CORRECTION) { sh->Branch("Theta6", &sSP.Theta6, "Theta6/D"); sh->Branch("Phi6", &sSP.Phi6, "Phi6/D"); sh->Branch("ThetaEven6", &sSP.ThetaEven6, "ThetaEven6/D"); sh->Branch("PhiEven6", &sSP.PhiEven6, "PhiEven6/D"); sh->Branch("ThetaOdd6", &sSP.ThetaOdd6, "ThetOdd6/D"); sh->Branch("PhiOdd6", &sSP.PhiOdd6, "PhiOdd6/D"); sh->Branch("ThetaLeft6", &sSP.ThetaLeft6, "ThetaLeftt6/D"); sh->Branch("PhiLeft6", &sSP.PhiLeft6, "PhiLeft6/D"); sh->Branch("ThetaRight6", &sSP.ThetaRight6, "ThetaRight6/D"); sh->Branch("PhiRight6", &sSP.PhiRight6, "PhiRight6/D"); } } // ::CreateShParBranches() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Filling reconstructed parameters tree //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::FillShTree() { if(inputs.dictionary.RECONSTRUCTION_RECONSTRUCT) { sSP.RunNo = sEVTH.runno; sSP.EventNo = sEVTH.eventno; sSP.Trigger = sEVTH.trigger; sSP.EvDate = sEVTH.evdate; sSP.EvTime1 = sEVTH.evtime1; sSP.EvTime2 = sEVTH.evtime2; sSP.EvTimeFlag = p_ScTimeFlag; sSP.nHitDetectors = sSD.nHitDetectors; sSP.PSum = sSD.ParticleSum; sSP.PSumIn = PSumIn; sSP.PSumOut = PSumOut; sSP.PSumRatio = PSumOut/PSumIn; sSP.COGX = sSD.COGX; sSP.COGY = sSD.COGY; sSP.NKGStatus = sSD.NKGStatus; sSP.NKGFitFlag = sSD.NKGFitFlag; sSP.NKGLoop = sSD.NKGLoop; sSP.NKGX = sSD.NKGX; sSP.NKGY = sSD.NKGY; sSP.NKGZ = sSD.NKGZ; sSP.NKGXerr = sSD.NKGXerr; sSP.NKGYerr = sSD.NKGYerr; sSP.NKGZerr = sSD.NKGZerr; sSP.LnCerr = sSD.LnCerr; sSP.Ageerr = sSD.Ageerr; sSP.LnNKGP = sSD.LnNKGP; sSP.EDM = sSD.EDM; sSP.Energy = sSD.Energy; sSP.Size = sSD.Size ; sSP.NKGSize = sSD.NKGSize ; sSP.Age = sSD.Age; sSP.Theta1 = sSD.Angle1.Theta; sSP.Phi1 = sSD.Angle1.Phi; sSP.StdDev1 = sSD.Angle1.StdDev; sSP.ChiSq1 = sSD.Angle1.ChiSq; sSP.TZero1 = sSD.Angle1.TZero; sSP.Theta2 = sSD.Angle2.Theta; sSP.Phi2 = sSD.Angle2.Phi; sSP.StdDev2 = sSD.Angle2.StdDev; sSP.Theta3 = sSD.Angle3.Theta; sSP.Phi3 = sSD.Angle3.Phi; sSP.StdDev3 = sSD.Angle3.StdDev; sSP.Theta4 = sSD.Angle4.Theta; sSP.Phi4 = sSD.Angle4.Phi; sSP.StdDev4 = sSD.Angle4.StdDev; // Odd method sSP.nHitDetectorsOdd = sSDOdd.nHitDetectors; sSP.ThetaOdd1 = sSDOdd.Angle1.Theta; sSP.PhiOdd1 = sSDOdd.Angle1.Phi; sSP.ChiSqOdd1 = sSDOdd.Angle1.ChiSq; sSP.TZeroOdd1 = sSDOdd.Angle1.TZero; sSP.ThetaOdd2 = sSDOdd.Angle2.Theta; sSP.PhiOdd2 = sSDOdd.Angle2.Phi; sSP.ThetaOdd3 = sSDOdd.Angle3.Theta; sSP.PhiOdd3 = sSDOdd.Angle3.Phi; sSP.ThetaOdd4 = sSDOdd.Angle4.Theta; sSP.PhiOdd4 = sSDOdd.Angle4.Phi; // Even method sSP.nHitDetectorsEven = sSDEven.nHitDetectors; sSP.ThetaEven1 = sSDEven.Angle1.Theta; sSP.PhiEven1 = sSDEven.Angle1.Phi; sSP.ChiSqEven1 = sSDEven.Angle1.ChiSq; sSP.TZeroEven1 = sSDEven.Angle1.TZero; sSP.ThetaEven2 = sSDEven.Angle2.Theta; sSP.PhiEven2 = sSDEven.Angle2.Phi; sSP.ThetaEven3 = sSDEven.Angle3.Theta; sSP.PhiEven3 = sSDEven.Angle3.Phi; sSP.ThetaEven4 = sSDEven.Angle4.Theta; sSP.PhiEven4 = sSDEven.Angle4.Phi; // Left method sSP.nHitDetectorsLeft = sSDLeft.nHitDetectors; sSP.ThetaLeft1 = sSDLeft.Angle1.Theta; sSP.PhiLeft1 = sSDLeft.Angle1.Phi; sSP.ChiSqLeft1 = sSDLeft.Angle1.ChiSq; sSP.TZeroLeft1 = sSDLeft.Angle1.TZero; sSP.ThetaLeft2 = sSDLeft.Angle2.Theta; sSP.PhiLeft2 = sSDLeft.Angle2.Phi; sSP.ThetaLeft3 = sSDLeft.Angle3.Theta; sSP.PhiLeft3 = sSDLeft.Angle3.Phi; sSP.ThetaLeft4 = sSDLeft.Angle4.Theta; sSP.PhiLeft4 = sSDLeft.Angle4.Phi; // Right method sSP.nHitDetectorsRight = sSDRight.nHitDetectors; sSP.ThetaRight1 = sSDRight.Angle1.Theta; sSP.PhiRight1 = sSDRight.Angle1.Phi; sSP.ChiSqRight1 = sSDRight.Angle1.ChiSq; sSP.TZeroRight1 = sSDRight.Angle1.TZero; sSP.ThetaRight2 = sSDRight.Angle2.Theta; sSP.PhiRight2 = sSDRight.Angle2.Phi; sSP.ThetaRight3 = sSDRight.Angle3.Theta; sSP.PhiRight3 = sSDRight.Angle3.Phi; sSP.ThetaRight4 = sSDRight.Angle4.Theta; sSP.PhiRight4 = sSDRight.Angle4.Phi; } // RECONSTRUCTION_RECONSTRUCT if(inputs.dictionary.RECONSTRUCTION_CURVATURE_CORRECTION) { sSP = cShP->Data; sSP.Theta6 = sSD.AngleCone.Theta; sSP.Phi6 = sSD.AngleCone.Phi; sSP.TZero6 = sSD.AngleCone.TZero; sSP.ChiSq6 = sSD.AngleCone.ChiSq; sSP.ThetaOdd6 = sSDOdd.AngleCone.Theta; sSP.PhiOdd6 = sSDOdd.AngleCone.Phi; sSP.TZeroOdd6 = sSDOdd.AngleCone.TZero; sSP.ChiSqOdd6 = sSDOdd.AngleCone.ChiSq; sSP.ThetaEven6 = sSDEven.AngleCone.Theta; sSP.PhiEven6 = sSDEven.AngleCone.Phi; sSP.TZeroEven6 = sSDEven.AngleCone.TZero; sSP.ChiSqEven6 = sSDEven.AngleCone.ChiSq; sSP.ThetaRight6 = sSDRight.AngleCone.Theta; sSP.PhiRight6 = sSDRight.AngleCone.Phi; sSP.TZeroRight6 = sSDRight.AngleCone.TZero; sSP.ChiSqRight6 = sSDRight.AngleCone.ChiSq; sSP.ThetaLeft6 = sSDLeft.AngleCone.Theta; sSP.PhiLeft6 = sSDLeft.AngleCone.Phi; sSP.TZeroLeft6 = sSDLeft.AngleCone.TZero; sSP.ChiSqLeft6 = sSDLeft.AngleCone.ChiSq; sSP.Theta6 = sSD.Angle6[2].Theta; sSP.Phi6 = sSD.Angle6[2].Phi; sSP.TZero6 = sSD.Angle6[2].TZero; sSP.ChiSq6 = sSD.Angle6[2].ChiSq; sSP.ThetaOdd6 = sSDOdd.Angle6[2].Theta; sSP.PhiOdd6 = sSDOdd.Angle6[2].Phi; sSP.TZeroOdd6 = sSDOdd.Angle6[2].TZero; sSP.ChiSqOdd6 = sSDOdd.Angle6[2].ChiSq; sSP.ThetaEven6 = sSDEven.Angle6[2].Theta; sSP.PhiEven6 = sSDEven.Angle6[2].Phi; sSP.TZeroEven6 = sSDEven.Angle6[2].TZero; sSP.ChiSqEven6 = sSDEven.Angle6[2].ChiSq; sSP.ThetaRight6 = sSDRight.Angle6[2].Theta; sSP.PhiRight6 = sSDRight.Angle6[2].Phi; sSP.TZeroRight6 = sSDRight.Angle6[2].TZero; sSP.ChiSqRight6 = sSDRight.Angle6[2].ChiSq; sSP.ThetaLeft6 = sSDLeft.Angle6[2].Theta; sSP.PhiLeft6 = sSDLeft.Angle6[2].Phi; sSP.TZeroLeft6 = sSDLeft.Angle6[2].TZero; sSP.ChiSqLeft6 = sSDLeft.Angle6[2].ChiSq; } // RECONSTRUCTION_CURVATURE_CORRECTION sh->Fill(); } // ::FillShTree() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Create particle density and time of flight branches //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::CreatePTBranches() { pt = new TTree("pt", "Density & TOF Parameters Tree"); sPT.nHitDetectors = 0; pt->Branch("RunNo", &sPT.RunNo, "RunNo/I"); pt->Branch("EventNo", &sPT.EventNo, "EventNo/I"); pt->Branch("EvDate", &sPT.EvDate, "EvDate/I"); pt->Branch("EvTime1", &sPT.EvTime1, "EvTime1/I"); pt->Branch("EvTime2", &sPT.EvTime2, "EvTime2/I"); pt->Branch("nHitDetectors", &sPT.nHitDetectors, "nHitDetectors/I"); pt->Branch("DetNo", &sPT.DetNo, "DetNo[nHitDetectors]/S"); pt->Branch("Density", &sPT.Density, "Density[nHitDetectors]/F"); pt->Branch("TOF", &sPT.TOF, "TOF[nHitDetectors]/F"); } // ::CreatePTBranches() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Filling shower parameters tree //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::FillPTTree() { sPT.RunNo = sEVTH.runno; sPT.EventNo = sEVTH.eventno; sPT.EvDate = sEVTH.evdate; sPT.EvTime1 = sEVTH.evtime1; sPT.EvTime2 = sEVTH.evtime2; pt->Fill(); } // ::FillPTTree() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Get ScHeader branches //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::GetScHeaderBranches() { // ScHeader tree schtree->SetBranchAddress("RunNo", &p_ScRunNo); schtree->SetBranchAddress("EventNo", &p_ScEvNo); schtree->SetBranchAddress("EvDate", &p_ScEvDate); schtree->SetBranchAddress("EvTime1", &p_ScEvTime1); schtree->SetBranchAddress("EvTime2", &p_ScEvTime2); schtree->SetBranchAddress("TimeFlag", &p_ScTimeFlag); } // ::GetScHeaderBranches() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Reset shower parameters data //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::ResetPTData(struct DensityTOF *sdt) { sdt->RunNo = -1; sdt->EventNo = -1; sdt->EvDate = -1; sdt->EvTime1 = -1; sdt->EvTime2 = -1; sdt->nHitDetectors = 0; for(int iDet(0); iDet < MAXD; iDet++) { sdt->Density[iDet] = 0; sdt->TOF[iDet] = 1E10; } // for iDet < MAXD } // ::ResetPTData() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// // Update detector weights //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// inline void MyReconstruction::UpdateWeights(struct ShowerData *ssd, int flag) { if(inputs.dictionary.RECONSTRUCTION_CURVATURE_CORRECTION == 1) { ssd->Age_Reference = cShP->Data.Age; ssd->Size_Reference = cShP->Data.NKGSize; ssd->NKGFlag_Reference = cShP->Data.NKGFitFlag; if(cShP->Data.NKGFitFlag > 0) { ssd->NKGX_Reference = cShP->Data.NKGX; ssd->NKGY_Reference = cShP->Data.NKGY; } else { ssd->NKGX_Reference = cShP->Data.COGX; ssd->NKGY_Reference = cShP->Data.COGY; } ssd->Th1_Reference = cShP->Data.Theta1; ssd->Ph1_Reference = cShP->Data.Phi1; } for(int iDet(1); iDet < MAXD; iDet++) { ssd->ADCWeight[iDet] = 0; ssd->TDCWeight[iDet] = 0; ssd->X[iDet] = -1000; ssd->Y[iDet] = -1000; ssd->Z[iDet] = -1000; ssd->DetStatus[iDet] = -1; if(cGDB->dDBDI.status[iDet] == 1) { if(inputs.dictionary.RECONSTRUCTION_RECONSTRUCT_PMT == 1 || (inputs.dictionary.RECONSTRUCTION_RECONSTRUCT_PMT > 1 && cGDB->dDBDI.npmt[iDet] == 2) || inputs.dictionary.RECONSTRUCTION_CURVATURE_CORRECTION == 1) { switch(flag) { // Full array case 0: { ssd->ADCWeight[iDet] = 1; ssd->TDCWeight[iDet] = 1; ssd->X[iDet] = cGDB->dDBDI.detX[iDet]; ssd->Y[iDet] = cGDB->dDBDI.detY[iDet]; ssd->Z[iDet] = cGDB->dDBDI.detZ[iDet]; ssd->DetStatus[iDet] = 1; break; } // Odd array case 1: { if(iDet % 2 == 1) { ssd->ADCWeight[iDet] = 1; ssd->TDCWeight[iDet] = 1; ssd->X[iDet] = cGDB->dDBDI.detX[iDet]; ssd->Y[iDet] = cGDB->dDBDI.detY[iDet]; ssd->Z[iDet] = cGDB->dDBDI.detZ[iDet]; ssd->DetStatus[iDet] = 1; } break; } // Even array case 2: { if(iDet % 2 == 0) { ssd->ADCWeight[iDet] = 1; ssd->TDCWeight[iDet] = 1; ssd->X[iDet] = cGDB->dDBDI.detX[iDet]; ssd->Y[iDet] = cGDB->dDBDI.detY[iDet]; ssd->Z[iDet] = cGDB->dDBDI.detZ[iDet]; ssd->DetStatus[iDet] = 1; } break; } }; } } // if(status) } // for iDet < MAXD } // ::UpdateWeights() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Read detector coordinates //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::ReadDetectorCoordinates() { // Scintillator detector coordinates ScPos = ReadScCoordinates((char *)"misc/scdet_pos.txt"); // Copy this data to ShowerData struct for(unsigned short iDet(1); iDet < MAXD; iDet++) { sSD.X[iDet] = -1000; sSDOdd.X[iDet] = -1000; sSDEven.X[iDet] = -1000; sSD.Y[iDet] = -1000; sSDOdd.Y[iDet] = -1000; sSDEven.Y[iDet] = -1000; sSD.Z[iDet] = -1000; sSDOdd.Z[iDet] = -1000; sSDEven.Z[iDet] = -1000; if(ScPos.status[iDet] != 1) { continue; } sSD.X[iDet] = ScPos.x[iDet]; sSDOdd.X[iDet] = ScPos.x[iDet]; sSDEven.X[iDet] = ScPos.x[iDet]; sSD.Y[iDet] = ScPos.y[iDet]; sSDOdd.Y[iDet] = ScPos.y[iDet]; sSDEven.Y[iDet] = ScPos.y[iDet]; sSD.Z[iDet] = ScPos.z[iDet]; sSDOdd.Z[iDet] = ScPos.z[iDet]; sSDEven.Z[iDet] = ScPos.z[iDet]; } // for iDet < MAXD } // ::ReadDetectorCoordinates() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Print shower data of one event //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::PrintShowerData() { printf("\n nHitDetectors = %3d", sSD.nHitDetectors); printf("\n ParticleSum = %7.2f", sSD.ParticleSum); for(unsigned short iDet(1); iDet < MAXD; iDet++) { if(sSD.HitFlag[iDet] != 1) { continue; } printf("\n %3d %d %7.2f %7.2f %7.2f %7.2f %7.2f", iDet, sSD.HitFlag[iDet], sSD.X[iDet], sSD.Y[iDet], sSD.Z[iDet], sSD.T[iDet], sSD.Density[iDet]); } // for iDet < MAXD } // ::PrintShowerData() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Compare sc run entries and db entries. Both should be same /// RETURN STATUS /// 0 : FAILURE /// 1 : SUCCESS //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// bool MyReconstruction::ValidateEntries() { if(RunEntries != DBRunEntries) { printf(COLOR_WARNING "\n **** WARNING : RUN & DB ENTRIES ARE NOT SAME ****\n\n" COLOR_RESET); return 0; } if(EventEntries != ScHEntries) { printf(COLOR_WARNING "\n **** WARNING : SC (%d) & SCHEADER (%d) ENTRIES ARE NOT SAME ****\n\n" COLOR_RESET, EventEntries, ScHEntries); return 0; } for(int i(0); i < RunEntries; i++) { } // for i < RunEntries return 1; } // ::ValidateEntries() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Validate runno in both stream //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::ValidateRunno(int runno) { if(!(runno == sSRP.runno && runno == cGDB->dDBRP.runno)) { printf(COLOR_ERROR "\n IN %s::%s(%d)" COLOR_RESET, __FILE__, __func__, __LINE__); printf(COLOR_ERROR "\n **** ERROR : RUN NOT MATCHING IN BOTH STREAM %6d %6d %6d ****\n\n" COLOR_RESET, runno, sSRP.runno, cGDB->dDBRP.runno); } } // ::ValidateRunno() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Calculate array center using center of gravity method //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::CalculateArrayCenter() { ArrayCenterX = 0.0; ArrayCenterY = 0.0; Det = 0; //printf("\n cGDB->dDBPED.ndet = %3d", cGDB->dDBPED.ndet); //printf("\n cGDB->dDBDI.ndet = %3d", cGDB->dDBDI.ndet); for(int iDet(1); iDet < MAXD; iDet++) { if(cGDB->dDBDI.status[iDet] != 1) { continue; } if(cGDB->dDBPED.detno[iDet] == -1) { continue; } //printf("\n %04d %1d %10.1f %10.1f %10.1f", iDet, cGDB->dDBDI.status[iDet], cGDB->dDBDI.detX[iDet], cGDB->dDBDI.detY[iDet], cGDB->dDBDI.detZ[iDet]); Det++; ArrayCenterX += cGDB->dDBDI.detX[iDet]; ArrayCenterY += cGDB->dDBDI.detY[iDet]; } // for iDet < MAXD if(Det > 0) { ArrayCenterX /= Det; ArrayCenterY /= Det; } else { ArrayCenterX = 0.0; ArrayCenterY = 0.0; } printf("\n Det = %04d", Det); printf("\n Array Center X = %10.2f", ArrayCenterX); printf("\n Array Center Y = %10.2f", ArrayCenterY); getchar(); } // ::CalculateArrayCenter() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Reset PSum values //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::ResetPSum() { PSumIn = 0; PSumOut = 0; PSumRatio = 0; } // ::ResetPSum() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Process detectorwise //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// inline void MyReconstruction::ProcessDetectors(int flag) { for (int i = 0; i < sSC.ndet; ++i) { int Det = sSC.detno[i]; TString histName = Form("hist_%d", Det); if (hParticleDensity.find(Det) == hParticleDensity.end()) { hParticleDensity[Det] = new TH1F(histName, "Particle Density", 1010, 0.0, 1000.0); std::cout << "Created histogram: " << hParticleDensity[Det]->GetName() << std::endl; } // When PMT option is selected for LG or HGofLG then if npmt is not equal to 2 then ignore this detector if(inputs.dictionary.RECONSTRUCTION_RECONSTRUCT_PMT > 1 && cGDB->dDBDI.npmt[Det] != 2) { continue; } //printf("\n %3d %3d tdc0DQF %d adcgrDQF %d", i, Det, cGDB->dDBTZ.tdc0DQF[Det], cGDB->dDBGRHP.adcgrDQF[Det]); if(!ValidateDetector(Det)) { continue; } //printf("\n %3d %3d = ", i, Det); if(!ConvertToPhysicalQuantities(i, Det, flag)) { continue; } //hParticleDensity[Det]->Fill(density); //std::cout << "Filled histogram: " << hParticleDensity[Det]->GetName() << " with density: " << density << std::endl; // Get PSum ratio inside and outside the fiducial area if(flag == 1) { GetPSumRatio(Det); } //printf("\n %3d %3d ==", i, Det); //printf(" %4d %7.2f %7.2f %7.2f ", ADCHH, cGDB->dDBPED.pedhhpeak[Det], cGDB->dDBPED.pedhhmean[Det], fADCHH); //printf(" %4d %7.2f %7.2f %7.2f ", ADCHL, cGDB->dDBPED.pedhlpeak[Det], cGDB->dDBPED.pedhlmean[Det], fADCHL); //printf(" %7.2f %7.2f %7.2f", cGDB->dDBCG.computedadc[Det], cGDB->dDBGRHP.adcgr[Det], density); //printf("\n %4d %7.2f %7.2f %d", TDC, cGDB->dDBTZ.tdc0[Det], tof, cGDB->sDBTZ.tdc0DQF[Det]); if(inputs.dictionary.RECONSTRUCTION_RECONSTRUCT || inputs.dictionary.RECONSTRUCTION_CURVATURE_CORRECTION) { if(flag == 1) { ProcessFullArray(Det); // if(inputs.dictionary.RECONSTRUCTION_RECONSTRUCT_PMT == 1 && inputs.dictionary.RECONSTRUCTION_RECONSTRUCT_EVENODD) // { // ProcessEvenOddArray(Det); //} // EVENODD } else { //if(inputs.dictionary.RECONSTRUCTION_RECONSTRUCT_PMT == 1 && inputs.dictionary.RECONSTRUCTION_RECONSTRUCT_LEFTRIGHT) // { // If full array shower core is not estimated then ignore // if(sSD.COGX == 0.0 && sSD.COGY == 0.0) // { // continue; // } // ProcessLeftRightArray(Det); } // LEFTRIGHT // } // else } // RECONSTRUCTION_RECONSTRUCT } // for i < sSC.ndet //getchar(); } // ::ProcessDetectors() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Calculate the Particle sum inside and outside the fiducial area //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// void MyReconstruction::GetPSumRatio(int det) { if(density < inputs.dictionary.RECONSTRUCTION_PARTICLE_CUT) { return; } if(cGDB->dDBDI.status[det] != 1) { return; } // Check if the detector is inside or outside if(IsCoreInside(cGDB->dDBDI.detX[Det],cGDB->dDBDI.detY[Det])==0) { PSumOut += density; } else { PSumIn += density; } // PSumRatio = PSumOut/PSumIn; } // ::GetPSumRatio() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Process for full array //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// inline void MyReconstruction::ProcessFullArray(int det) { // All detectors method sSD.T[det] = tof; sSD.Density[det] = density; sSD.HitFlag[det] = 1; //sSD.ADCWeight[det] = 1; //sSD.TDCWeight[det] = 1; sSD.ParticleSum += density; sSD.nHitDetectors++; } // ::ProcessFullArray() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Process for Even-Odd array //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /* inline void MyReconstruction::ProcessEvenOddArray(int det) { // Even-Odd method if(det % 2 == 0) { sSDEven.T[det] = tof; sSDEven.Density[det] = density; sSDEven.HitFlag[det] = 1; //sSDEven.ADCWeight[det] = 1; //sSDEven.TDCWeight[det] = 1; sSDEven.ParticleSum += density; sSDEven.nHitDetectors++; } else { sSDOdd.T[det] = tof; sSDOdd.Density[det] = density; sSDOdd.HitFlag[det] = 1; //sSDOdd.ADCWeight[det] = 1; //sSDOdd.TDCWeight[det] = 1; sSDOdd.ParticleSum += density; sSDOdd.nHitDetectors++; } } */ // ::ProcessEvenOddArray() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Process for Left-Right array //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /* inline void MyReconstruction::ProcessLeftRightArray(int det) { printf("\n Distance = %f", WhichSideOfLine(sSD.COGX, sSD.COGY, ArrayCenterX, ArrayCenterY, sSD.X[det], sSD.Y[det])); getchar(); // Left side if(ReconstructionCore::WhichSideOfLine(sSD.COGX, sSD.COGY, ArrayCenterX, ArrayCenterY, sSD.X[det], sSD.Y[det]) <= 0) { sSDLeft.T[det] = tof; sSDLeft.Density[det] = density; sSDLeft.HitFlag[det] = 1; //sSDLeft.ADCWeight[det] = 1; //sSDLeft.TDCWeight[det] = 1; sSDLeft.ParticleSum += density; sSDLeft.nHitDetectors++; } // Right side else { sSDRight.T[det] = tof; sSDRight.Density[det] = density; sSDRight.HitFlag[det] = 1; //sSDRight.ADCWeight[det] = 1; //sSDRight.TDCWeight[det] = 1; sSDRight.ParticleSum += density; sSDRight.nHitDetectors++; } } */ // ::ProcessLeftRightArray() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Print reconstructed data //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /*void MyReconstruction::PrintReconstructedData() { printf("\n nHitDetectors = %3d %3d %3d", sSD.nHitDetectors, sSDOdd.nHitDetectors, sSDEven.nHitDetectors); printf("\n Hits = %3d", sSD.nHitDetectors); printf("\n Theta = %f\t%f\t%f\t%f\t%f", sSD.Angle1.Theta, sSD.Angle2.Theta, sSD.Angle3.Theta, sSD.Angle4.Theta, sSD.Angle5.Theta); printf("\n Phi = %f\t%f\t%f\t%f\t%f", sSD.Angle1.Phi, sSD.Angle2.Phi, sSD.Angle3.Phi, sSD.Angle4.Phi, sSD.Angle5.Phi); printf("\n Size = %f\t%f", sSD.Size, sSD.NKGSize); printf("\n Energy = %f", sSD.Energy); printf("\n CoreX = %f\t%f", sSD.COGX, sSD.NKGX); printf("\n CoreY = %f\t%f", sSD.COGY, sSD.NKGY); printf("\n Age = %f", sSD.Age); getchar(); } */ // ::PrintReconstructedData() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Process PT data here //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// inline void MyReconstruction::ProcessPTData(int det) { sPT.DetNo[sPT.nHitDetectors] = det; sPT.TOF[sPT.nHitDetectors] = tof; sPT.Density[sPT.nHitDetectors] = density; sPT.nHitDetectors++; } // ::ProcessPTData() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Convert to physical quantities /// RETURN STATUS /// 0 : INVALID /// 1 : VALID /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// inline bool MyReconstruction::ConvertToPhysicalQuantities(int index, int det, int flag) { std::cout << "Calling ConvertToPhysicalQuantities for detector " << det << std::endl; // Calculation of Time Of Flight TDC = sSC.tdc[index]; // If TDC information is invalid skip the detector if(TDC == -1) { //printf("\n %3d = %3d, %5d", index, det, TDC); return 0; } tof = (TDC / 5.119604); tof = (TDC / 5.119604) - cGDB->dDBTZ.tdc0[det]; //for (int det = 0; det < Det; ++det) { // Calulation of particle density ADCHH = sSC.adchh[index]; /* if(cGDB->dDBPED.pedhhfitflag[det] == 1) { fADCHH = ADCHH - cGDB->dDBPED.pedhhpeak[det]; } else { fADCHH = ADCHH - cGDB->dDBPED.pedhhmean[det]; } */ //fADCHH = ADCHH - cGDB->dDBPED.pedhhmean[det]; // Old pedestal data fADCHH = GetPedestalSubtractedADC(ADCHH, det, 11); //if(cGDB->dDBPED2.pedhhflag[det] == 1) //{ // cout << fADCHH << " " << ADCHH-cGDB->dDBPED2.pedhhmean[det] << " " << ADCHH-cGDB->dDBPED.pedhhmean[det] << endl; // getchar(); //} ADCHL = sSC.adchl[index]; /* if(cGDB->dDBPED.pedhlfitflag[det] == 1) { fADCHL = ADCHL - cGDB->dDBPED.pedhlpeak[det]; } else { fADCHL = ADCHL - cGDB->dDBPED.pedhlmean[det]; } */ //fADCHL = ADCHL - cGDB->dDBPED.pedhlmean[det]; // Old pedestal data fADCHL = GetPedestalSubtractedADC(ADCHL, det, 12); // Before digital saturation if(ADCHH < 3500) { density = fADCHH / cGDB->dDBCG.computedadc[det]; // For Hourly Gain //density = fADCHH / cGDB->dDBMC.muadcpeak[det]; // For Calibration Gain } // Digital saturation else { density = fADCHL * cGDB->dDBGRHP.adcgr[det] / cGDB->dDBCG.computedadc[det]; //density = fADCHL * cGDB->dDBGRHP.adcgr[det] / cGDB->dDBMC.muadcpeak[det]; } // Only low-gain PMTs if(inputs.dictionary.RECONSTRUCTION_RECONSTRUCT_PMT == 2 && density > 50.0) { ADCLH = sSC.adclh[index]; /* if(cGDB->dDBPED.pedlhfitflag[det] == 1) { fADCLH = ADCLH - cGDB->dDBPED.pedlhpeak[det]; } else { fADCLH = ADCLH - cGDB->dDBPED.pedlhmean[det]; } */ //fADCLH = ADCLH - cGDB->dDBPED.pedlhmean[det]; fADCLH = GetPedestalSubtractedADC(ADCLH, det, 21); ADCLL = sSC.adcll[index]; /* if(cGDB->dDBPED.pedllfitflag[det] == 1) { fADCLL = ADCLL - cGDB->dDBPED.pedllpeak[det]; } else { fADCLL = ADCLL - cGDB->dDBPED.pedllmean[det]; } */ //fADCLL = ADCLL - cGDB->dDBPED.pedllmean[det]; fADCLL = GetPedestalSubtractedADC(ADCLL, det, 22); // Before digital saturation if(ADCLH < 3500) { density = fADCLH / (cGDB->dDBCG.computedadc[det] / cGDB->dDBLGR.ratio[det]); //density = fADCLH / (cGDB->dDBMC.muadcpeak[det] / cGDB->dDBLGR.ratio[det]); } // Digital saturation else { density = fADCLL * cGDB->dDBGRLP.adcgr[det] / (cGDB->dDBCG.computedadc[det] / cGDB->dDBLGR.ratio[det]); //density = fADCLL * cGDB->dDBGRLP.adcgr[det] / (cGDB->dDBMC.muadcpeak[det] / cGDB->dDBLGR.ratio[det]); } } if(inputs.dictionary.RECONSTRUCTION_GENERATE_PARTICLE_TOF && flag == 1) { ProcessPTData(det); } // RECONSTRUCTION_GENERATE_PARTICLE_TOF //printf(", Density = %10.3f, TOF = %10.3f", density, tof); // Reject signals beyond certain bandwidth (300.0 nS <= TOF <= 1200.0) if(tof < 300.0 || tof > 1200.0) { return 0; } // Reject low energy noise if(density < inputs.dictionary.RECONSTRUCTION_PARTICLE_CUT || tof > MAXTOF) { //printf("\n %3d = %3d, %7.2f", index, Det, density); if (hParticleDensity.find(det) != hParticleDensity.end()) { std::cout << "Filling histogram for detector " << det << " with density: " << density << std::endl; hParticleDensity[det]->Fill(density); } else { std::cout << "Histogram for detector " << det << " not found in the map" << std::endl; } return 0; } return 1; } // ::ConvertToPhysicalQuantities() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Substract the pedestal and return the pedestal substracted ADC /// /// dDBPED : Use this structure to use old pedestal values (using fitting technique) /// dDBPED2: Use this structure to use the new Pedestal values (Running Average technique) /// dDBPEDAvg: When problematic pedestal appears //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// double MyReconstruction::GetPedestalSubtractedADC(double adc, int det, int flag) { double pedADC = -1; // Low gains are forced to use the mean of the running average (until properly arranged) cGDB->dDBPED2.pedlhflag[det] = 0; cGDB->dDBPED2.pedllflag[det] = 0; switch(flag) { case 11: if(cGDB->dDBPED2.pedhhflag[det] == 0) { pedADC = adc - cGDB->dDBPED2.pedhhmean[det]; // New pedestal data } else if(cGDB->dDBPED2.pedhhflag[det] == 1) { for(int iped(0); iped < cGDB->dDBPEDAvg.nbadpedhh; iped++) { if(det == cGDB->dDBPEDAvg.hhbaddet[iped]) { if(CurrentEntry < cGDB->dDBPEDAvg.hhpedentry[iped]) { if(iped == 0 || det != cGDB->dDBPEDAvg.hhbaddet[iped-1]) { return adc - cGDB->dDBPEDAvg.hhpedestal[iped]; // New pedestal data //cout << cGDB->dDBPEDAvg.hhbaddet[iped] << " " << CurrentEntry << " "<< cGDB->dDBPEDAvg.hhpedentry[iped] << endl; } else if(CurrentEntry > cGDB->dDBPEDAvg.hhpedentry[iped-1]) { return adc - cGDB->dDBPEDAvg.hhpedestal[iped]; // New pedestal data //cout << cGDB->dDBPEDAvg.hhbaddet[iped] << " "<< cGDB->dDBPEDAvg.hhpedentry[iped-1] <<" '" << CurrentEntry << "' "<< cGDB->dDBPEDAvg.hhpedentry[iped]<< endl; } } // if( CurrentEntry < PedEntry) // When the events are beyond last pedestal trigger but in the same run if(det != cGDB->dDBPEDAvg.hhbaddet[iped+1] && CurrentEntry > cGDB->dDBPEDAvg.hhpedentry[iped]) { //cout << cGDB->dDBPEDAvg.hhbaddet[iped] << " "<< cGDB->dDBPEDAvg.hhpedentry[iped-1] <<" '" << CurrentEntry << "' "<< cGDB->dDBPEDAvg.hhpedentry[iped]<< endl; return adc - cGDB->dDBPEDAvg.hhpedestal[iped]; // New pedestal data } } // if(det == pedet) } // for(iPed) } // if(flag != 0) break; case 12: if(cGDB->dDBPED2.pedhlflag[det] == 0) { pedADC = adc - cGDB->dDBPED2.pedhlmean[det]; // New pedestal data } else if(cGDB->dDBPED2.pedhlflag[det] == 1) { for(int iped(0); iped < cGDB->dDBPEDAvg.nbadpedhl; iped++) { if(det == cGDB->dDBPEDAvg.hlbaddet[iped]) { if(CurrentEntry < cGDB->dDBPEDAvg.hlpedentry[iped]) { if(iped == 0 || det != cGDB->dDBPEDAvg.hlbaddet[iped-1]) { return adc - cGDB->dDBPEDAvg.hlpedestal[iped]; // New pedestal data //cout << cGDB->dDBPEDAvg.hlbaddet[iped] << " " << CurrentEntry << " "<< cGDB->dDBPEDAvg.hlpedentry[iped] << endl; } else if(CurrentEntry > cGDB->dDBPEDAvg.hlpedentry[iped-1]) { return adc - cGDB->dDBPEDAvg.hlpedestal[iped]; // New pedestal data //cout << cGDB->dDBPEDAvg.hlbaddet[iped] << " "<< cGDB->dDBPEDAvg.hlpedentry[iped-1] <<" '" << CurrentEntry << "' "<< cGDB->dDBPEDAvg.hlpedentry[iped]<< endl; } } // if( CurrentEntry < PedEntry) // When the events are beyond last pedestal trigger but in the same run if(det != cGDB->dDBPEDAvg.hlbaddet[iped+1] && CurrentEntry > cGDB->dDBPEDAvg.hlpedentry[iped]) { //cout << cGDB->dDBPEDAvg.hlbaddet[iped] << " "<< cGDB->dDBPEDAvg.hlpedentry[iped-1] <<" '" << CurrentEntry << "' "<< cGDB->dDBPEDAvg.hlpedentry[iped]<< endl; return adc - cGDB->dDBPEDAvg.hlpedestal[iped]; // New pedestal data } } // if(det == pedet) } // for(iPed) } // if(flag != 0) break; case 21: if(cGDB->dDBPED2.pedlhflag[det] == 0) { pedADC = adc - cGDB->dDBPED2.pedlhmean[det]; // New pedestal data } else if(cGDB->dDBPED2.pedlhflag[det] == 1) { for(int iped(0); iped < cGDB->dDBPEDAvg.nbadpedlh; iped++) { if(det == cGDB->dDBPEDAvg.lhbaddet[iped]) { if(CurrentEntry < cGDB->dDBPEDAvg.lhpedentry[iped]) { if(iped == 0 || det != cGDB->dDBPEDAvg.lhbaddet[iped-1]) { return adc - cGDB->dDBPEDAvg.lhpedestal[iped]; // New pedestal data //cout << cGDB->dDBPEDAvg.lhbaddet[iped] << " " << CurrentEntry << " "<< cGDB->dDBPEDAvg.lhpedentry[iped] << endl; } else if(CurrentEntry > cGDB->dDBPEDAvg.lhpedentry[iped-1]) { return adc - cGDB->dDBPEDAvg.lhpedestal[iped]; // New pedestal data //cout << cGDB->dDBPEDAvg.lhbaddet[iped] << " "<< cGDB->dDBPEDAvg.lhpedentry[iped-1] <<" '" << CurrentEntry << "' "<< cGDB->dDBPEDAvg.lhpedentry[iped]<< endl; } } // if( CurrentEntry < PedEntry) // When the events are beyond last pedestal trigger but in the same run if(det != cGDB->dDBPEDAvg.lhbaddet[iped+1] && CurrentEntry > cGDB->dDBPEDAvg.lhpedentry[iped]) { //cout << cGDB->dDBPEDAvg.lhbaddet[iped] << " "<< cGDB->dDBPEDAvg.lhpedentry[iped-1] <<" '" << CurrentEntry << "' "<< cGDB->dDBPEDAvg.lhpedentry[iped]<< endl; return adc - cGDB->dDBPEDAvg.lhpedestal[iped]; // New pedestal data } } // if(det == pedet) } // for(iPed) } // if(flag != 0) break; case 22: if(cGDB->dDBPED2.pedllflag[det] == 0) { pedADC = adc - cGDB->dDBPED2.pedllmean[det]; // New pedestal data } else if(cGDB->dDBPED2.pedllflag[det] == 1) { for(int iped(0); iped < cGDB->dDBPEDAvg.nbadpedll; iped++) { if(det == cGDB->dDBPEDAvg.llbaddet[iped]) { if(CurrentEntry < cGDB->dDBPEDAvg.llpedentry[iped]) { if(iped == 0 || det != cGDB->dDBPEDAvg.llbaddet[iped-1]) { return adc - cGDB->dDBPEDAvg.llpedestal[iped]; // New pedestal data //cout << cGDB->dDBPEDAvg.llbaddet[iped] << " " << CurrentEntry << " "<< cGDB->dDBPEDAvg.llpedentry[iped] << endl; } else if(CurrentEntry > cGDB->dDBPEDAvg.llpedentry[iped-1]) { return adc - cGDB->dDBPEDAvg.llpedestal[iped]; // New pedestal data //cout << cGDB->dDBPEDAvg.llbaddet[iped] << " "<< cGDB->dDBPEDAvg.llpedentry[iped-1] <<" '" << CurrentEntry << "' "<< cGDB->dDBPEDAvg.llpedentry[iped]<< endl; } } // if( CurrentEntry < PedEntry) // When the events are beyond last pedestal trigger but in the same run if(det != cGDB->dDBPEDAvg.llbaddet[iped+1] && CurrentEntry > cGDB->dDBPEDAvg.llpedentry[iped]) { //cout << cGDB->dDBPEDAvg.llbaddet[iped] << " "<< cGDB->dDBPEDAvg.llpedentry[iped-1] <<" '" << CurrentEntry << "' "<< cGDB->dDBPEDAvg.llpedentry[iped]<< endl; return adc - cGDB->dDBPEDAvg.llpedestal[iped]; // New pedestal data } } // if(det == pedet) } // for(iPed) } // if(flag != 0) break; }; return pedADC; } // ::GetPedestalSubtractedADC() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// /// Validate detector /// RETURN STATUS /// 0 : INVALID /// 1 : VALID //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// inline bool MyReconstruction::ValidateDetector(int det) { // Skip calibration detectors if(det > 800) { //printf("\n %3d = %3d", i, Det); return 0; } // Validate the data quality flags // 1) If detector coordinates are available then flag 1 // 2) If gain is not available then flag -1 (1 : within 2%, 2 : within 5%, 3 : within 10% of calibration gain) // 3) If gain ratio is not available then flag -1 // 4) If TDC zero is not available then flag -1 if(cGDB->dDBDI.status[det] != 1) { return 0; } /* if(cGDB->dDBPED.pedhhfitflag[det] == 1) { //if(cGDB->dDBPED.pedhhsigma[det] > 10) if(cGDB->dDBPED.pedhhsigma[det] > 20) // temproray { return 0; } } else { //if(cGDB->dDBPED.pedhhrms[det] > 10) if(cGDB->dDBPED.pedhhrms[det] > 20) // temporary { return 0; } } //if(cGDB->dDBPED.pedhhrms[det] > 10) */ // If new pedestal parameters are not proper if(cGDB->dDBPED2.pedhhflag[det] == -1) { return 0; } if(cGDB->dDBPED2.pedhhmean[det] == -1) { return 0; } // RMS of new pedestal is high if(cGDB->dDBPED2.pedhhrms[det] > 10) { return 0; } // If computed gain is not proper if(cGDB->dDBCG.adcDQF[det] <= 0) { return 0; } // If Gain ratio is improper if(cGDB->dDBGRHP.adcgrDQF[det] != 1) { return 0; } // Id TDC0 is not present if(cGDB->dDBTZ.tdc0DQF[det] != 1) { return 0; } /* if(cGDB->dDBDQF.ADCFlag[det] != 1) { return 0; } */ if(cGDB->dDBDQF.TDCFlag[det] != 1) { return 0; } // For low-gain PMTs if(inputs.dictionary.RECONSTRUCTION_RECONSTRUCT_PMT == 2) { if(cGDB->dDBPED.pedlhfitflag[det] == 1) { if(cGDB->dDBPED.pedlhsigma[det] > 10) { //return 0; } } else { if(cGDB->dDBPED.pedlhrms[det] > 10) { //return 0; } } if(cGDB->dDBLGR.lgDQF[det] <= 0) { return 0; } if(cGDB->dDBGRLP.adcgrDQF[det] != 1) { return 0; } } // RECONSTRUCTION_RECONSTRUCT_PMT return 1; } // ::ValidateDetector() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ///To check if coordinates are inside the fiducial area, added by Medha //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// bool MyReconstruction::IsCoreInside(float x, float y) { float p_x1 = -91.5; float p_y1 = 16.0; float p_x2 = -64.0; float p_y2 = 64.5; float p_x3 = 40.8; float p_y3 = 64.5; float p_x4 = 64.0; float p_y4 = 23.5; float p_x5 = 16.0; float p_y5 = -59.0; float p_x6 = -51.0; float p_y6 = -59.0; if(x < p_x1 || x > p_x4 || y < p_y6 || y > p_y2) { return 0; } // region 2 if(y > 0) { if(x > p_x2 && x < p_x3) { return 1; } if(x > p_x1 && x < p_x2) { float xx = p_x1 - x; float yy = fabs((p_y2 - p_y1) * xx / (p_x2 - p_x1)); // region 1 if(fabs(y - p_y1) <= yy) { return 1; } else { return 0; } } if(x > p_x3 && x < p_x4) { float xx = p_x4 - x; float yy = fabs((p_y4 - p_y3) * xx / (p_x4 - p_x3)); // region 3 if(fabs(y - p_y4) <= yy) { return 1; } else { return 0; } } } if(y <= 0) { // region 4 if(x > p_x6 && x < p_x5) { return 1; } if(x > p_x5 && x < p_x4) { float xx = p_x4 - x; float yy = fabs((p_y5 - p_y4) * xx / (p_x5 - p_x4)); // region5 if(fabs(y - p_y4) <= yy) { return 1; } else { return 0; } } if(x > p_x1 && x < p_x6) { float xx = p_x1 - x; float yy = fabs((p_y1 - p_y6) * xx / (p_x1 - p_x6)); // region 6 if(fabs(y - p_y1) <= yy) { return 1; } else { return 0; } } } return 0; } //::IsCoreInside() //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////