StRoot: StEmcCalibrationMaker/StEmcCalibrationMaker.cxx Source File

00001 #include "StEmcCalibrationMaker.h"
00002 
00003 #include <math.h>
00004 
00005 #include <TFile.h>
00006 
00007 #include "StEventTypes.h"
00008 #include "StEvent.h"
00009 #include "Stiostream.h"
00010 #include "Stsstream.h"
00011 #include "StEmcUtil/geometry/StEmcGeom.h"
00012 #include "StEmcUtil/filters/StEmcFilter.h"
00013 #include "StEmcUtil/projection/StEmcPosition.h"
00014 #include "StarClassLibrary/StThreeVectorF.hh"
00015 #include "StarClassLibrary/StElectron.hh"
00016 #include "tables/St_emcStatus_Table.h" 
00017 #include "PhysicalConstants.h"
00018 #include "tables/St_emcPed_Table.h"
00019 #include "tables/St_smdPed_Table.h"
00020 #include "tables/St_emcStatus_Table.h"
00021 #include "tables/St_smdStatus_Table.h"
00022 #include "StEmcADCtoEMaker/StBemcData.h"
00023 #include "StDbLib/StDbManager.hh"
00024 #include "StDbLib/StDbTable.h"
00025 #include "StDbLib/StDbConfigNode.hh"
00026 #include "StDaqLib/EMC/StEmcDecoder.h"
00027 
00028 #define CAP1 124
00029 #define CAP2 125
00030 
00031 ClassImp(StEmcCalibrationMaker)
00032 
00033 //_____________________________________________________________________________
00034 StEmcCalibrationMaker::StEmcCalibrationMaker(const char *name):StMaker(name)
00035 {
00036   this->Clear();
00037   TString n[MAXBEMC]={"bemc","bprs","bsmde","bsmdp"};
00038   for(int i=0;i<MAXBEMC;i++) mGeom[i] = StEmcGeom::instance(n[i].Data());
00039   mPosition = new StEmcPosition();
00040   mNChannel[0] = 4800;
00041   mNChannel[1] = 4800;
00042   mNChannel[2] = 18000;
00043   mNChannel[3] = 18000;
00044 }
00045 //_____________________________________________________________________________
00046 StEmcCalibrationMaker::~StEmcCalibrationMaker()
00047 {
00048 }
00049 //_____________________________________________________________________________
00050 Int_t StEmcCalibrationMaker::Init()
00051 {
00052   return StMaker::Init();
00053 }
00054 //_____________________________________________________________________________
00055 void StEmcCalibrationMaker::zeroAll()              
00056 {
00057   mStEvent = NULL;
00058   mEmcCol = NULL;
00059   mField = 0;
00060   mL3Tracks = false;
00061   mVx = mVy = mVz = 0;
00062   mNTracks = 0;
00063   for(int i=0;i<MAXTRACK;i++)
00064   {
00065     mTrack[i] = NULL;
00066     mTrackP[i] = mTrackPt[i] = mTrackTower[i] = mTrackNPoints[i] = mTrackTowerExit[i] =0;
00067   }
00068   for(int i=0;i<MAXBEMC;i++) for(int j=0;j<mNChannel[i];j++)
00069   {
00070     mADC[i][j] = mADCPedSub[i][j] = 0;
00071     mCap[i][j] = mNTracksTower[j] = 0;
00072     mPedRms[i][j] = mPed[i][j] =  0;
00073     if(i==0) mIsIsolatedTower[j] = true;
00074                 mHasDetector[i] = false;
00075   }
00076         mZDCSum = 0;
00077         mCTBSum = 0;
00078 }
00079 //_____________________________________________________________________________
00080 Int_t StEmcCalibrationMaker::Make()
00081 {  
00082   zeroAll();
00083         mStEvent = (StEvent*)GetInputDS("StEvent");
00084         if(!mStEvent) return kStOk;
00085   mEmcCol=(StEmcCollection*)mStEvent->emcCollection();  
00086   if(!mEmcCol) return kStOk;
00087   
00088   StPrimaryVertex *v = mStEvent->primaryVertex();  
00089   if(v)
00090   {
00091     StThreeVectorF position = v->position();
00092     mVx-=position.x(); mVy=position.y(); mVz=position.z();
00093   }
00094     
00095         mIsTriggerOk = true;
00096                 
00097         //need to set BField ////////////////
00098   StEventSummary *evtSummary = mStEvent->summary();
00099   if (evtSummary) mField = evtSummary->magneticField()/10;
00100   
00101   fillEmcArrays();
00102   fillTrackArrays();
00103   
00104   return kStOK;
00105 }
00106 //_____________________________________________________________________________
00107 Int_t StEmcCalibrationMaker::Finish()
00108 {
00109   return StMaker::Finish();
00110 }
00111 //_____________________________________________________________________________
00112 void StEmcCalibrationMaker::fillEmcArrays()              
00113 {    
00114   TString n[MAXBEMC]={"bemc","bprs","bsmde","bsmdp"};
00115   
00116   for(int i = 0;i<MAXBEMC;i++)
00117   {
00118     mHasDetector[i] = false;
00119                 TString calibDb="Calibrations/emc/y3";
00120     calibDb+=n[i];
00121     TDataSet *emcDb=GetInputDB(calibDb.Data());
00122     TString table=n[i];
00123     table+="Ped";
00124     if(i<2 && emcDb) // tower and preShower
00125     {
00126       St_emcPed *ped = (St_emcPed*)emcDb->Find(table.Data());
00127       emcPed_st *pedst=NULL;
00128       if(ped) pedst=ped->GetTable();
00129       if(pedst) for(int j=0;j<mNChannel[i];j++) 
00130       { 
00131         mPed[i][j] = (float)pedst[0].AdcPedestal[j]/100.;
00132         mPedRms[i][j] = (float)pedst[0].AdcPedestalRMS[j]/100.;
00133         //cout <<"Pedestal "<<i<<"  id = "<<j+1<<" Ped = "<<mPed[i][j]<<"  rms = "<<mPedRms[i][j]<<endl;
00134       }
00135     }
00136     else if(i>=2 && emcDb)
00137     {
00138       St_smdPed *ped = (St_smdPed*)emcDb->Find(table.Data());
00139       smdPed_st *pedst=NULL;
00140       if(ped) pedst=ped->GetTable();
00141       if(pedst) for(int j=0;j<mNChannel[i];j++) // uses only capacitor index 0
00142       { 
00143         mPed[i][j] = (float)pedst[0].AdcPedestal[j][0]/100.;
00144         mPedRms[i][j] = (float)pedst[0].AdcPedestalRMS[j][0]/100.;
00145         //cout <<"Pedestal "<<i<<"  id = "<<j+1<<" Ped = "<<mPed[i][j]<<"  rms = "<<mPedRms[i][j]<<endl;
00146       }           
00147     }
00148     
00149     StDetectorId id = static_cast<StDetectorId>(i+kBarrelEmcTowerId);
00150     StEmcDetector* detector=mEmcCol->detector(id);
00151     if(detector) for(int m=1;m<=120;m++)
00152     {
00153                         StEmcModule* module = detector->module(m);
00154       if(module)
00155       {
00156         StSPtrVecEmcRawHit& rawHit=module->hits();
00157         for(unsigned int k=0;k<rawHit.size();k++) if(rawHit[k])
00158         {
00159           mHasDetector[i] = true;
00160                                         int did;
00161           int mod=rawHit[k]->module();
00162           int e=rawHit[k]->eta(); 
00163           int s=abs(rawHit[k]->sub());
00164           int ADC = rawHit[k]->adc();
00165           int stat = mGeom[i]->getId(mod,e,s,did);
00166           if(stat ==0) 
00167           {
00168             mADC[i][did-1] = ADC;
00169             mADCPedSub[i][did-1] = (short)(ADC - mPed[i][did-1]);
00170             mCap[i][did-1] = (char) rawHit[k]->calibrationType();
00171                                                 if(mCap[i][did-1]>127) mCap[i][did-1]-=128;
00172             if(mPed[i][did-1]==0 || mPedRms[i][did-1]==0) mADCPedSub[i][did-1] = 0;
00173                                                 if(i>=2)
00174                                                 {
00175                                                         if(mCap[i][did-1]==CAP1 || mCap[i][did-1]==CAP1) //uses only capacitor 0 data
00176                                                         {
00177                                                           mADCPedSub[i][did-1] = 0;
00178                                                                 mADC[i][did-1] = 0;
00179                                                         }
00180                                                 }
00181           }
00182         }
00183       }
00184     }
00185   }
00186 }
00187 //_____________________________________________________________________________
00188 void StEmcCalibrationMaker::fillTrackArrays()              
00189 {
00190         StSPtrVecTrackNode& tracks=mStEvent->trackNodes();
00191   mNTracks =  tracks.size();
00192   mL3Tracks = false;
00193         if(mNTracks<=0) 
00194         {
00195                 if(mStEvent->l3Trigger())
00196                 {
00197       tracks =mStEvent->l3Trigger()->trackNodes();
00198                   mNTracks = tracks.size();     
00199                         mL3Tracks = true;
00200                 }
00201         }
00202   if(mNTracks > 0)
00203   {
00204           for(int t=0;t<mNTracks;t++)
00205           { 
00206                   StTrack *track = tracks[t]->track(0);
00207       if(track)
00208       {
00209         if(track->geometry())
00210         {
00211           mTrackPt[t] = track->geometry()->momentum().perp();
00212           mTrackP[t]  = track->geometry()->momentum().mag();
00213           mTrackNPoints[t] = track->fitTraits().numberOfFitPoints();
00214           StThreeVectorD momentum,position;
00215           bool tok = mPosition->trackOnEmc(&position,&momentum,(StTrack*)track,(Double_t)mField,(Double_t)mGeom[0]->Radius());
00216           if(tok)
00217           {
00218             int m,e,s,id=0;
00219             float eta=position.pseudoRapidity();
00220             float phi=position.phi();
00221             int stat = mGeom[0]->getBin(phi,eta,m,e,s);
00222             stat = mGeom[0]->getId(m,e,s,id); 
00223             if(stat==0) 
00224             {
00225               mTrackTower[t] = id;
00226               mNTracksTower[id-1]++;
00227               StThreeVectorD momentum1,position1;
00228               bool tok1 = mPosition->trackOnEmc(&position1,&momentum1,(StTrack*)track,(Double_t)mField,(Double_t)mGeom[0]->Radius()+30.0);
00229               if(tok1)
00230               {
00231                 int id2=0;
00232                 eta=position1.pseudoRapidity();
00233                 phi=position1.phi();
00234                 stat = mGeom[0]->getBin(phi,eta,m,e,s);
00235                 stat = mGeom[0]->getId(m,e,s,id2);
00236                 if(stat==0)
00237                 {
00238                   mTrackTowerExit[t] = id2;
00239                   if(id2!=id) mNTracksTower[id2-1]++;
00240                 }
00241               }
00242             }
00243           }
00244         }
00245       }
00246     }
00247   }
00248   
00249   for(int i=0;i<mNChannel[0];i++)
00250   {
00251     int id = i+1;
00252     mIsIsolatedTower[id-1] = true;
00253     for(int j=-1;j<=1;j++) for(int k=-1;k<=1;k++) 
00254     {
00255       int id2 = mPosition->getNextTowerId(id,j,k);
00256       if(id2!=id) if(id2>=1 && id2<=4800) if(mNTracksTower[id2-1]>0) { mIsIsolatedTower[id-1] = false; break; }
00257     }
00258   }
00259 }
00260 void StEmcCalibrationMaker::makeStatus(bool t, bool p, bool se, bool sp)
00261 {
00262   // for the towers
00264         int towers[4800];
00265         for(int i =0;i<4800;i++) towers[i] = 1; // include all towers
00266         
00267         int crates[] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
00268                         1,1,1,1,0,1,1,1,1,1,0,1,1,1,1}; // exclude some crates
00269         
00271         // for the PSD
00272         int psd[4800];
00273         for(int i =0;i<4800;i++) psd[i] = 1; // include all PSD towers
00274         for(int i =2400;i<4800;i++) psd[i] = 0; // include all PSD towers
00275         
00277         // for the SMD-eta
00278         int smde[18000];
00279         for(int i =0;i<18000;i++) smde[i] = 1; // include all SMD-eta towers
00280         
00281         int smdeMod[120];
00282         for(int i =0;i<120;i++)  smdeMod[i] = 1; // include all SMD-eta modules
00283         for(int i =60;i<120;i++) smdeMod[i] = 0; // include all SMD-eta modules
00284         
00286         // for the SMD-phi
00287         int smdp[18000];
00288         for(int i =0;i<18000;i++) smdp[i] = 1; // include all SMD-phi towers
00289         
00290         int smdpMod[120];
00291         for(int i =0;i<120;i++)  smdpMod[i] = 1; // include all SMD-phi modules
00292         for(int i =60;i<120;i++) smdpMod[i] = 0; // include all SMD-eta modules
00293         
00295         StEmcDecoder *dec = new StEmcDecoder(20300101,0);
00296         
00297         const char *timestamp = "2004-01-01 00:00:00";
00298         
00299   StDbManager* mgr=StDbManager::Instance();
00300         StDbConfigNode* node=mgr->initConfig(dbCalibrations,dbEmc);
00301         
00302         // for the towers
00303         
00304         if(t)
00305         {
00306           St_emcStatus *mBemc = new St_emcStatus("emcStatus",1);
00307                 emcStatus_st *bemc=mBemc->GetTable();
00308         
00309                 for(int i=0;i<4800;i++) 
00310                 {
00311                 bemc[0].Status[i] = towers[i];
00312                         int daq;
00313                         dec->GetDaqIdFromTowerId(i+1,daq);
00314                         int crate, pos;
00315                         dec->GetTowerCrateFromDaqId(daq,crate,pos);
00316                         bemc[0].Status[i] = crates[crate-1];
00317                 }
00318         
00319                 StDbTable* tab1=node->addDbTable("bemcStatus");
00320                 tab1->SetTable((char*)bemc,1);
00321                 mgr->setStoreTime(timestamp);
00322                 mgr->storeDbTable(tab1);
00323         }
00324         
00325         // for the PSD
00326         
00327   if(p)
00328         {
00329                 St_emcStatus *mPSD = new St_emcStatus("emcStatus",1);
00330                 emcStatus_st *psdt = mPSD->GetTable();
00331         
00332                 for(int i=0;i<4800;i++) 
00333                 {
00334                 psdt[0].Status[i] = psd[i];
00335                 }
00336         
00337                 StDbTable* tab2=node->addDbTable("bprsStatus");
00338                 tab2->SetTable((char*)psdt,1);
00339                 mgr->setStoreTime(timestamp);
00340                 mgr->storeDbTable(tab2);
00341         }
00342         
00343         // for the SMD-eta
00344         
00345         if(se)
00346         {
00347     St_smdStatus *mSMDE = new St_smdStatus("smdStatus",1);
00348                 smdStatus_st *smdet = mSMDE->GetTable();
00349                 StEmcGeom *geo3 = StEmcGeom::instance("bsmde");
00350         
00351                 for(int i=0;i<18000;i++) 
00352                 {
00353                 smdet[0].Status[i] = smde[i];
00354                         int m,e,s;
00355                         geo3->getBin(i+1,m,e,s);
00356                 smdet[0].Status[i] = smdeMod[m-1];
00357                 }
00358         
00359                 StDbTable* tab3=node->addDbTable("bsmdeStatus");
00360                 tab3->SetTable((char*)smdet,1);
00361                 mgr->setStoreTime(timestamp);
00362                 mgr->storeDbTable(tab3);
00363         }
00364         
00365         // for the SMD-phi
00366         if(sp)
00367         {
00368         St_smdStatus *mSMDP = new St_smdStatus("smdStatus",1);
00369                 smdStatus_st *smdpt = mSMDP->GetTable();
00370                 StEmcGeom *geo4 = StEmcGeom::instance("bsmdp");
00371         
00372                 for(int i=0;i<18000;i++) 
00373                 {
00374                 smdpt[0].Status[i] = smdp[i];
00375                         int m,e,s;
00376                         geo4->getBin(i+1,m,e,s);
00377                 smdpt[0].Status[i] = smdpMod[m-1];
00378                 }
00379         
00380                 StDbTable* tab4=node->addDbTable("bsmdpStatus");
00381                 tab4->SetTable((char*)smdpt,1);
00382                 mgr->setStoreTime(timestamp);
00383                 mgr->storeDbTable(tab4);
00384         }
00385         
00386         delete dec;
00387         dec = 0;
00388   
00389 }