StRoot: StEEmcPool/FgtSandbox1/StFgtClustFindMaker.cxx Source File

00001 // *-- Author : J.Balewski
00002 // 
00003 // $Id: StFgtClustFindMaker.cxx,v 1.2 2011/04/11 19:35:38 fisyak Exp $
00004 
00005 #include <TVector3.h>
00006 #include <TH2.h>
00007 #include <TF1.h>
00008 #include <TFile.h>
00009 #include <TLine.h>
00010 #include <TPolyLine.h>
00011 #include <TCrown.h>
00012 #include <TRandom3.h>
00013 
00014 #include "StFgtClustFindMaker.h"
00015 #include "StFgtSlowSimuMaker.h"
00016 
00017 ClassImp(StFgtClustFindMaker)
00018   
00019 
00020 //--------------------------------------------
00021 StFgtClustFindMaker::StFgtClustFindMaker(const char *name):StMaker(name){
00023   setHList(0);
00024   memset(hA,0,sizeof(hA));
00025   geom=new StFgtGeom();
00026   mRnd = new TRandom3(); // general use random generator
00027   //  mRnd->SetSeed(0); // activate, assure every set of data is different
00028 
00029   par_bx_valley=3; // min separation between 2 clusters
00030   par_cl_width=3; //min cluster width
00031   par_seedStripThres=0; // a.u. for valid cluster
00032   par_clusterMinAmpl=0.0; // a.u., min cluster ampl 
00033   par_stripNoiseSigma=0.0  ; // a.u.,added as gauss noise to every strip
00034 
00035 }
00036 
00037 //--------------------------------------------
00038 void 
00039 StFgtClustFindMaker::Clear(Option_t *) {
00040   int i;
00041   for(i=0;i<kFgtMxDisk;i++) {
00042     mRadClustList[i].clear();
00043     mPhiClustList[i].clear();
00044   }
00045   StMaker::Clear();
00046 }
00047 
00048 //--------------------------------------------
00049 StFgtClustFindMaker::~StFgtClustFindMaker(){
00050 
00051 }
00052 
00053 //_______________________________________________
00054 //________________________________________________
00055 void
00056 StFgtClustFindMaker::saveHisto(TString fname){
00057   TString outName=fname+".hist.root";
00058   TFile f( outName,"recreate");
00059   assert(f.IsOpen());
00060   printf("%d histos are written  to '%s' ...\n",HList->GetEntries(),outName.Data());
00061 
00062   HList->Write();
00063   f.Close();
00064 
00065 }
00066 
00067 
00068 //--------------------------------------------
00069 //--------------------------------------------
00070 //--------------------------------------------
00071 Int_t 
00072 StFgtClustFindMaker::Finish(){
00073   LOG_INFO<<"::Finish() \n"<<  endm; 
00074 
00075 
00076   return StMaker::Finish();
00077 }
00078 
00079 //--------------------------------------------
00080 //--------------------------------------------
00081 //--------------------------------------------
00082 Int_t 
00083 StFgtClustFindMaker::Make(){
00084   mInpEve++;
00085   
00086   LOG_INFO <<"::Make() inpEve="<<  mInpEve<<endm;
00087   StFgtSlowSimuMaker *ssMk=(StFgtSlowSimuMaker *) GetMaker("FgtSlowSimu");
00088   assert(ssMk);
00089 
00090   printf(" INPUT:  strips\n disk  Rad: #strip -> #clust   Phi: #strip -> #clust \n");
00091 
00092   // .... find 1D clusters in every plane of all disks
00093   for(int iDisk=0;iDisk<kFgtMxDisk;iDisk++) {
00094     int nRCl=findClust1D(ssMk->mRadAdcList[iDisk],mRadClustList[iDisk]);
00095     int nPCl=findClust1D(ssMk->mPhiAdcList[iDisk],mPhiClustList[iDisk]);
00096 
00097     printf("disk=%d  %d -->%d   %d -->%d \n",iDisk+1,ssMk->mRadAdcList[iDisk].size(),nRCl,ssMk->mPhiAdcList[iDisk].size(),nPCl);
00098     
00099     UInt_t i;
00100     int j;
00101 
00102     //::::::::::::::::process Rad-strips::::::::::::::::
00103     int nRQuad[kFgtMxQuad]; memset(nRQuad,0,sizeof(nRQuad));
00104 // delete the last cluster if strip ID out of boundary
00105     if(mRadClustList[iDisk].size() > 0) {  
00106       if(mRadClustList[iDisk][mRadClustList[iDisk].size()-1].fBin_mean >
00107                   geom->radStripLOCId_number() * kFgtMxQuad - 1)  {
00108           cout << "A Rad cluster deleted with fBin_mean = " << 
00109           mRadClustList[iDisk][mRadClustList[iDisk].size()-1].fBin_mean << endl;
00110 
00111           mRadClustList[iDisk].erase(mRadClustList[iDisk].end());
00112       }
00113     }
00114 
00115     for(i=0; i<mRadClustList[iDisk].size(); i++) {
00116       fgt_cluster1D *cl=&mRadClustList[iDisk][i];
00117       cl->position=geom->stripID2Rxy(cl->fBin_mean);
00118       int iRadID=(int) cl->fBin_mean; // get rid of fractional bin accuracy
00119       cl->iQuad=iRadID / geom->radStripLOCId_number(); 
00120        printf("icl=%d rRad=%.3f totAmpl=%.1f  totStrip=%d\n",i,cl->position,cl->totAmpl,cl->nbin);
00121 //      assert(  cl->iQuad <kFgtMxQuad);
00122       if(cl->iQuad >= 0 && cl->iQuad < kFgtMxQuad) {
00123         nRQuad[cl->iQuad]++;
00124         hA[5]->Fill(cl->peakAmpl);
00125         hA[7]->Fill(cl->nbin);
00126         hA[9]->Fill(cl->peakAmpl/cl->totAmpl);
00127       }
00128       else 
00129        cout << "ClusterFinder (Rad) error: iQuad = " << cl->iQuad << endl;
00130     }
00131 
00132 
00133     //::::::::::::::::process Phi-strips::::::::::::::::
00134     int nPQuad[kFgtMxQuad]; memset(nPQuad,0,sizeof(nPQuad));
00135 // delete the last cluster if strip ID out of boundary
00136     if(mPhiClustList[iDisk].size() > 0) {  
00137       if(mPhiClustList[iDisk][mPhiClustList[iDisk].size()-1].fBin_mean >
00138                   geom->phiStripLOCId_number() * kFgtMxQuad - 1 ) {
00139           cout << "A Phi cluster deleted with fBin_mean = " << 
00140           mPhiClustList[iDisk][mPhiClustList[iDisk].size()-1].fBin_mean << endl;
00141 
00142           mPhiClustList[iDisk].erase(mPhiClustList[iDisk].end());
00143       }
00144     }
00145 
00146     for(i=0; i<mPhiClustList[iDisk].size(); i++) {
00147       fgt_cluster1D *cl=&mPhiClustList[iDisk][i];
00148       cl->position=geom->stripID2PhiLab(cl->fBin_mean);
00149       int iPhiID=(int) cl->fBin_mean; // get rid of fractional bin accuracy
00150       cl->iQuad=iPhiID / geom->phiStripLOCId_number(); 
00151       printf("icl=%d rPhi/deg=%.3f totAmpl=%.1f totStrip=%d\n",i,cl->position/3.1417*180.,cl->totAmpl,cl->nbin);
00152       if(cl->iQuad >= 0 && cl->iQuad < kFgtMxQuad) {
00153         nPQuad[cl->iQuad]++;
00154 
00155         hA[6]->Fill(cl->peakAmpl);
00156         hA[8]->Fill(cl->nbin);
00157         hA[10]->Fill(cl->peakAmpl/cl->totAmpl);
00158       }
00159       else 
00160        cout << "ClusterFinder (Phi) error: iQuad = " << cl->iQuad << endl;
00161     }
00162 
00163 
00164     //... QA results for both planes
00165     hA[0]->Fill(20*iDisk+0,nRCl);
00166     hA[0]->Fill(20*iDisk+1,nPCl);
00167 
00168     for(j=0;j<kFgtMxQuad;j++) {
00169       hA[1]->Fill(nRQuad[j]);
00170       hA[2]->Fill(nPQuad[j]);
00171       hA[3]->Fill(nRQuad[j],nPQuad[j]);
00172       cout << "iquad nRQuad, nPQuad = " << j << " " << nRQuad[j] 
00173                                              << " " << nPQuad[j] << endl;
00174     }
00175    
00176   } // end of given disk
00177 
00178  return kStOK;
00179 }
00180 
00181  
00182 
00183 //--------------------------------------------
00184 //--------------------------------------------
00185 //--------------------------------------------
00186 Int_t 
00187 StFgtClustFindMaker::Init(){
00188   LOG_INFO<<"::Init() "<<  endm; 
00189   assert(HList);
00190   mInpEve=0;
00191   int i;
00192   LOG_INFO<<Form("::Init params minPeakWidth=%d, minPeakSepar=%d (# of strips), seedStripThres=%.1f a.u., clusterMinAmpl=%.1f a.u., stripNoiseSigma=%.1f a.u.",par_cl_width,par_bx_valley,par_seedStripThres,par_clusterMinAmpl,par_stripNoiseSigma)<<  endm;
00193   assert(par_clusterMinAmpl>=0);
00194   assert(par_stripNoiseSigma>=0);
00195 
00196   hA[0]=new TH1F("cl_Stat1D","Found 1D clusters, odd=Rad, even=Phi; x=20*disk",180, -0.5,179.5); 
00197   
00198   hA[1]=new TH1F("cl_rMul","Mult Rad-clusters, 1D , any disks; multiplicity/Quad/event",6,0.5,6.5); // zero is in underflow
00199   hA[2]=new TH1F("cl_pMul","Mult Phi-clusters, 1D , any disks; multiplicity/Quad/event",6,0.5,6.5); // zero is in underflow
00200   hA[3]=new TH2F("cl_rpMul","Mult per quadrant, 2x1D, any disks;Phi-Rad multiplicity/Quad/event; rad-mult",6,0.5,6.5,6,0.5,6.5); // zero is in underflow
00201   
00202 
00203   hA[5]= new TH1F("cl_RmxAmp"," Rad-strips 1D clust ; Max strip Amplitude (a.u.)",100,0,1000);
00204   hA[6]= new TH1F("cl_PmxAmp"," phi-strips 1D clust ; Max strip Amplitude  (a.u.)",100,0,1000);
00205 
00206   hA[7]= new TH1F("cl_Rwid"," Rad-strips 1D clust ;total cluster width ",10,0.5,10.5);
00207   hA[8]= new TH1F("cl_Pwid"," Phi-strips 1D clust ;total cluster width ",20,0.5,20.5);
00208 
00209   hA[9] = new TH1F("cl_Rpf"," Rad-strips 1D clust ; fraction of peak strip energy ",50,0.,1.);
00210   hA[10]= new TH1F("cl_Ppf"," Phi-strips 1D clust ; fraction of peak strip energy ",50,0.,1.);
00211 
00212  
00213 
00214   for(i=0;i<mxH;i++) 
00215     if(hA[i]) HList->Add(hA[i]);
00216   
00217   return StMaker::Init();
00218 }
00219 
00220 //--------------------------------------------
00221 //--------------------------------------------
00222 //--------------------------------------------
00223 int
00224 StFgtClustFindMaker::findClust1D(vector<fgt_strip> &sL, vector<fgt_cluster1D> &clL) {
00225   //  printf("AA %d\n",sL.size());
00226   if(sL.size()<=0) return 0;
00227   fgt_strip fake; // assures the last cluster gets saved by the code below
00228   fake.id=99999999; // set above highest strip
00229   sL.push_back(fake);
00230 
00231   int bx0=-999,bx1=bx0; // first and last bin of current cluster
00232   double sum=0, sumx=0, peakA=0;
00233   for(UInt_t i=0;i<sL.size();i++) {
00234     assert(sL[i].id>=bx1); //assumes strips are ordered  increasingly
00235     if(sL[i].id>=bx1+par_bx_valley) { // new clusters starts
00236       int nbin=bx1-bx0+1;
00237       if(nbin>=par_cl_width && peakA>par_seedStripThres) { 
00238         if(par_stripNoiseSigma>0) { // fold in noise
00239           // printf(" before noise sum=%.3f sumx=%.3f  meanX=%.3f\n",sum,sumx,sumx/sum);
00240           int k;
00241           for(k=bx0-1;k<=bx1+1;k++){
00242             double rndAdc=mRnd->Gaus(0,par_stripNoiseSigma);
00243             sum+=rndAdc;
00244             sumx+=rndAdc*(k+0.5); // use center of the bin
00245           }
00246           // printf("  after noise sum=%.3f sumx=%.3f  meanX=%.3f\n",sum,sumx,sumx/sum);
00247 
00248         } 
00249         if(sum > par_clusterMinAmpl) { // record old cluster if surviving noise
00250           fgt_cluster1D cl;
00251           cl.nbin=nbin;
00252           cl.fBin_mean=sumx/sum;
00253           cl.totAmpl=sum;
00254           cl.peakAmpl=peakA;
00255           printf("add cl=%d  fBin=%.3f, sum=%.1f peakA=%.1f width=%d\n",clL.size(),cl.fBin_mean,sum,peakA,nbin);
00256           clL.push_back(cl);
00257         }
00258       }
00259       sum=sumx=peakA=0;
00260       bx1=bx0=sL[i].id; 
00261     }
00262     sum+=sL[i].adc;
00263     sumx+=sL[i].adc*(sL[i].id+0.5); // use center of thebin
00264     if(peakA<sL[i].adc) peakA=sL[i].adc;
00265     bx1=sL[i].id;
00266     //printf("\nnc=%d i=%d bx=%d adc=%.1f sum=%f sumx=%.1f bx0=%d bx1=%d\n",clL.size(),i,sL[i].id,sL[i].adc,sum,sumx,bx0,bx1);
00267   }
00268   return clL.size(); 
00269 }
00270 
00271 
00274 
00275 // $Log: StFgtClustFindMaker.cxx,v $
00276 // Revision 1.2  2011/04/11 19:35:38  fisyak
00277 // Replace uint by UInt_t, use TMath
00278 //
00279 // Revision 1.1  2011/04/07 19:31:22  balewski
00280 // start
00281 //
00282 
00283 
00284  
00285 
00286