StRoot: StPmdCalibrationMaker/StPmdCalibConstMaker.cxx Source File

00001 /*************************************************
00002  *
00003  * StPmdCalibConstMaker.cxx$
00004  * Author: Raghunath Sahoo, IOP, Bhubaneswar, Subhasis Chattopadhyay, VECC, Kolkata.
00005  *************************************************
00006  *
00007  * Description: Base class for Pmd Calibration Constant Maker
00008  *
00009  *************************************************
00010  *
00011  * Initial version : 28th aug 2003
00012  *  The code is tested to work in simulation chain, and tested with 
00013  *  data taken with PMD last year. Numbers are written to DB, and tested
00014  *  to read them back. 
00015  *
00016  *   Some changes are forseen in data storing once a large amount 
00017  *   PMD data will be run to generate the numbers.
00018  *
00019  *************************************************/
00020 
00021 /* !
00022  
00023 This Maker obtains calibrations constants for all cells in PMD and CPV.
00024 This collects MIP spectra for each cell and at the end it fits them with 
00025 Landau function and then generates the calibration constants. these calibration 
00026 constants are then stored in calibration database
00027 */
00028 
00029 #include"Stsstream.h"
00030 #include<assert.h>
00031 #include<math.h>
00032 #include"TROOT.h"
00033 #include"TMath.h"
00034 #include<TRandom.h>
00035 #include<StMessMgr.h>
00036 #include<TFile.h>
00037 #include<TF1.h>
00038 #include "StPmdUtil/StPmdGeom.h"
00039 #include "StPmdCalibConstMaker.h"
00040 #include "StPmdUtil/StPmdDetector.h" 
00041 #include "StPmdUtil/StPmdModule.h" 
00042 #include "StPmdUtil/StPmdCollection.h" 
00043 #include "StPmdUtil/StPmdHit.h"         
00044 //...........................................................
00045 #include "TStopwatch.h"
00046 #include "StarClassLibrary/SystemOfUnits.h"
00047 #include "time.h"
00048 #include "TDatime.h"
00049 
00050 #include "StDbLib/StDbManager.hh"
00051 #include "StDbLib/StDbTable.h"
00052 #include "StDbLib/StDbConfigNode.hh"
00053 
00054 #ifndef ST_NO_NAMESPACES
00055 using units::tesla;
00056 #endif 
00057 
00058 
00059 
00060 ClassImp(StPmdCalibConstMaker)
00061 
00062   
00069   Double_t MipArray[PMD_CRAMS_MAX*2][PMD_ROW_MAX][PMD_COL_MAX][MIP_CH_MAX];
00070 
00071   //___________________------------
00073 Int_t StPmdCalibConstMaker::neibx[PMD_CELL_NEIGHBOUR] = {1,0,-1,-1,0,1};
00074 
00075 Int_t StPmdCalibConstMaker::neiby[PMD_CELL_NEIGHBOUR] = {0,1,1,0,-1,-1};
00076 
00077 Int_t StPmdCalibConstMaker::imax[2*PMD_CRAMS_MAX] = {48,72,48,72,72,72,48,48,48,72,72,96,48,72,48,72,72,72,48,48,48,72,72,96};
00078 //for col
00079 Int_t StPmdCalibConstMaker::jmax[2*PMD_CRAMS_MAX] = {48,48,72,72,48,48,48,72,48,72,48,72,48,48,72,72,48,48,48,72,48,72,48,72};
00080 //for row
00081 
00082 //-------------------- 
00083 StPmdCalibConstMaker::StPmdCalibConstMaker(const char *name):StMaker(name)
00084 {
00085   // mSaveCalibToDB = kFALSE;  //To be checked
00086   mSaveCalibToDB = kTRUE; 
00087   mOptHist = kTRUE; 
00088   mDate=0;
00089   mTime=0;
00090   InitMipParams();
00091   mPmdGeom=new StPmdGeom();
00092   mPmdDbUtil=new StPmdDBUtil();
00093 }
00094 //----------------------
00095 
00096 StPmdCalibConstMaker::~StPmdCalibConstMaker()
00097 {
00098   cout << "**** I am in StPmdCalibConstMaker::~StPmdCalibConstMaker()"<< endl;
00099   
00100   if(mPmdGeom){mPmdGeom=0;delete mPmdGeom;}
00101   if(mPmdDbUtil){mPmdDbUtil=0;delete mPmdDbUtil;}
00102   if(mOptHist)ClearHists();
00103 }
00104 
00105 //---------------------
00106 Int_t StPmdCalibConstMaker::Init()
00107 {
00108   //Initialize date and time
00109   if(mOptHist)BookHistograms();
00110   
00111   //---------------------------------------------------
00112   // Newly included from this
00113   
00121   ClearMipArray();
00122 
00123   return StMaker::Init();
00124 }
00125 
00126 //--------------------------------------
00127 void StPmdCalibConstMaker::BookHistograms()
00128 {
00129 
00130   for(Int_t sm=0;sm<2*PMD_CRAMS_MAX;sm++){
00131     for(Int_t i =0;i<imax[sm];i++){
00132       for(Int_t j =0;j<jmax[sm];j++){
00133         char text[40],title[80];
00134         sprintf(text,"calib_%02d_%02d_%02d",sm,i,j);
00135         sprintf(title,"MIP plot for %02d_%02d_%2d",sm,i,j);
00136         mMipEnergy[sm][i][j]=new TH1F(text,title,200,0.,150.);
00137       }
00138     }
00139   }
00140 }
00141 
00142 
00143 //-------------------------------------------
00144 
00145 Int_t StPmdCalibConstMaker::Make() 
00146 {
00147   if (mDate==0) {mDate = GetDate();}
00148   if (mTime==0) {mTime = GetTime();}
00149 
00150   TDataSet* CalibIn = GetDataSet("pmdReader");
00151   if (!CalibIn){CalibIn = GetDataSet("PmdSimulator");}
00152   if (CalibIn){
00153     StPmdCollection *CalibHit = (StPmdCollection*)CalibIn->Find("PmdCollection");
00154     if(CalibHit){
00155       StPmdDetector * cpv_det = CalibHit->detector(Int_t(1));
00156       StPmdDetector * pmd_det = CalibHit->detector(Int_t(0));
00157       GetIsoHit(pmd_det,cpv_det);
00158     }
00159   }
00160   else{cout<<"Could not find Simulator/Reader"<<endl; return kStOK;}
00161   return kStOK;
00162 }
00163 
00164 
00165 void StPmdCalibConstMaker::GetIsoHit(StPmdDetector* pmd_det, 
00166                                      StPmdDetector* cpv_det)
00167 {
00170   Int_t xpad,ypad,gsuper,idet;
00171   Float_t edep;
00172   Float_t d0[PMD_ROW_MAX][PMD_COL_MAX];
00173   
00174   if(cpv_det){  
00175     for(Int_t id=0;id<PMD_CRAMS_MAX;id++){   
00176 
00177       for(Int_t j=0;j<PMD_ROW_MAX;j++){
00178         for(Int_t k=0;k<PMD_COL_MAX;k++){
00179           d0[j][k]=0;  
00180         }
00181       }
00182       
00183       
00184       StPmdModule * cpv_mod=cpv_det->module(id);
00185       if(cpv_det->module_hit(id)>0){
00186         Int_t nmh=cpv_det->module_hit(id);  
00187         TIter next(cpv_mod->Hits());
00188         StPmdHit *spmhit0;   
00189         for(Int_t im=0; im<nmh; im++){
00190           spmhit0 = (StPmdHit*)next();
00191           if(spmhit0)
00192             {
00193               ypad=spmhit0->Row();          
00194               xpad=spmhit0->Column();       
00195               edep=spmhit0->Edep();         
00196               Int_t adc=spmhit0->Adc();     
00197               gsuper = spmhit0->Gsuper();   
00198               idet=spmhit0->SubDetector();  
00199               if(idet==2){
00200 if((xpad>0 && xpad<=PMD_ROW_MAX) && (ypad>0 && ypad <=PMD_COL_MAX))d0[xpad-1][ypad-1]=adc;  
00201               }
00202             }
00203         }
00204         
00206         
00207         for(Int_t i=0;i<imax[id];i++){
00208           for(Int_t j=0;j<jmax[id];j++){
00209             if(d0[i][j] > 0.){
00210               Int_t count = 0;
00211               for(Int_t ii=0;ii<PMD_CELL_NEIGHBOUR;ii++){
00212                 Int_t id1 = i + neibx[ii];
00213                 Int_t jd1 = j + neiby[ii];
00214                 if(d0[id1][jd1]==0.){
00215                   count++;
00216                   if(count == PMD_CELL_NEIGHBOUR){//finding isolated cell
00217                     gsuper = spmhit0->Gsuper(); 
00218                     gsuper = gsuper-1;
00219                     if(mOptHist)mMipEnergy[gsuper][i][j]->Fill(d0[i][j]);
00220                     Int_t channel=(Int_t)d0[i][j];
00221                     if(channel<MIP_CH_MAX)MipArray[gsuper][i][j][channel]++;
00222                     
00223                   }
00224                 }
00225               }
00226             }
00227           }
00228         }
00229         
00230       }// if-- module hit
00231     }// for supermodule
00232   }// if---CPV detector
00233 
00234 
00235 
00236   Float_t d1[PMD_ROW_MAX][PMD_COL_MAX];
00237   
00238   if(pmd_det){  
00239     for(Int_t id=0;id<PMD_CRAMS_MAX;id++){   
00240 
00241       for(Int_t j=0;j<PMD_ROW_MAX;j++){
00242         for(Int_t k=0;k<PMD_COL_MAX;k++){
00243           d1[j][k]=0;  
00244         }
00245       }
00246       
00247       StPmdModule * pmd_mod=pmd_det->module(id);
00248       if(pmd_det->module_hit(id)>0){
00249         Int_t nmh=pmd_det->module_hit(id);  
00250         TIter next(pmd_mod->Hits());
00251         StPmdHit *spmhit1;   
00252         for(Int_t im=0; im<nmh; im++){
00253           spmhit1 = (StPmdHit*)next();
00254           if(spmhit1)
00255             {
00256               ypad=spmhit1->Row();          
00257               xpad=spmhit1->Column();       
00258               edep=spmhit1->Edep();         
00259               Int_t adc=spmhit1->Adc();     
00260               gsuper = spmhit1->Gsuper();   
00261               idet=spmhit1->SubDetector();  
00262               if(idet==1){
00263 if((xpad>0 && xpad<=PMD_ROW_MAX) && (ypad>0 && ypad <=PMD_COL_MAX))d1[xpad-1][ypad-1]=adc;  
00264               }
00265             }
00266         }
00267 
00269         
00270         for(Int_t i=0;i<imax[id];i++){
00271           for(Int_t j=0;j<jmax[id];j++){
00272             if(d1[i][j] > 0.){
00273               Int_t count = 0;
00274               for(Int_t ii=0;ii<PMD_CELL_NEIGHBOUR;ii++){
00275                 Int_t id1 = i + neibx[ii];
00276                 Int_t jd1 = j + neiby[ii];
00277                 if(d1[id1][jd1]==0.){
00278                   count++;
00279                   if(count == PMD_CELL_NEIGHBOUR){//finding isolated cell
00280                     gsuper = spmhit1->Gsuper();
00281                     gsuper = gsuper + 11;
00282                     if(mOptHist)mMipEnergy[gsuper][i][j]->Fill(d1[i][j]);
00283                     Int_t channel=(Int_t)d1[i][j];
00284                     if(channel<MIP_CH_MAX)MipArray[gsuper][i][j][channel]++;;
00285                  
00286                   }
00287                 }
00288               }
00289             }
00290           }
00291         }
00292         
00293       }// if-- module hit
00294     }// for supermodule
00295   }// if---PMD detector
00296 }//void
00297 
00298 Int_t StPmdCalibConstMaker::FindMipParameters()
00299 {
00300   Float_t MipFitParam[3];
00301   Float_t MipFitChiSqr, MipPeakPosition, MipPeakWidth;
00302   TF1 LandauFunction("LandauFunction","landau",0,100);
00303   LandauFunction.SetParameters(1,1,1);
00304   LandauFunction.SetParNames("Constant","MPV","Sigma");
00305   Float_t MPV_Sum   = 0;
00306   Int_t   MPV_Count = 0;
00307   TH1F * miphist;
00308   Stat_t entryFlag=0;
00309   for(Int_t sm=0;sm<2*PMD_CRAMS_MAX;sm++){
00310     for(Int_t i =0;i<imax[sm];i++){
00311       for(Int_t j =0;j<jmax[sm];j++){
00312         if(!mOptHist){
00313          miphist=new TH1F("miphist","MipHist",200,0.,140.);
00314        for(Int_t k =0;k<10;k++){
00315               miphist->SetBinContent(k+1,(Stat_t)MipArray[sm][i][j][k]);
00316        }
00317         entryFlag = miphist->GetEntries();
00318         }
00319         if(mOptHist)entryFlag = mMipEnergy[sm][i][j]->GetEntries();
00320 
00321         //if(entryFlag > 5){
00322          if(entryFlag > MIP_MIN_ENTRY){ // new 
00323          if(mOptHist)mMipEnergy[sm][i][j]->Fit("LandauFunction","q","r");
00324          if(!mOptHist)miphist->Fit("LandauFunction","q","r");
00325           for(Int_t m=0; m < 3; m++){
00326             MipFitParam[m]    = LandauFunction.GetParameter(m);
00327           }
00328 
00329           if(MipFitParam[1] > 0.){
00330             MPV_Sum             = MPV_Sum + MipFitParam[1];
00331             MPV_Entry[sm][i][j] = MipFitParam[1]; 
00332             MPV_Count++;
00333           }
00334           MipFitChiSqr        = LandauFunction.GetChisquare();
00335           //Float_t ndf       = func->GetNDF();
00336           if(mOptHist){
00337                 MipPeakPosition     = mMipEnergy[sm][i][j]->GetMean();
00338                 MipPeakWidth        = mMipEnergy[sm][i][j]->GetRMS();
00339           }
00340           else{
00341                 MipPeakPosition     = miphist->GetMean();
00342                 MipPeakWidth        = miphist->GetRMS();
00343           }
00344           
00345           Int_t BrdNo=0;
00346           mPmdDbUtil->BoardNumber(sm,i,j,BrdNo);
00347           Int_t BrdCh=0;
00348           mPmdDbUtil->ChannelInBoard(sm,i,j,BrdCh);
00349           
00350           //mMipPeak[BrdNo][BrdCh]=MipPeakPosition;//!new
00351           mMipWidth[BrdNo][BrdCh]=MipPeakWidth;
00352           
00353         }
00354         if(!mOptHist)delete miphist;
00355       } // ! col
00356     }  
00357   } 
00358   
00359   Float_t MPV_Av = MPV_Sum/MPV_Count;
00360   //  Float_t MPV_Av = 4.5;
00361   
00362   for(Int_t sm=0;sm<2*PMD_CRAMS_MAX;sm++){
00363     for(Int_t i =0;i<imax[sm];i++){
00364       for(Int_t j =0;j<jmax[sm];j++){
00365         if(MPV_Entry[sm][i][j] > 0.)
00366           normFactor[sm][i][j] = (Float_t) MPV_Av/(Float_t) MPV_Entry[sm][i][j];
00367         Int_t BrdNo=0;
00368         mPmdDbUtil->BoardNumber(sm,i,j,BrdNo);
00369         Int_t BrdCh=0;
00370         mPmdDbUtil->ChannelInBoard(sm,i,j,BrdCh);
00371         
00372         mMipPeak[BrdNo][BrdCh]=normFactor[sm][i][j];
00373       }
00374     }
00375   }
00376   return kStOK;
00377 }
00378 
00379 Int_t StPmdCalibConstMaker::Finish()
00380 {
00381   cout << " *** I am in StPmdCalibConstMaker::Finish() " << endl;
00382   
00383   TStopwatch clock;
00384   clock.Start();
00385   
00386   FindMipParameters(); 
00387   if(mSaveCalibToDB)SaveCalibration(); 
00388   
00389   gMessMgr->Info("StPmdCalibConstMaker::Finish()");
00390   return StMaker::Finish();
00391   //return kStOK;
00392 }
00393 
00394 //_____________________________________________________________________________
00395 void StPmdCalibConstMaker::SaveCalibration()
00396 {
00397   
00398   cout <<" ***** I am now saving the Calibration to DB : Wait !!" << endl;
00399   //  TDatime *tt = new TDatime(date,time);
00400   //  TString timestamp = tt->AsSQLString();
00401   //  delete tt;
00402   
00403   //  if(!mSaveCalibToDB) return;
00404   /* Presently it is storing some numbers , coming from simulated files,
00405      it is tested with daq file, and time stamp works., so once we have good 
00406      amount fo data , it should work*/
00407   
00408   
00409   StDbManager* mgr=StDbManager::Instance();
00410   StDbConfigNode* node=mgr->initConfig("Calibrations_pmd");
00411   StDbTable* tab1=node->addDbTable("pmdBrdMipCalib");
00412   StDbTable* tab2=node->addDbTable("pmdCalSummary");
00413   pmdBrdMipCalib_st pbd[PMD_BOARD_MAX];
00414   memset(pbd,0,PMD_BOARD_MAX*sizeof(pmdBrdMipCalib_st));
00415   for(Int_t ism=0;ism<2*PMD_CRAMS_MAX;ism++){
00416     for(Int_t irow=0;irow<jmax[ism];irow++){
00417       for(Int_t icol=0;icol<imax[ism];icol++){
00418         Int_t brd=0;
00419         Int_t brd_ch=0;
00420         mPmdDbUtil->BoardNumber(ism,irow,icol,brd);
00421         mPmdDbUtil->ChannelInBoard(ism,irow,icol,brd_ch);
00422         pbd[brd].FEEBoardNumber = brd;
00423         //filling peak and width from arrays
00424         
00425         pbd[brd].MipPeakPosition[brd_ch] = mMipPeak[brd][brd_ch];
00426         pbd[brd].MipPeakWidth[brd_ch] = mMipWidth[brd][brd_ch];
00427       }
00428     }
00429   }
00430   tab1->SetTable((char*)pbd,PMD_BOARD_MAX);
00431   
00432   pmdCalSummary_st pcs[1];
00433   memset(pcs,0,sizeof(pmdCalSummary_st));
00434   // Summary table will be filled here. 
00435   
00436   tab2->SetTable((char*)pcs,1);
00437   
00438   //storing timestamp from the first event in the chain
00439   char timestamp[15]; 
00440   sprintf(timestamp,"%08d%06d",mDate,mTime);
00441   mgr->setStoreTime(timestamp);
00442   //mgr->setStoreTime("2003-12-29 11:11:15");
00443   cout<<"PMD Calib: storing in DB"<<endl;
00444   
00445   mgr->storeDbTable(tab1);
00446   mgr->storeDbTable(tab2);
00447   cout<<"PMD Calib: Stored in DB "<<endl;
00448   
00449 }
00450 
00451 
00452 void StPmdCalibConstMaker::InitMipParams()
00453 {
00454   for(Int_t i=0;i<PMD_BOARD_MAX;i++){
00455     for(Int_t j=0;j<PMD_BOARD_CH_MAX;j++){
00456       mMipPeak[i][j]=0.;
00457       mMipWidth[i][j]=0.;
00458       
00459     }
00460   }
00461 }
00462 
00463 void StPmdCalibConstMaker::ClearHists()
00464 {
00465   cout << " *** Let me to clear the Histograms in PmdCalib" << endl;
00466   
00467   for(Int_t sm=0;sm<2*PMD_CRAMS_MAX;sm++){
00468     for(Int_t i =0;i<imax[sm];i++){
00469       for(Int_t j =0;j<jmax[sm];j++){
00470         
00471         if(mMipEnergy[sm][i][j]){
00472           mMipEnergy[sm][i][j]=0;
00473           delete mMipEnergy[sm][i][j];
00474         }
00475       }
00476     }
00477   }
00478 }
00479 void StPmdCalibConstMaker::ClearMipArray()
00480 {
00481   for(Int_t sm=0;sm<2*PMD_CRAMS_MAX;sm++){
00482     for(Int_t i =0;i<imax[sm];i++){
00483       for(Int_t j =0;j<jmax[sm];j++){
00484        for(Int_t k =0;k<MIP_CH_MAX;k++){
00485           MipArray[sm][i][j][k]=0;
00486         }
00487       }
00488     }
00489   }
00490 }
00491 
00492 
00493 
00494 
00495 
00496 
00497 
00498 
00499 
00500 
00501