StRoot: StEmcPool/StBarrelCalib2008/StBarrelMonitorMaker.cxx Source File

00001 // *-- Author : Jan Balewski
00002 // 
00003 // $Id: StBarrelMonitorMaker.cxx,v 1.4 2009/08/25 16:17:48 fine Exp $
00004 
00005 #include <TFile.h>
00006 #include <TH1.h>
00007 #include <TH2.h>
00008 #include <TH3.h>
00009 #include <StMessMgr.h>
00010  
00011 #include "StBarrelMonitorMaker.h"
00012 #include "StJanBarrelDbMaker.h"
00013 #include "JanBprsEveA.h"
00014 #include "BprsCapPolygraph.h"
00015 #include "BarrelMipCalib.h"
00016 
00017 #include "StMuDSTMaker/COMMON/StMuEvent.h"
00018 #include "StMuDSTMaker/COMMON/StMuDst.h"
00019 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
00020 
00021 #include "StEmcRawMaker/defines.h"
00022 #include "StEmcUtil/geometry/StEmcGeom.h"
00023 
00024 #include "StDetectorDbMaker/St_tpcGasC.h"
00025 #include "St_db_Maker/St_db_Maker.h"
00026 #include "StEmcUtil/database/StEmcDecoder.h"
00027 
00028 
00029 #include "StEventTypes.h"
00030 #include "StMcEventTypes.hh"
00031 #include "StMcEvent.hh"
00032 
00033 
00034 ClassImp(StBarrelMonitorMaker)
00035 
00036 //________________________________________________
00037 //________________________________________________
00038 StBarrelMonitorMaker::StBarrelMonitorMaker( const char* self ): StMaker(self){
00039 
00040   cTile[0]='T';  cTile[1]='P';
00041   cTile4[0]="BTOW";  cTile4[1]="BPRS";
00042 
00043   nInpEve=nAcceptEve=nTrigEve=nCorrEve= mGeantEveInp=0;
00044  
00045   setHList(0);
00046   setTrigIdFilter(0);
00047   setMC(0);
00048 
00049   mMappB=0;
00050   hBprs3D=0;// special very large 3D histo for BPRS cap monitoring
00051 
00052   setBprsHisto(1);
00053   setCalibPass(0);
00054 }
00055  
00056 //________________________________________________
00057 //________________________________________________
00058 Int_t 
00059 StBarrelMonitorMaker::Init(){
00060 
00061   assert(HList);
00062   initHistos();
00063 
00064   mBtowGeom = StEmcGeom::instance("bemc");
00065   mBprsGeom = StEmcGeom::instance("bprs");
00066   mSmdEGeom = StEmcGeom::instance("bsmde");
00067   mSmdPGeom = StEmcGeom::instance("bsmdp"); // needed only to map m-s-e --> softID
00068 
00069   mMuDstMaker = (StMuDstMaker*)GetMaker("MuDst"); assert(mMuDstMaker);
00070   mJanDbMaker = (StJanBarrelDbMaker*)GetMaker("janBarrelDb"); assert(mJanDbMaker);
00071 
00072   bprsPolygraph=new BprsCapPolygraph(HList, mJanDbMaker) ;
00073   bprsPolygraph->setCut(50, 0.98,1,0.1);
00074   bprsPolygraph->print(); 
00075 
00076   mipCalib=new BarrelMipCalib(HList, mJanDbMaker,  mMuDstMaker) ;
00077   mipCalib->setCut(60.,0.35,1.3,0.51, 3.3, 1.0, 180); 
00078   //               zV, pT, eta, nFit, dedx, zMargin, Rxy
00079   mipCalib->print(); 
00080 
00081   LOG_INFO<<Form("::Init() filter trigID=%d done", trigID)<<endm;  
00082   return StMaker::Init();
00083 }
00084 
00085 //________________________________________________
00086 //________________________________________________
00087 Int_t 
00088 StBarrelMonitorMaker::InitRun(int runNo){
00089   LOG_INFO<<Form("::InitRun(%d) start",runNo)<<endm;
00090   St_db_Maker* mydb = (St_db_Maker*) StMaker::GetChain()->GetMaker("StarDb");
00091   assert(mydb);
00092    
00093   // this is how BTOW mapping is accesible
00094   mMappB = new StEmcDecoder(mydb->GetDateTime().GetDate(),mydb->GetDateTime().GetTime());
00095   float airPres=St_tpcGasC::instance()-> barometricPressure();
00096   LOG_INFO<<Form("::InitRun(%d) AirPressure=%.2f",runNo, airPres)<<endm;
00097   if(runNo==1000000) {
00098     LOG_WARN<<Form("::InitRun(%d) ??? , it is OK for M-C ",runNo)<<endm;
00099   }  
00100 
00101   LOG_INFO<<Form("::InitRun() algo params:  bprsHisto=%d calibPass=%d (bits)",
00102                  par_bprsHisto,par_calibPass
00103                  )<<endm;
00104   return kStOK;
00105 }
00106 
00107 //________________________________________________
00108 //________________________________________________
00109 Int_t 
00110 StBarrelMonitorMaker::Finish(){
00111  
00112   gMessMgr->Message("","I") <<GetName()<<"::Finish()\n    inputEve="<<nInpEve<<" trigFilterEve="<<nTrigEve<<" nCorrEve="<<nCorrEve<<" nAcceptEve="<<nAcceptEve<<endm;
00113   return kStOK;
00114 }
00115 
00116 //________________________________________________
00117 //________________________________________________
00118 void 
00119 StBarrelMonitorMaker::Clear(const Option_t*){
00120   eventID=-999;
00121   janEve.clear();
00122   for( int icr=0; icr <mxBprsCrate;icr++) janBprsEveA[icr].clear();
00123 }
00124 
00125 //________________________________________________
00126 //________________________________________________
00127 //________________________________________________
00128 Int_t 
00129 StBarrelMonitorMaker::Make(){
00130   nInpEve++;
00131   StMuEvent* muEve = mMuDstMaker->muDst()->event();
00132 
00133   eventID=muEve->eventId();
00134   janEve.id=eventID;
00135 
00136   int nPrimV=mMuDstMaker->muDst()->numberOfPrimaryVertices();
00137  
00138   LOG_INFO <<GetName()<<"\n\n::Make()================ isMC="<<isMC<<" eveID="<<eventID<<" nInpEve="<<nInpEve<<" nPrimV="<<nPrimV<<endm;
00139   
00140 
00141   vector<unsigned int> trgL=mMuDstMaker->muDst()->event()->triggerIdCollection().nominal().triggerIds();
00142   printf("trigL len=%d\n",int(trgL.size()));
00143   uint ii;
00144   for( ii=0;ii<trgL.size();ii++) printf("ii=%d trigID=%d\n",ii,trgL[ii]);
00145   if (mMuDstMaker->muDst()->event()->triggerIdCollection().nominal().isTrigger(19)) return kStOK; // drop bbc-fast events
00146   
00147 
00148   //  printf("trigOK\n"); 
00149   // nTrigEve++;
00150 
00151 
00152   unpackStTiles(kBTow);
00153   unpackStTiles(kBPrs);  
00154 
00155   populateBprsEveA();
00156   
00157   // ...... check for BPRS capID corruption .....
00158   for(int bprsCrateID=0; bprsCrateID<mxBprsCrate;  bprsCrateID++) 
00159   // int bprsCrateID=2;
00160    {    
00161     JanBprsEveA &bprsEve=janBprsEveA[bprsCrateID];
00162     bprsPolygraph->doBaseline(bprsEve,janEve);
00163     bprsPolygraph->findBestCap(bprsEve,janEve);
00164     bprsPolygraph->doPedResidua(bprsEve); 
00165     if(par_calibPass & kPassCapFix)  // do correct  capID corruption
00166       janEve.bprsCap[bprsCrateID]=bprsEve.getBestCapID();   
00167   }
00168   //if(nInpEve<3)  janEve.print(2);  
00169   calibrateTiles(kBPrs);
00170   calibrateTiles(kBTow);  
00171  
00172 
00173  // 
00174 
00175   mipCalib->search(janEve); // full analysis
00176   // mipCalib->searchEtaBin20(janEve); // just last eta bin location
00177 
00178   return kStOK;
00179 } 
00180 
00181 
00182 
00183 //___________________ _____________________________
00184 //________________________________________________
00185 StBarrelMonitorMaker::~StBarrelMonitorMaker(){
00186   
00187 }
00188 
00189 //___________________ _____________________________
00190 //________________________________________________
00191 void 
00192 StBarrelMonitorMaker::saveHisto(TString fname){
00193   TString outName=fname+".hist.root";
00194   TFile f( outName,"recreate");
00195   assert(f.IsOpen());
00196   printf("HHH %d histos are written  to '%s' ...\n",HList->GetEntries(),outName.Data());
00197 
00198   HList->Write();
00199   f.Close();
00200 
00201 }
00202 
00203 //________________________________________________
00204 //________________________________________________
00205 void
00206 StBarrelMonitorMaker::unpackStTiles(int ibp){
00207  StEvent *mEvent = (StEvent*)StMaker::GetChain()-> GetInputDS("StEvent");
00208   assert(mEvent);
00209 
00210   //.........................  BTOW  or  BPRS  ....................
00211   //use StEvent as default  to get simulation right in BEMC
00212   StEmcCollection *emc = mEvent->emcCollection();
00213   assert (emc);
00214 
00215   int jBP=BTOW; 
00216   StEmcDetector* detector=emc->detector(kBarrelEmcTowerId);
00217   StEmcGeom  *geomB=mBtowGeom;
00218   if( ibp==kBPrs) {// change pointers for  preshower
00219     jBP=BPRS; 
00220     geomB=mBprsGeom;
00221     detector=emc->detector(kBarrelEmcPreShowerId);
00222   }
00223 
00224   if(!detector) { 
00225     printf("no %s  data, nInpEve=%d\n",cTile4[ibp],nInpEve);
00226     return;
00227   }
00228 
00229   // for BPRS unpack caps and store locally
00230   if(ibp==kBPrs ) {
00231     for(int icr=0;icr<mxBprsCrate;icr++) {
00232       StEmcModule* module = detector->module(2+icr*30);assert(module);
00233       StSPtrVecEmcRawHit& rawHit=module->hits();
00234       if(rawHit.size()<=0) {
00235          printf("ss icr=%d n=%dL, ABORT BpRS for this event\n",icr,int(rawHit.size())); 
00236         return ;
00237 }       
00238         assert(rawHit.size()>0);
00239         janEve.bprsCap[icr]= rawHit[10]->calibrationType();
00240       // tmp
00241       if(icr==3)    janEve.bprsCap[icr]= rawHit[32]->calibrationType();
00242 
00243     }
00244   }
00245 
00246   // unpack raw data and store it locally in an array
00247   for(Int_t m = 1; m <= 120; ++m) {
00248     StEmcModule* module = detector->module(m);
00249     //printf("MMM %s m=%p mod=%d\n",cTile4[ibp],module,m);
00250     if(!module) continue;
00251     
00252     StSPtrVecEmcRawHit& rawHit=module->hits();
00253     // printf("HHH n=%d\n",rawHit.size());
00254     for(UInt_t k = 0; k < rawHit.size(); ++k){
00255       int id;
00256       
00257       Int_t m=rawHit[k]->module();
00258       Int_t e=rawHit[k]->eta();
00259       Int_t s=abs(rawHit[k]->sub());
00260       float rawAdc=rawHit[k]->adc()+0.5;// return middle of the bin
00261       
00262       //Get software tower id to get DaqID
00263       geomB->getId(m,e,s,id);
00264 
00265 #if 0       // tmp for Matt, to QA new mapping
00266       if(ibp==kBPrs ){
00267         int   id1=(int)mJanDbMaker->bprsReMap()->GetBinContent(id);
00268         printf("BPRSX %d %d %d %d\n",nInpEve-1,id1,rawHit[k]->adc(),rawHit[k]->calibrationType());
00269       }
00270       // end tmp
00271 #endif
00272 
00273       assert(id>=1);
00274       assert(id<=mxBtow);
00275       janEve.rawAdcTile[ibp][id-1]=rawAdc;
00276       // if(id%20==7) 
00277       // if(ibp==kBPrs && (id>650 && id<710)) printf("unpack %s  hit=%d softid=%d, m=%d rawADC=%.1f\n",cTile4[ibp],k,id,m,rawAdc);
00278     }
00279   }// end of module loop
00280 
00281   janEve.tileIn[ibp]=1; // tag data block as delivered
00282 
00283   assert(!isMC); // verify peds for MC
00284 }
00285 //________________________________________________
00286 //________________________________________________
00287 void
00288 StBarrelMonitorMaker::calibrateTiles(int ibp){
00289   // printf("ZZ %d %d\n",ibp, janEve.tileIn[ibp]);
00290   if(janEve.tileIn[ibp]==0) return; // no data in this event
00291 
00292   for(int id0=0; id0<mxBtow; id0++) {
00293     int   id=id0+1;
00294     int   crateID=0;
00295 
00296     int   capID=0;
00297     int   stat =mJanDbMaker->statTile(ibp,id);
00298 
00299     if(ibp==kBPrs ) {
00300      crateID= mJanDbMaker->bprsCrate(id);
00301      capID=janEve.bprsCap[crateID]; //fix capID corruption in 2008 
00302     }
00303     
00304     float ped     =mJanDbMaker->pedTile(ibp,id,capID);
00305     float rawAdc  =janEve.rawAdcTile[ibp][id0];
00306 
00307     // if(ibp==kBPrs) printf("ss id=%d stat=%d  rawAdc=%.1f\n", id,stat,rawAdc);
00308     if(ibp==kBPrs  && rawAdc>100 && rawAdc<300){
00309       //accumulate pedestals vs. capID similat to Tonko's on-line approach
00310       hTonko0->Fill(id,capID,1.);
00311       hTonko1->Fill(id,capID,rawAdc);
00312       hTonko2->Fill(id,capID,rawAdc*rawAdc);
00313     }
00314 
00315     if(!(par_calibPass & kPassPedSub)) {
00316       hTile[ibp]->Fill(id,rawAdc);
00317       if(par_bprsHisto==2) hBprs3D->Fill(id,rawAdc,capID);
00318     }
00319     
00320     if(par_calibPass==0) continue; // no swap info loaded in this mode
00321 
00322     //------ TMP ----- apply BPRS & BTOW swap only for ped-corrected ADC data
00323     if(ibp==kBTow) {
00324       id=(int)mJanDbMaker->btowReMap()->GetBinContent(id);
00325     } else  if(ibp==kBPrs) {
00326       id=(int)mJanDbMaker->bprsReMap()->GetBinContent(id);
00327     } else assert(1==2);
00328 
00329     // below  use only 'id' for indexing to pick up this swaping
00330 
00331 
00332     if(stat) {
00333       janEve.statTile[ibp][id-1]=stat; // is bad
00334       continue;
00335     } else {
00336       janEve.statTile[ibp][id-1]=0; // is good
00337     }
00338 
00339     float adc=rawAdc-ped;    
00340     janEve.adcTile[ibp][id-1]= adc;
00341     
00342     if(par_calibPass & kPassPedSub){
00343       hTile[ibp]->Fill(id,adc);
00344       if(ibp==kBPrs){ // later add alos BTOW histos per crate, fix it
00345         hBprsA[crateID]->Fill(adc);
00346         if(par_bprsHisto==2) hBprs3D->Fill(id,adc,capID);
00347       }
00348     }
00349     
00350     //..... energy calibration goes here, later
00351     //...
00352 
00353   } // end of data loop 
00354 
00355 
00356 }
00357 
00358 //________________________________________________
00359 //________________________________________________
00360 void
00361 StBarrelMonitorMaker::populateBprsEveA(){
00362   
00363   if(janEve.tileIn[kBPrs]==0) return; // no data in this event
00364 
00365   for( int icr=0; icr <mxBprsCrate;icr++) {
00366     int   capID=janEve.bprsCap[icr];
00367     janBprsEveA[icr].set(capID,icr,eventID);
00368   }
00369          
00370   int ibp=kBPrs;
00371   for(int id0=0; id0<mxBtow; id0++) {
00372     int id=id0+1;
00373     int   icr=mJanDbMaker->bprsCrate(id);
00374     float rawAdc  =janEve.rawAdcTile[ibp][id0];
00375     janBprsEveA[icr].addRawValue(id,rawAdc);
00376     //printf("id=%d, icr=%d \n",id,icr);
00377   }
00378   
00379 }
00380     
00381 
00382 //---------------------------------------------------
00383 // $Log: StBarrelMonitorMaker.cxx,v $
00384 // Revision 1.4  2009/08/25 16:17:48  fine
00385 // fix the compilation issues under SL5_64_bits  gcc 4.3.2
00386 //
00387 // Revision 1.3  2009/08/25 16:08:04  fine
00388 // fix the compilation issues under SL5_64_bits  gcc 4.3.2
00389 //
00390 // Revision 1.2  2009/02/04 20:33:32  ogrebeny
00391 // Moved the EEMC database functionality from StEEmcDbMaker to StEEmcUtil/database. See ticket http://www.star.bnl.gov/rt2/Ticket/Display.html?id=1388
00392 //
00393 // Revision 1.1  2008/11/24 23:06:35  balewski
00394 // start
00395 //