StRoot: StEEmcPool/StEzSmdCal/StEEsmdCalMaker.cxx Source File

00001 // *-- Author : Jan Balewski
00002 // 
00003 // $Id: StEEsmdCalMaker.cxx,v 1.7 2009/02/04 20:33:23 ogrebeny Exp $
00004 
00005 #include <TFile.h>
00006 #include <TH2.h>
00007 
00008 #include "StEEsmdCalMaker.h"
00009 
00010 #include "StMuDSTMaker/COMMON/StMuEvent.h"
00011 #include "StMuDSTMaker/COMMON/StMuDst.h"
00012 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
00013 
00014 #include <StMessMgr.h>
00015 
00016 #include "StEEmcUtil/database/EEmcDbItem.h"
00017 #include "StEEmcUtil/database/StEEmcDb.h"
00018 
00019 
00020 ClassImp(StEEsmdCalMaker)
00021 
00022 //________________________________________________
00023 //________________________________________________
00024 StEEsmdCalMaker::StEEsmdCalMaker( const char* self ,const char* muDstMakerName) : StMaker(self){
00025   mMuDstMaker = (StMuDstMaker*)GetMaker(muDstMakerName);
00026   assert(mMuDstMaker);
00027   MCflag=0;// default off
00028   //  printf("constr of calib =%s=\n",GetName());
00029  
00030   //tmp, JB
00031   // ideal calibration used by Fast simulator has eta dependent gains for towers and plain energy deposit for SMD,pre,post
00032   
00033   // towers are gain matched to fixed E_T
00034   float maxAdc=4095;
00035   float maxEtot=60;  // in GeV
00036   const float feta[MaxEtaBins]= {1.95,1.855,1.765,1.675,1.59,1.51,1.435,1.365,1.3,1.235,1.17,1.115};
00037   
00038   int i;
00039   
00040   mfixEmTgain=new float [MaxEtaBins];
00041   for (i=0;i<MaxEtaBins;i++) {
00042     mfixEmTgain[i]=maxAdc/maxEtot/cosh(feta[i]); 
00043   }
00044   mfixSMDgain=23000;
00045   mfixPgain=23000;
00046   
00047 }
00048 
00049 
00050 //___________________ _____________________________
00051 //________________________________________________
00052 StEEsmdCalMaker::~StEEsmdCalMaker(){
00053 }
00054 
00055 //________________________________________________
00056 //________________________________________________
00057 void StEEsmdCalMaker::SetSector(int sec){
00058   setSector(sec);
00059   TString name=GetName();
00060   name+="-";
00061   name+=sec;
00062   //  printf("change name to %s sec=%d\n", name.Data(),sec);
00063   SetName(name);
00064 }
00065 
00066  
00067 //________________________________________________
00068 //________________________________________________
00069 Int_t StEEsmdCalMaker::Init(){
00070   eeDb = (StEEmcDb*)this->GetDataSet("StEEmcDb");
00071   //printf("%s got eeDb=%p\n",GetName(),eeDb);
00072   EEsmdCal::init();
00073   printf("%s has MCflag=%d\n",GetName(),MCflag);
00074   return StMaker::Init();
00075 }
00076 
00077 //________________________________________________
00078 //________________________________________________
00079 Int_t StEEsmdCalMaker::InitRun(int runNo){
00080   if(runNo==0) {
00081     gMessMgr->Message("","W")<<GetName()<<"::InitRun("<<runNo<<") ??? changed to 555, it s OK for M-C - perhaps, JB"<<endm;
00082     runNo=555;
00083   }
00084   initRun(runNo);
00085   return kStOK;
00086 }
00087 
00088 
00089 
00090 //________________________________________________
00091 //________________________________________________
00092 Int_t StEEsmdCalMaker::Finish(){
00093   finish(0); // do not draw
00094   return kStOK;
00095 }
00096 
00097 
00098 //________________________________________________
00099 //________________________________________________
00100 Int_t StEEsmdCalMaker::Make(){
00101   
00102   clear();
00103   // ............. acuire EEMC data 
00104   if(unpackMuDst()<0)   return kStOK;
00105 
00106   findSectorMip();// do real analysis
00107  
00108   return kStOK;
00109 }
00110 
00111 
00112 //________________________________________________
00113 //________________________________________________
00114 Int_t StEEsmdCalMaker::unpackMuDst(){
00115   
00116   nInpEve++;
00117   gMessMgr->Message("","D") <<GetName()<<"::::getAdc() is called "<<endm;
00118   
00119   // Access to muDst .......................
00120   StMuEmcCollection* emc = mMuDstMaker->muDst()->muEmcCollection();
00121   if (!emc) {
00122     gMessMgr->Message("","W") <<"No EMC data for this event"<<endm;    return kStOK;
00123   }
00124     
00125 
00126   int i, n1=0,n2=0,n3=0;
00127 
00128   //.........................  T O W E R S .....................
00129   for (i=0; i< emc->getNEndcapTowerADC(); i++) {
00130     int sec,eta,sub,val;
00131     //muDst  ranges:sec:1-12, sub:1-5, eta:1-12
00132     emc->getEndcapTowerADC(i,val,sec,sub,eta);
00133     assert(sec>0 && sec<=MaxSectors);// total corruption of muDst
00134 
00135     //tmp, for fasted analysis use only hits from sectors init in DB
00136     if(sec<eeDb->getFirstSector() || sec>eeDb->getLastSector()) continue;
00137  
00138     const EEmcDbItem *x=eeDb->getTile(sec,'A'+sub-1,eta,'T');
00139     assert(x); // it should never happened for muDst
00140     if(x->fail ) continue; // drop broken channels
00141 
00142     int iphi=(x->sec-1)*MaxSubSec+(x->sub-'A');
00143     int ieta=x->eta-1;
00144     assert(iphi>=0 && iphi<MaxPhiBins);
00145     assert(ieta>=0 && ieta<MaxEtaBins);
00146  
00147     float adc=-100, rawAdc=-101, ene=-102;
00148     
00149     if(MCflag) {  // M-C & real data needs different handling
00150       adc=val; //  this is stored in muDst
00151       rawAdc=adc+x->ped; // ped noise could be added
00152       ene=adc/ mfixEmTgain[ieta];
00153     } else {
00154       rawAdc=val;
00155       adc=rawAdc-x->ped;
00156       if(x->gain) ene=adc/x->gain;
00157     }
00158     
00159     int iT=0;// for towers
00160     
00161     // if(adc>0) printf("adc=%f ene/GeV=%f rawAdc=%f idG=%f,hSec=%d %s\n",adc,ene,rawAdc,mfixEmTgain[ieta],sectID,x->name); 
00162    
00163     tileAdc[iT][ieta][iphi]=adc;// store towers
00164     tileThr[iT][ieta][iphi]=rawAdc>x->thr;
00165     if(rawAdc>x->thr)  n1++;
00166     if(x->gain<=0) continue;// drop channels w/o gains
00167     tileEne[iT][ieta][iphi]=ene;
00168     
00169   }
00170 
00171 
00172   //.........................  P R E - P O S T .....................  
00173   int pNh= emc->getNEndcapPrsHits();
00174   for (i=0; i<pNh; i++) {
00175     int pre;
00176     int sec,eta,sub;
00177     //muDst  ranges: sec:1-12, sub:1-5, eta:1-12 ,pre:1-3==>pre1/pre2/post
00178     StMuEmcHit *hit=emc->getEndcapPrsHit(i,sec,sub,eta,pre);
00179 
00180     //tmp, for fasted analysis use only hits from sectors init in DB
00181     if(sec<eeDb->getFirstSector() || sec>eeDb->getLastSector()) continue;
00182  
00183     
00184     //Db ranges: sec=1-12,sub=A-E,eta=1-12,type=T,P-R ; slow method
00185     const EEmcDbItem *x=eeDb-> getTile(sec,sub-1+'A', eta, pre-1+'P'); 
00186     if(x==0) continue;
00187     if(x->fail ) continue; // drop broken channels
00188     
00189     // accept this hit
00190     int iphi=(x->sec-1)*MaxSubSec+(x->sub-'A');
00191     int ieta=x->eta-1;
00192     assert(iphi>=0 && iphi<MaxPhiBins);
00193     assert(ieta>=0 && ieta<MaxEtaBins);
00194     int iT=pre; // P,Q,R fall in to iT=1,2,3 <== there is no '+1' error, JB
00195     assert(iT>0 && iT<mxTile);
00196     
00197     float val=hit->getAdc();
00198     float adc=-100, rawAdc=-101, ene=-102;
00199   
00200     if(MCflag) {// val is Geant energy * const
00201       adc=val;
00202       rawAdc=adc+x->ped; // ped noise could be added
00203       ene=val/mfixPgain; // (GeV) Geant energy deposit 
00204     } else {
00205       rawAdc=val;
00206       adc=rawAdc-x->ped;
00207       if(x->gain) ene=adc/x->gain;
00208     }
00209     if(rawAdc>x->thr)  n2++;
00210     //if(adc>0) printf("adc=%f ene/GeV=%f  hSec=%d %s g=%f\n",adc,ene,sectID,x->name, x->gain); 
00211       
00212     tileAdc[iT][ieta][iphi]=adc;// // store P,Q,R depending on 'iT'
00213     tileThr[iT][ieta][iphi]=rawAdc>x->thr;
00214     
00215     if(x->gain<=0) continue;// drop channels w/o gains
00216     tileEne[iT][ieta][iphi]=ene;
00217 
00218   }
00219   
00220 
00221   //.......................  S M D ................................
00222   
00223   char uv='U';
00224   for(uv='U'; uv<='V'; uv++) {
00225     int sec,strip;
00226     int nh= emc->getNEndcapSmdHits(uv);
00227     //printf("pl=%c nH=%d  secLim=%d %d\n",uv,nh,eeDb->mfirstSecID,eeDb->mlastSecID);
00228     for (i=0; i<nh; i++) {
00229       StMuEmcHit *hit=emc->getEndcapSmdHit(uv,i,sec,strip);
00230       assert(sec>0 && sec<=MaxSectors);// total corruption of muDst
00231       
00232       //tmp, for fasted analysis use only hits from sectors init in DB
00233       if(sec<eeDb->getFirstSector() || sec>eeDb->getLastSector()) continue;
00234       
00235       const EEmcDbItem *x=eeDb->getByStrip(sec,uv,strip);
00236       assert(x); // it should never happened for muDst
00237       if(x->fail ) continue; // drop broken channels
00238 
00239       float val=hit->getAdc();
00240       float adc=-100, rawAdc=-101, ene=-102;
00241 
00242       // M-C & real data needs different handling
00243       if(MCflag) {
00244         adc=val;
00245         rawAdc=adc+x->ped; // ped noise could be added
00246         ene=val/mfixSMDgain; // (GeV) Geant energy deposit      
00247       }else {
00248         rawAdc=val;
00249         adc=rawAdc-x->ped;
00250         if(x->gain) ene=adc/x->gain;
00251       }
00252       if(rawAdc>x->thr)  n3++;
00253       //   if(adc>0)    printf("adc=%f ene/GeV=%f  hSec=%d %s\n",adc,ene,sectID,x->name); 
00254 
00255       smdAdc[x->plane-'U'][x->strip-1]=adc;
00256       if(x->gain<=0)continue; // drop channels w/o gains
00257       smdEne[x->plane-'U'][x->strip-1]=ene;
00258     }
00259   }
00260   //  printf("nTw=%d nPQR=%d,nSmd=%d\n",n1,n2,n3);
00261   
00262   return n1;
00263 }
00264 
00265 
00266 // $Log: StEEsmdCalMaker.cxx,v $
00267 // Revision 1.7  2009/02/04 20:33:23  ogrebeny
00268 // Moved the EEMC database functionality from StEEmcDbMaker to StEEmcUtil/database. See ticket http://www.star.bnl.gov/rt2/Ticket/Display.html?id=1388
00269 //
00270 // Revision 1.6  2005/03/11 15:44:25  balewski
00271 // works with muEzt, cucu200
00272 //
00273 // Revision 1.5  2004/10/21 13:31:31  balewski
00274 // to match new name of emcCollection in muDst
00275 //
00276 // Revision 1.4  2004/07/27 21:59:47  balewski
00277 // now runs on muDst as well
00278 //
00279 // Revision 1.3  2004/06/29 16:37:41  balewski
00280 // towards SMD calib
00281 //
00282 // Revision 1.2  2004/06/22 23:31:11  balewski
00283 // few more gadgets added
00284 //
00285 // Revision 1.1  2004/06/12 04:09:25  balewski
00286 // start
00287 //
00288 // Revision 1.1  2004/06/06 04:54:08  balewski
00289 // dual analyzis
00290 //