diff --git a/slsDetectorCalibration/moenchReadData.C b/slsDetectorCalibration/moenchReadData.C index 473a78f31..e98fea494 100644 --- a/slsDetectorCalibration/moenchReadData.C +++ b/slsDetectorCalibration/moenchReadData.C @@ -83,7 +83,6 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb THStack *hs=new THStack("hs",fformat); - int iev=0; TH2F *h1=new TH2F("h1",tit,nbins,hmin-0.5,hmax-0.5,NC*NR,-0.5,NC*NR-0.5); hs->Add(h1); @@ -113,7 +112,6 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb tall=filter->initEventTree(tit, &iFrame); - double tot, tl, tr, bl, br, v; #ifdef MY_DEBUG @@ -152,16 +150,16 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb if (hitfinder) { - filter->newFrame(); + filter->newFrame(); //calculate pedestals and common modes - if (cmsub) { + if (cmsub) { // cout << "cm" << endl; - for (ix=xmin-1; ixgetEventType(buff, ix, iy,0); - } - } + for (ix=xmin-1; ixgetEventType(buff, ix, iy,0); + } + } } // cout << "new frame " << endl; @@ -179,87 +177,14 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb #endif if (nf>1000) { - h1->Fill(filter->getClusterElement(0,0), iy+NR*ix); + h1->Fill(filter->getClusterTotal(1), iy+NR*ix); if (hitfinder) { - tot=0; - tl=0; - tr=0; - bl=0; - br=0; - // if (nf%1000==0 && ix==20 && iy==20) cout << " val="<< decoder->getClusterElement(0,0)<< endl; - if (thisEvent==PHOTON_MAX ) { nph++; - for (ir=-1; ir<2; ir++) { - for (ic=-1; ic<2; ic++) { - v=filter->getClusterElement(ic,ir); - // data[ic+1][ir+1]=v; - - - - tot+=v; - if (ir<1) { - - if (ic<1) { - bl+=v; - - } - if (ic>-1) { - br+=v; - } - } - - if (ir>-1) { - if (ic<1) { - tl+=v; - } - if (ic>-1) { - tr+=v; - } - } - } - } - - if (bl>br && bl>tl && bl>tr) { - h2->Fill(bl, iy+NR*ix); - if (bl>0) { - hetaX->Fill((filter->getClusterElement(0,0)+filter->getClusterElement(0,-1))/bl,iy+NR*ix); - hetaY->Fill((filter->getClusterElement(0,0)+filter->getClusterElement(-1,0))/bl,iy+NR*ix); - hetaX->Fill(1.-(filter->getClusterElement(0,0)+filter->getClusterElement(0,-1))/bl,iy+NR*(ix-1)); - hetaY->Fill(1.-(filter->getClusterElement(0,0)+filter->getClusterElement(-1,0))/bl,(iy-1)+NR*ix); - } - } else if (br>bl && br>tl && br>tr) { - h2->Fill(br, iy+NR*ix); - if (br>0) { - hetaX->Fill((filter->getClusterElement(0,0)+filter->getClusterElement(0,-1))/br,iy+NR*ix); - hetaY->Fill((filter->getClusterElement(0,0)+filter->getClusterElement(1,0))/br,iy+NR*ix); - hetaX->Fill(1.-(filter->getClusterElement(0,0)+filter->getClusterElement(0,-1))/br,iy+NR*(ix+1)); - hetaY->Fill(1.-(filter->getClusterElement(0,0)+filter->getClusterElement(1,0))/br,iy-1+NR*ix); - } - } else if (tl>br && tl>bl && tl>tr) { - h2->Fill(tl, iy+NR*ix); - if (tl>0) { - hetaX->Fill((filter->getClusterElement(0,0)+filter->getClusterElement(0,1))/tl,iy+NR*ix); - hetaY->Fill((filter->getClusterElement(0,0)+filter->getClusterElement(-1,0))/tl,iy+NR*ix); - hetaX->Fill(1.-(filter->getClusterElement(0,0)+filter->getClusterElement(0,1))/tl,iy+NR*(ix-1)); - hetaY->Fill(1.-(filter->getClusterElement(0,0)+filter->getClusterElement(-1,0))/tl,iy+1+NR*ix); - } - } else if (tr>bl && tr>tl && tr>br) { - h2->Fill(tr, iy+NR*ix); - if (tr>0) { - hetaX->Fill((filter->getClusterElement(0,0)+filter->getClusterElement(0,1))/tr,iy+NR*ix); - hetaY->Fill((filter->getClusterElement(0,0)+filter->getClusterElement(1,0))/tr,iy+NR*ix); - hetaX->Fill(1.-(filter->getClusterElement(0,0)+filter->getClusterElement(0,1))/tr,iy+NR*(ix+1)); - hetaY->Fill(1.-(filter->getClusterElement(0,0)+filter->getClusterElement(1,0))/tr,iy+1+NR*ix); - } - - - - } - - h3->Fill(tot, iy+NR*ix); + h2->Fill(filter->getClusterTotal(2), iy+NR*ix); + h3->Fill(filter->getClusterTotal(3), iy+NR*ix); iFrame=decoder->getFrameNumber(buff); tall->Fill(); diff --git a/slsDetectorCalibration/singlePhotonDetector.h b/slsDetectorCalibration/singlePhotonDetector.h index 363879d87..1582030d4 100644 --- a/slsDetectorCalibration/singlePhotonDetector.h +++ b/slsDetectorCalibration/singlePhotonDetector.h @@ -28,7 +28,15 @@ using namespace std; PHOTON=2, PHOTON_MAX=3, NEGATIVE_PEDESTAL=4, - UNDEFINED=-1 + UNDEFINED_EVENT=-1 + }; + + enum quadrant { + TOP_LEFT=0, + TOP_RIGHT=1, + BOTTOM_LEFT=2, + BOTTOM_RIGHT=3, + UNDEFINED_QUADRANT=-1 }; @@ -62,7 +70,7 @@ class singlePhotonDetector { int sign=1, commonModeSubtraction *cm=NULL, int nped=1000, - int nd=100) : det(d), nx(0), ny(0), stat(NULL), cmSub(cm), nDark(nd), eventMask(NULL),nSigma (nsigma), clusterSize(csize), clusterSizeY(csize), cluster(NULL), iframe(-1), dataSign(sign) { + int nd=100) : det(d), nx(0), ny(0), stat(NULL), cmSub(cm), nDark(nd), eventMask(NULL),nSigma (nsigma), clusterSize(csize), clusterSizeY(csize), cluster(NULL), iframe(-1), dataSign(sign), quad(UNDEFINED_QUADRANT), tot(0), quadTot(0) { det->getDetectorSize(nx,ny); @@ -94,7 +102,7 @@ class singlePhotonDetector { void newDataSet(){iframe=-1; for (int iy=0; iyClear(); }; /** resets the eventMask to undefined and the commonModeSubtraction */ - void newFrame(){iframe++; for (int iy=0; iynewFrame();}; + void newFrame(){iframe++; for (int iy=0; iynewFrame();}; /** sets the commonModeSubtraction algorithm to be used @@ -184,8 +192,13 @@ class singlePhotonDetector { eventType getEventType(char *data, int ix, int iy, int cm=0) { // eventType ret=PEDESTAL; - double tot=0, max=0; + double max=0, tl=0, tr=0, bl=0,br=0, v; // cout << iframe << endl; + + tot=0; + quadTot=0; + quad=UNDEFINED_QUADRANT; + if (iframeset_data(dataSign*(det->getValue(data, ix, iy)-getPedestal(ix,iy,cm)), 0,0 ); // cout << "=" << endl; } - return UNDEFINED; + return UNDEFINED_EVENT; } @@ -214,12 +227,22 @@ class singlePhotonDetector { for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) { if ((iy+ir)>=0 && (iy+ir)=0 && (ix+ic)set_data(dataSign*(det->getValue(data, ix+ic, iy+ir)-getPedestal(ix+ic,iy+ir,cm)), ic, ir ); - tot+=cluster->get_data(ic,ir); + v=cluster->get_data(ic,ir); + tot+=v; + if (ir<=0 && ic<=0) + bl+=v; + if (ir<=0 && ic>=0) + br+=v; + if (ir>=0 && ic<=0) + tl+=v; + if (ir>=0 && ic>=0) + tr+=v; + if (cluster->get_data(ic,ir)>max) { - max=cluster->get_data(ic,ir); + max=v; } if (ir==0 && ic==0) { - if (cluster->get_data(ic,ir)>nSigma*cluster->rms) { + if (v>nSigma*cluster->rms) { eventMask[iy][ix]=PHOTON; } else if (cluster->get_data(ic,ir)<-nSigma*cluster->rms) eventMask[iy][ix]=NEGATIVE_PEDESTAL; @@ -233,31 +256,54 @@ class singlePhotonDetector { } else if (eventMask[iy][ix]==PHOTON) { if (cluster->get_data(0,0)>=max) { eventMask[iy][ix]=PHOTON_MAX; - /* if (iframe%1000==0) { */ -/* for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) { */ -/* for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) { */ -/* if ((iy+ir)>=0 && (iy+ir)=0 && (ix+ic)get_data(ic,ir) << " "; */ -/* } */ -/* } */ -/* } */ -/* cout << endl;; */ -/* } */ } } else if (eventMask[iy][ix]==PEDESTAL) { if (cm==0) addToPedestal(det->getValue(data, ix, iy),ix,iy); } - // } - + + + if (bl>=br && bl>=tl && bl>=tr) { + quad=BOTTOM_LEFT; + quadTot=bl; + } else if (br>=bl && br>=tl && br>=tr) { + quad=BOTTOM_RIGHT; + quadTot=br; + } else if (tl>=br && tl>=bl && tl>=tr) { + quad=TOP_LEFT; + quadTot=tl; + } else if (tr>=bl && tr>=tl && tr>=br) { + quad=TOP_RIGHT; + quadTot=tr; + } + + return eventMask[iy][ix]; }; + /**< + retrurns the total signal in a cluster + \param size cluser size should be 1,2 or 3 + \returns cluster center if size=1, sum of the maximum quadrant if size=2, total of the cluster if size=3 or anything else + */ + double getClusterTotal(int size) { + switch (size) { + case 1: + return getClusterElement(0,0); + case 2: + return quadTot; + default: + return tot; + }; + }; + /**< + retrurns the quadrant with maximum signal + \returns quadrant where the cluster is located */ + + quadrant getQuadrant() {return quad;}; /** sets/gets number of samples for moving average pedestal calculation \param i number of samples to be set (0 or negative gets) @@ -279,6 +325,7 @@ class singlePhotonDetector { */ eventType getEventMask(int ic, int ir=0){return eventMask[ir][ic];}; + #ifdef MYROOT1 /** generates a tree and maps the branches \param tname name for the tree @@ -298,8 +345,8 @@ class singlePhotonDetector { char tit[100]; sprintf(tit,"data[%d]/D",clusterSize*clusterSizeY); tall->Branch("data",cluster->data,tit); - tall->Branch("pedestal",&(cluster->ped),"pedestal/D"); - tall->Branch("rms",&(cluster->rms),"rms/D"); + // tall->Branch("pedestal",&(cluster->ped),"pedestal/D"); + // tall->Branch("rms",&(cluster->rms),"rms/D"); return tall; }; #endif @@ -322,7 +369,11 @@ class singlePhotonDetector { single_photon_hit *cluster; /**< single photon hit data structure */ int iframe; /**< frame number (not from file but incremented within the dataset every time newFrame is called */ int dataSign; /**< sign of the data i.e. 1 if photon is positive, -1 if negative */ - + quadrant quad; /**< quadrant where the photon is located */ + double tot; /**< sum of the 3x3 cluster */ + double quadTot; /**< sum of the maximum 2x2cluster */ + + };