61 _nsamples(iConfig.getUntrackedParameter<unsigned
int>(
"nSamples", 10)),
62 _presample(iConfig.getUntrackedParameter<unsigned
int>(
"nPresamples", 2)),
63 _firstsample(iConfig.getUntrackedParameter<unsigned
int>(
"firstSample", 1)),
64 _lastsample(iConfig.getUntrackedParameter<unsigned
int>(
"lastSample", 2)),
65 _nsamplesPN(iConfig.getUntrackedParameter<unsigned
int>(
"nSamplesPN", 50)),
66 _presamplePN(iConfig.getUntrackedParameter<unsigned
int>(
"nPresamplesPN", 6)),
67 _firstsamplePN(iConfig.getUntrackedParameter<unsigned
int>(
"firstSamplePN", 7)),
68 _lastsamplePN(iConfig.getUntrackedParameter<unsigned
int>(
"lastSamplePN", 8)),
69 _timingcutlow(iConfig.getUntrackedParameter<unsigned
int>(
"timingCutLow", 2)),
70 _timingcuthigh(iConfig.getUntrackedParameter<unsigned
int>(
"timingCutHigh", 9)),
71 _timingquallow(iConfig.getUntrackedParameter<unsigned
int>(
"timingQualLow", 3)),
72 _timingqualhigh(iConfig.getUntrackedParameter<unsigned
int>(
"timingQualHigh", 8)),
73 _ratiomincutlow(iConfig.getUntrackedParameter<double>(
"ratioMinCutLow", 0.4)),
74 _ratiomincuthigh(iConfig.getUntrackedParameter<double>(
"ratioMinCutHigh", 0.95)),
75 _ratiomaxcutlow(iConfig.getUntrackedParameter<double>(
"ratioMaxCutLow", 0.8)),
76 _presamplecut(iConfig.getUntrackedParameter<double>(
"presampleCut", 5.0)),
77 _niter(iConfig.getUntrackedParameter<unsigned
int>(
"nIter", 5)),
78 _noise(iConfig.getUntrackedParameter<double>(
"noise", 2.0)),
79 _ecalPart(iConfig.getUntrackedParameter<
std::
string>(
"ecalPart",
"EB")),
80 _saveshapes(iConfig.getUntrackedParameter<
bool>(
"saveShapes",
true)),
81 _docorpn(iConfig.getUntrackedParameter<
bool>(
"doCorPN",
false)),
82 _fedid(iConfig.getUntrackedParameter<
int>(
"fedID", -999)),
83 _saveallevents(iConfig.getUntrackedParameter<
bool>(
"saveAllEvents",
false)),
84 _qualpercent(iConfig.getUntrackedParameter<double>(
"qualPercent", 0.2)),
85 _debug(iConfig.getUntrackedParameter<
int>(
"debug", 0)),
113 channelIteratorEE(0),
145 for (
unsigned int j = 0;
j <
nCrys;
j++) {
157 for (
unsigned int j = 0;
j <
nMod;
j++) {
213 ADCfile +=
"/APDSamplesLaser.root";
216 APDfile +=
"/APDPNLaserAllEvents.root";
220 for (
unsigned int i = 0;
i <
nCrys;
i++) {
222 name <<
"ADCTree" <<
i + 1;
252 stringstream namefile1;
253 namefile1 <<
resdir_ <<
"/SHAPE_LASER.root";
256 stringstream namefile2;
257 namefile2 <<
resdir_ <<
"/APDPN_LASER.root";
260 stringstream namefile3;
261 namefile3 <<
resdir_ <<
"/MATACQ.root";
270 cout <<
" ERROR! No matacq shape available: analysis aborted !" << endl;
290 DCCHeader = pDCCHeader.
product();
316 cout <<
" Wrong ecalPart in cfg file " << endl;
336 TheMapping = ecalmapping.
product();
338 std::cerr <<
"Error! can't get the product EcalMappingRcd" << std::endl;
349 int fed = headerItr->fedId();
353 runType = headerItr->getRunType();
354 runNum = headerItr->getRunNumber();
355 event = headerItr->getLV1();
357 dccID = headerItr->getDccInTCCCommand();
358 fedID = headerItr->fedId();
380 if (iter ==
colors.end()) {
408 unsigned int samplemax = 0;
411 map<int, vector<double> > allPNAmpl;
412 map<int, vector<double> > allPNGain;
420 cout <<
"-- debug test -- Inside PNDigi - pnID=" << pnDetId.
iPnId() <<
", dccID=" << pnDetId.
iDCCId() << endl;
430 for (
int samId = 0; samId < (*pnItr).size(); samId++) {
431 pn[samId] = (*pnItr).sample(samId).adc();
432 pnG[samId] = (*pnItr).sample(samId).gainId();
440 cout <<
"PN gain different from 1" << endl;
448 if (chi2pn == 101 || chi2pn == 102 || chi2pn == 103)
465 cout <<
"-- debug -- Inside PNDigi - PNampl=" <<
pnAmpl <<
", PNgain=" << pnGain << endl;
484 EBDetId id_crystal(digiItr->id());
488 int etaG = id_crystal.ieta();
489 int phiG = id_crystal.iphi();
493 int etaL = LocalCoord.first;
494 int phiL = LocalCoord.second;
497 int xtal = elecid_crystal.
xtalId();
505 unsigned int MyPn0 = pnpair.first;
506 unsigned int MyPn1 = pnpair.second;
516 <<
" modules:" <<
modules.size() << endl;
523 for (
unsigned int i = 0;
i < (*digiItr).size(); ++
i) {
525 adc[
i] = samp_crystal.adc();
526 adcG[
i] = samp_crystal.gainId();
548 int mem0 =
Mem->
Mem(lmr, 0);
549 int mem1 =
Mem->
Mem(lmr, 1);
551 if (allPNAmpl[mem0].
size() > MyPn0)
552 pn0 = allPNAmpl[mem0][MyPn0];
555 if (allPNAmpl[mem1].
size() > MyPn1)
556 pn1 = allPNAmpl[mem1][MyPn1];
578 EEDetId id_crystal(digiItr->id());
582 int etaG = id_crystal.iy();
583 int phiG = id_crystal.ix();
585 int iX = (phiG - 1) / 5 + 1;
586 int iY = (etaG - 1) / 5 + 1;
599 unsigned int MyPn0 = pnpair.first;
600 unsigned int MyPn1 = pnpair.second;
610 setGeomEE(etaG, phiG, iX, iY,
iZ,
module,
tower, ch, apdRefTT, channel, lmr);
614 <<
" modules:" <<
modules.size() << endl;
619 if ((*digiItr).size() > 10)
620 cout <<
"SAMPLES SIZE > 10!" << (*digiItr).size() << endl;
624 for (
unsigned int i = 0;
i < (*digiItr).size(); ++
i) {
626 adc[
i] = samp_crystal.adc();
627 adcG[
i] = samp_crystal.gainId();
650 int mem0 =
Mem->
Mem(lmr, 0);
651 int mem1 =
Mem->
Mem(lmr, 1);
653 if (allPNAmpl[mem0].
size() > MyPn0)
654 pn0 = allPNAmpl[mem0][MyPn0];
657 if (allPNAmpl[mem1].
size() > MyPn1)
658 pn1 = allPNAmpl[mem1][MyPn1];
681 cout <<
"\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << endl;
682 cout <<
"\t+=+ WARNING! NO MATACQ +=+" << endl;
683 cout <<
"\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << endl;
714 system(del.str().c_str());
715 cout <<
" No Laser Events " << endl;
722 double BadGainEvtPercentage = 0.0;
723 double BadTimingEvtPercentage = 0.0;
725 int nChanBadGain = 0;
726 int nChanBadTiming = 0;
728 for (
unsigned int i = 0;
i <
nCrys;
i++) {
743 double BadGainChanPercentage = double(nChanBadGain) / double(
nCrys);
744 double BadTimingChanPercentage = double(nChanBadTiming) / double(
nCrys);
754 cout <<
"\n\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << endl;
755 cout <<
"\t+=+ Analyzing laser data: getting APD, PN, APD/PN, PN/PN +=+" << endl;
758 cout <<
"\t+=+ ............................ WARNING! APD GAIN WAS NOT 1 +=+" << endl;
760 cout <<
"\t+=+ ............................ WARNING! TIMING WAS BAD +=+" << endl;
767 for (
unsigned int i = 0;
i <
nCrys;
i++) {
769 name <<
"APDTree" <<
i + 1;
787 APDtrees[
i]->Branch(
"flagfit", &flagfit,
"flagfit/I");
803 APDtrees[
i]->SetBranchAddress(
"flagfit", &flagfit);
808 for (
unsigned int iref = 0; iref <
nRefChan; iref++) {
809 for (
unsigned int imod = 0; imod <
nMod; imod++) {
812 stringstream nameref;
813 nameref <<
"refAPDTree" << imod <<
"_" << iref;
815 RefAPDtrees[iref][jmod] =
new TTree(nameref.str().c_str(), nameref.str().c_str());
834 unsigned int nCol =
colors.size();
839 for (
unsigned int iM = 0; iM <
nMod; iM++) {
840 unsigned int iMod =
modules[iM] - 1;
841 for (
unsigned int ich = 0; ich <
nPNPerMod; ich++) {
842 for (
unsigned int icol = 0; icol < nCol; icol++) {
854 for (
unsigned int iCry = 0; iCry <
nCrys; iCry++) {
855 for (
unsigned int icol = 0; icol < nCol; icol++) {
863 if (
ADCtrees[iCry]->GetEntries(
cut.str().c_str()) < 10)
867 unsigned int iMod =
iModule[iCry] - 1;
876 Long64_t nbytes = 0, nb = 0;
877 for (Long64_t jentry = 0; jentry <
ADCtrees[iCry]->GetEntriesFast(); jentry++) {
878 nb =
ADCtrees[iCry]->GetEntry(jentry);
889 unsigned int iCol = 0;
890 for (
unsigned int i = 0;
i < nCol;
i++) {
914 cout <<
"-- debug test -- endJob -- apdAmpl:" <<
apdAmpl <<
" apdTime:" <<
apdTime << endl;
917 if (pn0 < 10 && pn1 > 10) {
919 }
else if (pn1 < 10 && pn0 > 10) {
922 pnmean = 0.5 * (
pn0 +
pn1);
925 cout <<
"-- debug test -- endJob -- pnMean:" << pnmean << endl;
931 for (
unsigned int ichan = 0; ichan <
nPNPerMod; ichan++) {
942 cout <<
"-- debug test -- endJob -- filling APDTree" << endl;
955 for (
unsigned int ir = 0; ir <
nRefChan; ir++) {
957 cout <<
"-- debug test -- ir:" << ir <<
" tt:" <<
towerID <<
" refmap:" <<
apdRefMap[ir][iMod + 1]
958 <<
" iCry:" << iCry << endl;
962 cout <<
"-- debug test -- cut passed " << endl;
968 cout <<
"-- debug test -- apdAmplA=" <<
apdAmplA << endl;
970 cout <<
"-- debug test -- apdAmplB=" <<
apdAmplB << endl;
972 cout <<
"-- debug test -- color=" <<
color <<
", event:" <<
event <<
", ir:" << ir
973 <<
" tt-1:" <<
towerID - 1 << endl;
978 cout <<
"-- debug test -- tree filled" <<
event << endl;
990 cout <<
"-- debug test -- endJob -- after apdAmpl Loop" << endl;
997 system(del.str().c_str());
1004 for (
unsigned int iColor = 0; iColor < nCol; iColor++) {
1005 stringstream nametree;
1006 nametree <<
"APDCol" <<
colors.at(iColor);
1007 stringstream nametree2;
1008 nametree2 <<
"PNCol" <<
colors.at(iColor);
1010 restrees[iColor] =
new TTree(nametree.str().c_str(), nametree.str().c_str());
1011 respntrees[iColor] =
new TTree(nametree2.str().c_str(), nametree2.str().c_str());
1020 restrees[iColor]->Branch(
"APD", &
APD,
"APD[6]/D");
1064 for (
unsigned int iM = 0; iM <
nMod; iM++) {
1065 unsigned int iMod =
modules[iM] - 1;
1067 for (
unsigned int ich = 0; ich <
nPNPerMod; ich++) {
1068 for (
unsigned int icol = 0; icol < nCol; icol++) {
1078 for (
unsigned int imod = 0; imod <
nMod; imod++) {
1089 for (
unsigned int iCry = 0; iCry <
nCrys; iCry++) {
1090 unsigned int iMod =
iModule[iCry] - 1;
1095 for (
unsigned int iCol = 0; iCol < nCol; iCol++) {
1096 std::vector<double> lowcut;
1097 std::vector<double> highcut;
1106 lowcut.push_back(cutMin);
1107 highcut.push_back(cutMax);
1111 lowcut.push_back(cutMin);
1112 highcut.push_back(cutMax);
1133 Long64_t nbytes = 0, nb = 0;
1134 for (Long64_t jentry = 0; jentry <
APDtrees[iCry]->GetEntriesFast(); jentry++) {
1135 nb =
APDtrees[iCry]->GetEntry(jentry);
1139 if (pn0 < 10 && pn1 > 10) {
1141 }
else if (pn1 < 10 && pn0 > 10) {
1144 pnmean = 0.5 * (
pn0 +
pn1);
1149 unsigned int iCol = 0;
1150 for (
unsigned int i = 0;
i < nCol;
i++) {
1161 for (
unsigned int ichan = 0; ichan <
nPNPerMod; ichan++) {
1170 cout <<
"-- debug test -- LastLoop event:" <<
event <<
" apdAmpl:" <<
apdAmpl << endl;
1174 for (
unsigned int iRef = 0; iRef <
nRefChan; iRef++) {
1179 cout <<
"-- debug test -- LastLoop apdAmplA:" <<
apdAmplA <<
" apdAmplB:" <<
apdAmplB <<
", event:" <<
event
1180 <<
", eventref:" <<
eventref << endl;
1196 for (
unsigned int iColor = 0; iColor < nCol; iColor++) {
1197 std::vector<double> apdvec =
APDAnal[iCry][iColor]->
getAPD();
1205 for (
unsigned int i = 0;
i < apdvec.size();
i++) {
1206 APD[
i] = apdvec.at(
i);
1232 cout <<
"-- debug test -- endJob -- APD[0]" <<
APD[0] <<
" APDoPN[0] " <<
APDoPN[0] <<
" APDoAPDA[0] "
1241 for (
unsigned int iM = 0; iM <
nMod; iM++) {
1242 unsigned int iMod =
modules[iM] - 1;
1246 for (
unsigned int ch = 0; ch <
nPNPerMod; ch++) {
1253 for (
unsigned int iColor = 0; iColor < nCol; iColor++) {
1254 std::vector<double> pnvec =
PNAnal[iMod][ch][iColor]->
getPN();
1255 std::vector<double> pnopnvec =
PNAnal[iMod][ch][iColor]->
getPNoPN();
1256 std::vector<double> pnopn0vec =
PNAnal[iMod][ch][iColor]->
getPNoPN0();
1257 std::vector<double> pnopn1vec =
PNAnal[iMod][ch][iColor]->
getPNoPN1();
1259 for (
unsigned int i = 0;
i < pnvec.size();
i++) {
1260 PN[
i] = pnvec.at(
i);
1267 cout <<
"-- debug test -- endJob -- filling pn results'tree: PN[0]:" <<
PN[0] <<
" iModule:" << iMod
1268 <<
" iColor:" << iColor <<
" ch:" << ch << endl;
1285 system(del2.str().c_str());
1297 for (
unsigned int i = 0;
i < nCol;
i++) {
1304 cout <<
"\t+=+ .................................................. done +=+" << endl;
1305 cout <<
"\t+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+" << endl;
1317 int doesMatFileExist = 0;
1318 int doesMatShapeExist = 0;
1320 TProfile* laserShape =
nullptr;
1321 test2 = fopen(
matfile.c_str(),
"r");
1323 doesMatFileExist = 1;
1325 TFile* MatShapeFile;
1326 if (doesMatFileExist == 1) {
1327 MatShapeFile =
new TFile(
matfile.c_str());
1328 laserShape = (TProfile*)MatShapeFile->Get(
"shapeLaser");
1330 doesMatShapeExist = 1;
1331 double y = laserShape->Integral(
"w");
1333 laserShape->Scale(1.0 /
y);
1336 cout <<
" ERROR! Matacq shape file not found !" << endl;
1338 if (doesMatShapeExist)
1344 int doesElecFileExist = 0;
1348 doesElecFileExist = 1;
1350 TFile* ElecShapesFile;
1351 TH1D* elecShape =
nullptr;
1353 if (doesElecFileExist == 1) {
1354 ElecShapesFile =
new TFile(
elecfile_.c_str());
1356 name <<
"MeanElecShape";
1357 elecShape = (TH1D*)ElecShapesFile->Get(
name.str().c_str());
1358 if (elecShape && doesMatShapeExist == 1) {
1359 double x = elecShape->GetMaximum();
1361 elecShape->Scale(1.0 /
x);
1368 cout <<
" ERROR! Elec shape file not found !" << endl;
1374 unsigned int nBins =
int(laserShape->GetEntries());
1376 double elec_jj, laser_iiMinusjj;
1380 unsigned int nBins2 =
int(elecShape->GetNbinsX());
1382 if (nBins2 <
nBins) {
1383 cout <<
"EcalLaserAnalyzer2::getShapes: wrong configuration of the shapes' number of bins" << std::endl;
1389 name <<
"PulseShape";
1395 for (
int ii = 0;
ii < 50;
ii++) {
1400 for (
unsigned int ii = 0;
ii <
nBins - 50;
ii++) {
1402 for (
unsigned int jj = 0;
jj <
ii;
jj++) {
1403 elec_jj = elecShape->GetBinContent(
jj + 1);
1404 laser_iiMinusjj = laserShape->GetBinContent(
ii -
jj + 1);
1405 sum_jj += elec_jj * laser_iiMinusjj;
1430 system(del.str().c_str());
1437 int etaG,
int phiG,
int module,
int tower,
int strip,
int xtal,
int apdRefTT,
int channel,
int lmr) {
1449 for (
unsigned int iref = 0; iref <
nRefChan; iref++) {
1470 int etaG,
int phiG,
int iX,
int iY,
int iZ,
int module,
int tower,
int ch,
int apdRefTT,
int channel,
int lmr) {
1482 for (
unsigned int iref = 0; iref <
nRefChan; iref++) {