StRoot: StEmcPool/StEmcOfflineCalibrationMaker/StEmcOfflineCalibrationMaker.cxx Source File

00001 //StEmcOfflineCalibrationMaker.cxx
00002 
00003 #include <math.h>
00004 #include <vector>
00005 //#include <algorithm>
00006 #include <stdio.h>
00007 
00008 #include "TFile.h"
00009 #include "TTree.h"
00010 #include "TH2.h"
00011 #include "TH3.h"
00012 #include "TChain.h"
00013 #include "TCanvas.h"
00014 
00015 //StMuDstMaker
00016 #include "StMuDSTMaker/COMMON/StMuDst.h"
00017 #include "StMuDSTMaker/COMMON/StMuEvent.h"
00018 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
00019 #include "StMuDSTMaker/COMMON/StMuEmcCollection.h"
00020 #include "StMuDSTMaker/COMMON/StMuEmcUtil.h"
00021 #include "StMuDSTMaker/COMMON/StMuTrack.h"
00022 #include "StMuDSTMaker/COMMON/StMuPrimaryVertex.h"
00023 
00024 #include "StEvent/StTriggerId.h"
00025 #include "StEvent/StEvent.h"
00026 #include "StEventTypes.h"
00027 
00028 #include "StDetectorDbMaker/StDetectorDbTriggerID.h"
00029 #include "StEvent/StL0Trigger.h"
00030 
00031 #include "StDaqLib/TRG/trgStructures2005.h"
00032 #include "StarClassLibrary/StPhysicalHelixD.hh"
00033 //StEmc
00034 #include "StEmcUtil/database/StBemcTables.h"
00035 #include "StEmcUtil/geometry/StEmcGeom.h"
00036 #include "StEmcUtil/projection/StEmcPosition.h"
00037 #include "StEmcADCtoEMaker/StBemcData.h"
00038 #include "StEmcADCtoEMaker/StEmcADCtoEMaker.h"
00039 #include "StEmcRawMaker/StBemcRaw.h"
00040 #include "StEmcRawMaker/defines.h"
00041 #include "StEmcRawMaker/StBemcTables.h"
00042 #include "StEmcTriggerMaker/StEmcTriggerMaker.h"
00043 #include "StTriggerUtilities/StTriggerSimuMaker.h"
00044 #include "StTriggerUtilities/Bemc/StBemcTriggerSimu.h"
00045 
00046 //logger
00047 #include "StMessMgr.h"
00048 
00049 //my files
00050 #include "StEmcOfflineCalibrationEvent.h"
00051 #include "StEmcOfflineCalibrationMaker.h"
00052 
00053 //use this to get kBarrelEmcTowerId defined
00054 #include "StEnumerations.h"
00055 
00056 using namespace std;
00057 
00058 ClassImp(StEmcOfflineCalibrationMaker)
00059 
00060 StEmcOfflineCalibrationMaker::StEmcOfflineCalibrationMaker(const char* name, const char* file)
00061 {
00062   filename = file;
00063         
00064   muDstMaker = NULL;
00065   mEmcPosition = new StEmcPosition();
00066         
00067   towerSlopes[0] = NULL;
00068   towerSlopes[1] = NULL;
00069   smdSlopes[0]=NULL;
00070   smdSlopes[1]=NULL;
00071   preshowerSlopes = NULL;
00072   mapcheck = NULL;
00073   mbTriggers.clear();
00074   htTriggers.clear();
00075 }
00076 
00077 StEmcOfflineCalibrationMaker::~StEmcOfflineCalibrationMaker() { }
00078 
00079 Int_t StEmcOfflineCalibrationMaker::Init()
00080 {
00081   TString mfname;
00082   mfname += filename;
00083   mfname.ReplaceAll(".calib.",".map.");
00084   mapFile = new TFile(mfname.Data(),"RECREATE");
00085   myFile = new TFile(filename,"RECREATE");
00086   calibTree = new TTree("calibTree","BTOW calibration tree");
00087   myEvent = new StEmcOfflineCalibrationEvent();
00088   myTrack = new StEmcOfflineCalibrationTrack();
00089   calibTree->Branch("event_branch","StEmcOfflineCalibrationEvent",&myEvent);
00090   calibTree->SetAutoSave(1000000000);
00091         
00092   muDstMaker = dynamic_cast<StMuDstMaker*>(GetMaker("MuDst")); assert(muDstMaker);
00093   mADCtoEMaker = dynamic_cast<StEmcADCtoEMaker*>(GetMaker("Eread")); assert(mADCtoEMaker);
00094   //emcTrigMaker = dynamic_cast<StEmcTriggerMaker*>(GetMaker("bemctrigger")); assert(emcTrigMaker);
00095   emcTrigMaker = dynamic_cast<StTriggerSimuMaker*>(GetMaker("StarTrigSimu")); assert(emcTrigMaker);
00096 
00097 
00098   mTables = mADCtoEMaker->getBemcData()->getTables();
00099   mEmcGeom = StEmcGeom::instance("bemc");
00100   mSmdEGeom=StEmcGeom::instance("bsmde");
00101   mSmdPGeom=StEmcGeom::instance("bsmdp");
00102 
00103   LOG_INFO << "StEmcOfflineCalibrationMaker::Init() == kStOk" << endm;
00104   return StMaker::Init();
00105 }
00106 
00107 Int_t StEmcOfflineCalibrationMaker::InitRun(int run)
00108 {
00109   LOG_DEBUG<<"StEmcOfflineCalibrationMaker::InitRun"<<endm;
00110         //switch histograms
00111   if(towerSlopes[0]){
00112     myFile->cd();
00113     towerSlopes[0]->Write();
00114     towerSlopes[1]->Write();
00115     preshowerSlopes->Write();
00116     smdSlopes[0]->Write();
00117     smdSlopes[1]->Write();
00118     mapcheck->Write();
00119   }
00120 
00121   mHT0threshold = emcTrigMaker->bemc->barrelHighTowerTh(0);
00122   mHT1threshold = emcTrigMaker->bemc->barrelHighTowerTh(1);
00123   mHT2threshold = emcTrigMaker->bemc->barrelHighTowerTh(2);
00124   mHT3threshold = emcTrigMaker->bemc->barrelHighTowerTh(3);
00125         
00126         //look for histograms from this run in the current file and switch to them if found, otherwise create them
00127   char name[200];
00128   sprintf(name,"towerSlopes_R%i",run);
00129   towerSlopes[0] = NULL;
00130   myFile->GetObject(name,towerSlopes[0]);
00131   if(towerSlopes[0]){
00132     sprintf(name,"towerSlopes_HT_R%i",run);
00133     myFile->GetObject(name,towerSlopes[1]);
00134                 
00135     sprintf(name,"preshowerSlopes_R%i",run);
00136     myFile->GetObject(name,preshowerSlopes);
00137 
00138     sprintf(name,"smdSlopesE_R%i",run);
00139     myFile->GetObject(name,smdSlopes[0]);
00140 
00141     sprintf(name,"smdSlopesP_R%i",run);
00142     myFile->GetObject(name,smdSlopes[1]);
00143   }else{
00144     towerSlopes[0] = new TH2F(name,"ADC vs. towerID",4800,0.5,4800.5,250,-49.5,200.5);
00145                
00146     sprintf(name,"towerSlopes_HT_R%i",run);
00147     towerSlopes[1] = new TH2F(name,"ADC vs. towerID",4800,0.5,4800.5,250,-49.5,200.5);
00148                 
00149     sprintf(name,"preshowerSlopes_R%i",run);            
00150     preshowerSlopes = new TH2F(name,"ADC vs. cap vs. towerID",4800,0.5,4800.5,500,-49.5,450.5);
00151 
00152     sprintf(name,"smdSlopesE_R%i",run);
00153     smdSlopes[0]=new TH2F(name,"ADC vs. stripID",18000,0.5,18000.5,1223,-199.5,1023.5);
00154 
00155     sprintf(name,"smdSlopesP_R%i",run);
00156     smdSlopes[1]=new TH2F(name,"ADC vs. stripID",18000,0.5,18000.5,1223,-199.5,1023.5);
00157   }
00158   if(!mapcheck){
00159     mapcheck = new TH2F("mapcheck","Track Projection vs EMC hit",4800,0.5,4800.5,4800,0.5,4800.5);
00160   }     
00161         //pedestals, rms, status
00162 
00163 
00164   float pedestal, rms;
00165   int status;
00166   int pedstatus;
00167   for (int id=1; id<=4800; id++){
00168 
00169     mTables->getPedestal(BTOW, id, 0, pedestal, rms);
00170     mTables->getStatus(BTOW, id, pedstatus, "pedestal");
00171     mTables->getStatus(BTOW, id, status);
00172     mPedestal[BTOW-1][id-1] = pedestal;
00173     mPedRMS[BTOW-1][id-1] = rms;
00174     mStatus[BTOW-1][id-1] = status * pedstatus;
00175                 
00176     mTables->getPedestal(BPRS, id, 0, pedestal, rms);
00177     mTables->getStatus(BPRS, id, status);
00178     mPedestal[BPRS-1][id-1] = pedestal;
00179     mPedRMS[BPRS-1][id-1]       = rms;
00180     mStatus[BPRS-1][id-1]       = status;               
00181   }
00182         
00183   for(int sid=1;sid<=18000;sid++){
00184                 
00185     mTables->getPedestal(BSMDE, sid, 0, pedestal, rms);
00186     mTables->getStatus(BSMDE, sid, status);
00187     mPedestalSmd[BSMDE-3][sid-1] = pedestal;
00188     mPedRMSSmd[BSMDE-3][sid-1]  = rms;
00189     mStatusSmd[BSMDE-3][sid-1]  = status;
00190                 
00191     mTables->getPedestal(BSMDP, sid, 0, pedestal, rms);
00192     mTables->getStatus(BSMDP, sid, status);
00193     mPedestalSmd[BSMDP-3][sid-1] = pedestal;
00194     mPedRMSSmd[BSMDP-3][sid-1]  = rms;
00195     mStatusSmd[BSMDP-3][sid-1]  = status;               
00196   }
00197         
00198 
00199   LOG_DEBUG << "got peds and status in InitRun()" << endm;
00200   LOG_DEBUG << "finished init run for " << run << endm;
00201   return StMaker::InitRun(run);
00202 }
00203 
00204 Int_t StEmcOfflineCalibrationMaker::Make()
00205 {               
00206   LOG_DEBUG<<"StEmcOfflineCalibrationMaker::Make()"<<endm;
00207         //get pointers to useful objects
00208   TChain* chain = muDstMaker->chain(); assert(chain);
00209   StMuDst* muDst = muDstMaker->muDst(); assert(muDst);
00210   StMuEvent* event = muDst->event(); assert(event);
00211   StRunInfo* runInfo = &(event->runInfo()); assert(runInfo);
00212 //      TrgDataType2005* trgData        = (TrgDataType2005*)event->triggerData()->getTriggerStructure(); assert(trgData);
00213         
00214   LOG_DEBUG << "general pointer assertions OK" << endm;
00215         
00216   const StTriggerId& trigs = event->triggerIdCollection().nominal();
00217         //toss out fpd1-tpcdead-fast triggers immediately
00218   for(unsigned int i = 0; i < fastTriggers.size(); i++){
00219     if(trigs.isTrigger(fastTriggers[i]))return kStOK;
00220   } 
00221   HT0towersAboveThreshold.clear();
00222   HT1towersAboveThreshold.clear();
00223   HT2towersAboveThreshold.clear();
00224   HT3towersAboveThreshold.clear();
00225 
00226   myEvent->triggerResult.clear();
00227   myEvent->towersAboveThreshold.clear();
00228         
00229   myEvent->triggerIds = trigs.triggerIds();
00230   myEvent->l2Result = event->L2Result();        
00231         //basic event info
00232   bool fillQA = (runInfo->beamFillNumber(blue)==runInfo->beamFillNumber(yellow));
00233   myEvent->fill = (fillQA) ? (unsigned short)runInfo->beamFillNumber(blue):0;
00234   myEvent->run = event->runNumber();
00235   myEvent->event = event->eventNumber();
00236   myEvent->date = GetDate();
00237   myEvent->time = GetTime();
00238         
00239         //filename manipulations
00240   TString inputfile(chain->GetFile()->GetName());
00241   int index1 = inputfile.Index(".MuDst");
00242   TString string1(inputfile(index1-19,7));
00243   TString string2(inputfile(index1-7,7));
00244   myEvent->fileid1 = string1.Atoi();
00245   myEvent->fileid2 = string2.Atoi();    
00246         
00247         //store all vertices from this event    
00248   myEvent->nVertices    = muDst->numberOfPrimaryVertices();
00249   for(unsigned int i=0; i<myEvent->nVertices; i++){
00250     if(i>9){
00251       LOG_WARN << "found more than 10 vertices for R"<<myEvent->run<< ", event "<<myEvent->event<<endm;
00252       continue;
00253     }
00254                 
00255     assert(muDst->primaryVertex(i));
00256     StThreeVectorF stvertex = muDst->primaryVertex(i)->position();
00257     myEvent->vx[i] = stvertex.x();
00258     myEvent->vy[i] = stvertex.y();
00259     myEvent->vz[i] = stvertex.z();
00260     myEvent->ranking[i] = muDst->primaryVertex(i)->ranking();
00261   }
00262                 
00263         //fill ADC values from StEmcCollection obtained from MuDst
00264   //mEmcCollection = muDst->emcCollection();
00265   mEmcCollection = mADCtoEMaker->getEmcCollection();
00266   for(int p=0;p<2;p++)
00267     for(int q=0;q<18000;q++) mADCSmd[p][q]=0;
00268   getADCs(BTOW);
00269   getADCs(BPRS);
00270   getADCs(BSMDP);
00271   getADCs(BSMDE);
00272   for(int id=1; id<=4800; id++){
00273     if((mADC[BTOW-1][id-1] != 0)){
00274       //if(myEvent->mbTrigger)
00275       towerSlopes[0]->Fill(id,mADC[BTOW-1][id-1]-mPedestal[BTOW-1][id-1]);
00276       //if(myEvent->htTrigger) towerSlopes[1]->Fill(id,mADC[BTOW-1][id-1]-mPedestal[BTOW-1][id-1]); 
00277     }
00278     if(mADC[BPRS-1][id-1] != 0 && mCapacitor[id-1]!= CAP1 && mCapacitor[id-1]!=CAP2){           
00279       //                LOG_DEBUG << "filling preshower slopes histogram" << endm;
00280       preshowerSlopes->Fill(id, mADC[BPRS-1][id-1]-mPedestal[BPRS-1][id-1]);
00281     }
00282   }
00283         
00284   for(int id=1; id<=18000; id++){
00285     if((mADCSmd[BSMDE-3][id-1] != 0) && mCapacitorSmd[BSMDE-3][id-1] != CAP1 && mCapacitorSmd[BSMDE-3][id-1] != CAP2){
00286                         //if(myEvent->mbTrigger)
00287       smdSlopes[0]->Fill(id,mADCSmd[BSMDE-3][id-1]);
00288                         //if(myEvent->htTrigger) towerSlopes[1]->Fill(id,mADC[BTOW-1][id-1]-mPedestal[BTOW-1][id-1]); 
00289     }
00290     if(mADCSmd[BSMDP-3][id-1] != 0 && mCapacitorSmd[BSMDP-3][id-1] != CAP1 && mCapacitorSmd[BSMDP-3][id-1] != CAP2){
00291         //              LOG_DEBUG << "filling preshower slopes histogram" << endm;
00292       smdSlopes[1]->Fill(id, mADCSmd[BSMDP-3][id-1]);
00293     }
00294   }
00295         
00296         
00297   LOG_DEBUG << "got ADCs" << endm;
00298 
00299 
00300         
00301         //trigger maker
00302   for(unsigned int i = 0; i < htTriggers.size(); i++){
00303     if(trigs.isTrigger(htTriggers[i])){
00304       myEvent->triggerResult[htTriggers[i]]=emcTrigMaker->isTrigger(htTriggers[i]);
00305       if(emcTrigMaker->bemc->getTowersAboveThreshold((int)htTriggers[i]).size() > 0){
00306       myEvent->towersAboveThreshold[htTriggers[i]] = emcTrigMaker->bemc->getTowersAboveThreshold((int)htTriggers[i]);
00307       }else{
00308         if(htTriggers[i] == 240540)myEvent->towersAboveThreshold[htTriggers[i]] = HT2towersAboveThreshold;
00309         if(htTriggers[i] == 240560)myEvent->towersAboveThreshold[htTriggers[i]] = HT1towersAboveThreshold;
00310         if(htTriggers[i] == 240570 || htTriggers[i] == 240550)myEvent->towersAboveThreshold[htTriggers[i]] = HT0towersAboveThreshold;
00311       }
00312     }
00313   }
00314   for(unsigned int i = 0; i < httpTriggers.size(); i++){
00315     myEvent->triggerResult[httpTriggers[i]]=emcTrigMaker->isTrigger(httpTriggers[i]);
00316     myEvent->towersAboveThreshold[httpTriggers[i]] = emcTrigMaker->bemc->getTriggerPatchesAboveThreshold(httpTriggers[i]);
00317   }
00318 
00319 
00320   /*
00321   for(unsigned int i = 0; i < mbTriggers.size(); i++){
00322     myEvent->triggerResult[mbTriggers[i]]=emcTrigMaker->isTrigger(mbTriggers[i]);
00323   }
00324   */
00325   myEvent->htTrigMaker[0] = emcTrigMaker->isTrigger(137213);
00326   //std::map<int,int>::const_iterator p = (emcTrigMaker->barrelTowersAboveThreshold(137213)).begin();
00327   //myEvent->htTrigMaker[1] = p->first;
00328   //myEvent->htTrigMaker[2] = p->second;
00329 
00330   const StEventSummary& evtSummary = muDstMaker->muDst()->event()->eventSummary();
00331   Double_t mField = evtSummary.magneticField()/10;
00332   //now for the tracks
00333 
00334   vector<int> track_ids; 
00335 
00336  for(unsigned int vertex_index=0; vertex_index<myEvent->nVertices; vertex_index++){
00337     muDst->setVertexIndex(vertex_index);
00338     TObjArray* primaryTracks = muDst->primaryTracks();
00339     const StMuTrack* track;
00340     const StMuTrack* primarytrack;
00341 
00342     int nentries = muDst->numberOfPrimaryTracks();
00343     //cout<<nentries<<" tracks in vertex "<<vertex_index<<endl;
00344     assert(nentries==primaryTracks->GetEntries());
00345     pair<unsigned int, pair<float,float> > smd_eta_center;
00346     pair<unsigned int, pair<float,float> > smd_phi_center;
00347 
00348                 //cout<<nentries<<endl;
00349     for(int i=0; i<nentries; i++){
00350       track = NULL;
00351       primarytrack = muDst->primaryTracks(i);
00352       myTrack->flag = primarytrack->flag();
00353       myTrack->bad = primarytrack->bad();
00354 
00355       //if(primarytrack->flag()<=0)continue;
00356       //if(primarytrack->bad())continue;
00357       //if((float)(primarytrack->nHitsFit()/primarytrack->nHitsPoss()) < 0.51)continue;
00358       //if(primarytrack->eta() > 2)continue;
00359       track = primarytrack->globalTrack();
00360       int primused = 0;
00361       if(!track)primused=1;
00362       if(!track)track = primarytrack;
00363       myTrack->charge = track->charge();
00364       myTrack->eta = track->eta();
00365       double p;
00366       StPhysicalHelixD outerhelix;
00367       pair<unsigned int, pair<float,float> > center_tower;
00368       outerhelix = track->outerHelix();
00369       p = outerhelix.momentum(mField*tesla).mag();
00370       //float testeta = outerhelix.momentum(mField*tesla).pseudoRapidity();
00371       myTrack->track = outerhelix.momentum(mField*tesla);
00372       //p = track->p().mag();
00373       center_tower = getTrackTower(track, false,1);
00374       smd_eta_center=getTrackTower(track,false,3);
00375       smd_phi_center=getTrackTower(track,false,4);
00376 
00377 
00378                         //project track to BEMC
00379 
00380       int id = center_tower.first;
00381       int exitid = (getTrackTower(track, true)).first;
00382 
00383                         
00384                         //cout<<p<<endl;
00385 
00386       //cout<<i<<" "<<id<<" "<<primused<<" "<<p<<" "<<testeta<<endl;
00387                         
00388       if(id > 0 && p > 1.0){
00389         if (id == exitid)track_ids.push_back(id);
00390         int etaid=smd_eta_center.first;
00391         int phiid=smd_phi_center.first;
00392         LOG_DEBUG<<"using track projection to the SMD, we get strips eta: "<<etaid<<" and phi: "<<phiid<<endm;
00393         if(etaid>0){
00394           int smdeids[11];
00395           smdeids[0] = etaid;
00396           smdeids[1] = mEmcPosition->getNextId(3,etaid,1,0);
00397           smdeids[2] = mEmcPosition->getNextId(3,etaid,2,0);
00398           smdeids[3] = mEmcPosition->getNextId(3,etaid,3,0);
00399           smdeids[4] = mEmcPosition->getNextId(3,etaid,4,0);
00400           smdeids[5] = mEmcPosition->getNextId(3,etaid,5,0);
00401           smdeids[6] = mEmcPosition->getNextId(3,etaid,-1,0);
00402           smdeids[7] = mEmcPosition->getNextId(3,etaid,-2,0);
00403           smdeids[8] = mEmcPosition->getNextId(3,etaid,-3,0);
00404           smdeids[9] = mEmcPosition->getNextId(3,etaid,-4,0);
00405           smdeids[10] = mEmcPosition->getNextId(3,etaid,-5,0);
00406 
00407           for(int i = 0; i < 11; i++){
00408             myTrack->smde_id[i] = smdeids[i];
00409             myTrack->smde_adc[i] = mADCSmd[0][smdeids[i]-1];
00410             myTrack->smde_pedestal[i] = mPedestalSmd[0][smdeids[i]-1];
00411             myTrack->smde_pedestal_rms[i] = mPedRMSSmd[0][smdeids[i]-1];
00412             myTrack->smde_status[i] = mStatusSmd[0][smdeids[i]-1];
00413           }
00414         }else{
00415           for(int i=0;i<11;i++){
00416             myTrack->smde_id[i]=0;
00417             myTrack->smde_adc[i]=0;
00418             myTrack->smde_pedestal[i]=0;
00419             myTrack->smde_pedestal_rms[i]=0;
00420             myTrack->smde_status[i]=0;
00421           }
00422         }
00423         if(phiid>0){
00424           int smdpids[11];
00425           smdpids[0] = phiid;
00426           smdpids[1] = mEmcPosition->getNextId(4,phiid,0,1);
00427           smdpids[2] = mEmcPosition->getNextId(4,phiid,0,2);
00428           smdpids[3] = mEmcPosition->getNextId(4,phiid,0,3);
00429           smdpids[4] = mEmcPosition->getNextId(4,phiid,0,4);
00430           smdpids[5] = mEmcPosition->getNextId(4,phiid,0,5);
00431           smdpids[6] = mEmcPosition->getNextId(4,phiid,0,-1);
00432           smdpids[7] = mEmcPosition->getNextId(4,phiid,0,-2);
00433           smdpids[8] = mEmcPosition->getNextId(4,phiid,0,-3);
00434           smdpids[9] = mEmcPosition->getNextId(4,phiid,0,-4);
00435           smdpids[10] = mEmcPosition->getNextId(4,phiid,0,-5);
00436 
00437           for(int i = 0; i < 11; i++){
00438             myTrack->smdp_id[i] = smdpids[i];
00439             myTrack->smdp_adc[i] = mADCSmd[1][smdpids[i]-1];
00440             myTrack->smdp_pedestal[i] = mPedestalSmd[1][smdpids[i]-1];
00441             myTrack->smdp_pedestal_rms[i] = mPedRMSSmd[1][smdpids[i]-1];
00442             myTrack->smdp_status[i] = mStatusSmd[1][smdpids[i]-1];
00443           }
00444         }else{
00445           for(int i=0;i<11;i++){
00446             myTrack->smdp_id[i]=0;
00447             myTrack->smdp_adc[i]=0;
00448             myTrack->smdp_pedestal[i]=0;
00449             myTrack->smdp_pedestal_rms[i]=0;
00450             myTrack->smdp_status[i]=0;
00451           }
00452         }
00453                           
00454 
00455                                 //find neighboring towers
00456         int softid[9];
00457         softid[0] = id;
00458         softid[1] = mEmcPosition->getNextTowerId(id,-1,-1);
00459         softid[2] = mEmcPosition->getNextTowerId(id,0,-1);
00460         softid[3] = mEmcPosition->getNextTowerId(id,1,-1);
00461         softid[4] = mEmcPosition->getNextTowerId(id,-1,0);
00462         softid[5] = mEmcPosition->getNextTowerId(id,1,0);
00463         softid[6] = mEmcPosition->getNextTowerId(id,-1,1);
00464         softid[7] = mEmcPosition->getNextTowerId(id,0,1);
00465         softid[8] = mEmcPosition->getNextTowerId(id,1,1);                                                               
00466                                 //save EMC info for 3x3 tower block around track
00467         for(int tower=0; tower<9; tower++){
00468           myTrack->tower_id[tower] = softid[tower];
00469           myTrack->tower_adc[tower] = mADC[BTOW-1][softid[tower]-1];
00470           myTrack->tower_pedestal[tower] = mPedestal[BTOW-1][softid[tower]-1];
00471           myTrack->tower_pedestal_rms[tower] = mPedRMS[BTOW-1][softid[tower]-1];
00472           myTrack->tower_status[tower] = mStatus[BTOW-1][softid[tower]-1];
00473         }
00474                                                         
00475         for(int tower=0; tower<9; tower++){
00476           myTrack->preshower_adc[tower] = mADC[BPRS-1][softid[tower]-1];
00477           myTrack->preshower_pedestal[tower] = mPedestal[BPRS-1][softid[tower]-1];
00478           myTrack->preshower_pedestal_rms[tower] = mPedRMS[BPRS-1][softid[tower]-1];
00479           myTrack->preshower_status[tower] = mStatus[BPRS-1][softid[tower]-1];
00480           myTrack->preshower_cap[tower] = mCapacitor[softid[tower]-1];
00481         }
00482                                 
00483         myTrack->p = p;
00484         myTrack->deta = center_tower.second.first;
00485         myTrack->dphi = center_tower.second.second;
00486         myTrack->tower_id_exit = (getTrackTower(track, true)).first;
00487         myTrack->highest_neighbor = highestNeighbor(myTrack->tower_id[0]);
00488                                 
00489         myTrack->nSigmaElectron = track->nSigmaElectron();
00490         myTrack->nHits = track->nHits();
00491         myTrack->nFitPoints = track->nHitsFit();
00492         myTrack->nDedxPoints = track->nHitsDedx();
00493         myTrack->nHitsPossible = track->nHitsPoss();
00494         myTrack->dEdx = track->dEdx();
00495                                 
00496         myEvent->addTrack(myTrack);
00497         //                      cout<<"I think I added a track"<<endl;
00498       }
00499     }
00500   }
00501 
00502  for(int i = 0; i < 4800; i++){
00503    if((mADC[BTOW-1][i]-mPedestal[BTOW-1][i]) < 5 || (mADC[BTOW-1][i]-mPedestal[BTOW-1][i]) > 100)continue;
00504    for(unsigned int j = 0; j < track_ids.size(); j++){
00505      mapcheck->Fill(track_ids[j],i+1);
00506    }
00507  }
00508 
00509   calibTree->Fill();
00510   myEvent->Clear();
00511   LOG_DEBUG << "StEmcOfflineCalibrationMaker::Make() == kStOk" << endm;
00512   return kStOK;
00513 }
00514 
00515 void StEmcOfflineCalibrationMaker::Clear(Option_t* option)
00516 {       
00517         myEvent->Clear();
00518         for(int i=0; i<2; i++){
00519                 for(int j=0; j<4800; j++){
00520                         mADC[i][j] = 0;
00521                 }
00522         }
00523         for(int i=0; i<2; i++){
00524           for(int j=0; j<18000; j++){
00525             mADCSmd[i][j] = 0;
00526           }
00527         }
00528 }
00529 
00530 Int_t StEmcOfflineCalibrationMaker::Finish()
00531 {
00532         LOG_INFO << "cd to output file" << endm;
00533         myFile->cd();
00534         LOG_INFO << "write calibration tree" << calibTree->Write() << endm;
00535         LOG_INFO << "write tower histogram" << towerSlopes[0]->Write() << endm;
00536         towerSlopes[1]->Write();
00537         LOG_INFO << "write preshower histogram" << preshowerSlopes->Write() << endm;
00538         LOG_INFO << "write smd histogram" <<smdSlopes[0]->Write() << endm;
00539         smdSlopes[1]->Write();
00540         LOG_INFO << "close the output file" << endm;
00541         myFile->Close();
00542         mapFile->cd();
00543         mapcheck->Write();
00544         mapFile->Close();
00545         LOG_INFO << "StEmcOfflineCalibrationMaker::Finish() == kStOk"<<endm;
00546         return kStOk;
00547 }
00548 
00549 void StEmcOfflineCalibrationMaker::getADCs(int det) //1==BTOW, 2==BPRS, 3=BSMDE, 4=BSMDP
00550 {
00551   det--;
00552   StDetectorId detectorid = static_cast<StDetectorId>(kBarrelEmcTowerId + det);
00553   StEmcDetector* detector=mEmcCollection->detector(detectorid);
00554   if(detector){
00555     for(int m=1;m<=120;m++){
00556       StEmcModule* module = detector->module(m);
00557       if(module)
00558         {
00559           StSPtrVecEmcRawHit& rawHit=module->hits();
00560           for(unsigned int k=0;k<rawHit.size();k++){
00561             if(rawHit[k]){
00562               int module = rawHit[k]->module();
00563               int eta = rawHit[k]->eta(); 
00564               int submodule = TMath::Abs(rawHit[k]->sub());
00565               int ADC = rawHit[k]->adc();
00566               int hitid;
00567               int stat=-9;            
00568               if(det<2) stat = mEmcGeom->getId(module,eta,submodule,hitid);
00569               else if(det==2) stat=mSmdEGeom->getId(module,eta,submodule,hitid);
00570               else if(det==3) stat=mSmdPGeom->getId(module,eta,submodule,hitid);
00571               if(det == 0 && stat == 0){
00572                 int adc6 = emcTrigMaker->bemc->getHT6bitAdc(hitid);
00573                 pair<int,int> httowerthreshold (hitid,adc6);
00574                 if(adc6 > mHT0threshold)HT0towersAboveThreshold.push_back(httowerthreshold);
00575                 if(adc6 > mHT1threshold)HT1towersAboveThreshold.push_back(httowerthreshold);
00576                 if(adc6 > mHT2threshold)HT2towersAboveThreshold.push_back(httowerthreshold);
00577                 if(adc6 > mHT3threshold)HT3towersAboveThreshold.push_back(httowerthreshold);
00578               }
00579               if(stat==0){
00580                 if(det<2) mADC[det][hitid-1] = ADC;
00581                 else{
00582                   mADCSmd[det-2][hitid-1]=ADC;
00583                   //cout<<"added SMD hit"<<endl;
00584                 }
00585               }
00586               if(det==1){
00587                 unsigned char CAP =  rawHit[k]->calibrationType();
00588                 if(CAP > 127) CAP -= 128;
00589                 mCapacitor[hitid-1] = CAP;
00590               }
00591               if(det>1){
00592                 unsigned char CAP=rawHit[k]->calibrationType();
00593                 if(CAP>127) CAP-=128;
00594                 mCapacitorSmd[det-2][hitid-1]=CAP;
00595               }
00596             }
00597           }
00598         }
00599       else { LOG_WARN<<"couldn't find StEmcModule "<<m<<" for detector "<<det<<endm; }
00600     }
00601   }
00602   else { LOG_WARN<<"couldn't find StEmcDetector for detector "<<det<<endm; }
00603 }
00604 
00605 /*
00606 //returns (tower, (deta,dphi))
00607 pair<unsigned short, pair<float,float> > StEmcOfflineCalibrationMaker::getTrackTower(StMuTrack* track, bool useExitRadius){
00608         pair<unsigned short, pair<float,float> > tower;
00609         tower.first = 0;
00610         tower.second.first = 1000.;
00611         tower.second.second = 1000.;
00612 
00613         StThreeVectorD momentum,position;
00614         Double_t radius = mEmcGeom->Radius();
00615                 
00616         const StEventSummary& evtSummary = muDstMaker->muDst()->event()->eventSummary();
00617         Double_t mField = evtSummary.magneticField()/10;
00618                 
00619         //add 30 cm to radius to find out if track left same tower
00620         if(useExitRadius) radius += 30.0;
00621         
00622         bool goodProjection = mEmcPosition->trackOnEmc(&position,&momentum,track,mField,radius);
00623         if(goodProjection){
00624                 int m,e,s,id=0;
00625                 float eta=position.pseudoRapidity();
00626                 float phi=position.phi();
00627                 mEmcGeom->getBin(phi,eta,m,e,s);
00628                 
00629                 if(mEmcGeom->getId(m,e,s,id)==0){
00630                         tower.first = id;
00631                         tower.second = getTrackDetaDphi(eta, phi, id);
00632                 }
00633         }
00634         
00635         return tower;
00636 }
00637 
00638 pair<float,float> StEmcOfflineCalibrationMaker::getTrackDetaDphi(float track_eta, float track_phi, int id){
00639         float eta_tower, phi_tower;
00640         eta_tower = phi_tower = -999.;
00641         mEmcGeom->getEtaPhi(id,eta_tower,phi_tower);
00642                                         
00643         //now calculate distance from center of tower:
00644         float dphi = phi_tower - track_phi;
00645         while(dphi >  M_PI)             {dphi -= 2.0 * M_PI;}
00646         while(dphi < -1.*M_PI)  {dphi += 2.0 * M_PI;}
00647 
00648         pair<float, float> dtow;
00649         dtow.first = eta_tower - track_eta;
00650         dtow.second = dphi;
00651         
00652         return dtow;
00653 }
00654 */
00655 
00656 pair<unsigned short, pair<float,float> > StEmcOfflineCalibrationMaker::getTrackTower(const StMuTrack* track, bool useExitRadius,int det){ //1=BTOW, 3=BSMDE, 4=BSMDP
00657         pair<unsigned short, pair<float,float> > tower;
00658         tower.first = 0;
00659         tower.second.first = 1000.;
00660         tower.second.second = 1000.;
00661 
00662         StThreeVectorD momentum,position;
00663         Double_t radius;
00664         if(det==1) radius = mEmcGeom->Radius();
00665         if(det==2) radius = mEmcGeom->Radius();
00666         if(det==3) radius = mSmdEGeom->Radius();
00667         if(det==4) radius = mSmdPGeom->Radius();
00668                 
00669         const StEventSummary& evtSummary = muDstMaker->muDst()->event()->eventSummary();
00670         Double_t mField = evtSummary.magneticField()/10;
00671                 
00672         //add 30 cm to radius to find out if track left same tower
00673         if(useExitRadius) radius += 30.0;
00674         
00675         bool goodProjection = mEmcPosition->trackOnEmc(&position,&momentum,track,mField,radius);
00676         if(goodProjection){
00677                 int m,e,s,id=0;
00678                 float eta=position.pseudoRapidity();
00679                 float phi=position.phi();
00680                 if(det==1){
00681                   mEmcGeom->getBin(phi,eta,m,e,s);
00682                   s = abs(s);
00683                 
00684                   if(mEmcGeom->getId(m,e,s,id)==0){
00685                     tower.first = id;
00686                     tower.second = getTrackDetaDphi(eta, phi, id, det);
00687                   }
00688                 }
00689                 else if(det==3){
00690                   int check=mSmdEGeom->getBin(phi,eta,m,e,s);
00691                   if(!check){
00692                     s = abs(s);
00693                     if(mSmdEGeom->getId(m,e,s,id)==0){
00694                       tower.first = id;
00695                       tower.second = getTrackDetaDphi(eta, phi, id, det);
00696                     }
00697                   }
00698                 }
00699                 else if(det==4){
00700                   int check=mSmdPGeom->getBin(phi,eta,m,e,s);
00701                   s = abs(s);
00702                   if(!check){
00703                     if(mSmdPGeom->getId(m,e,s,id)==0){
00704                       tower.first = id;
00705                       tower.second = getTrackDetaDphi(eta, phi, id, det);
00706                     }
00707                   }
00708                 }
00709         }
00710         
00711         return tower;
00712 }
00713 
00714 pair<float,float> StEmcOfflineCalibrationMaker::getTrackDetaDphi(float track_eta, float track_phi, int id, int det){
00715         float eta_tower, phi_tower;
00716         eta_tower = phi_tower = -999.;
00717         pair<float, float> dtow;
00718         if(det==1){
00719           mEmcGeom->getEtaPhi(id,eta_tower,phi_tower);
00720                                         
00721           //now calculate distance from center of tower:
00722           float dphi = phi_tower - track_phi;
00723           while(dphi >  M_PI)           {dphi -= 2.0 * M_PI;}
00724           while(dphi < -1.*M_PI)        {dphi += 2.0 * M_PI;}
00725 
00726           dtow.first = eta_tower - track_eta;
00727           dtow.second = dphi;
00728         }
00729         if(det==3){
00730           mSmdEGeom->getEtaPhi(id,eta_tower,phi_tower);
00731                                         
00732           //now calculate distance from center of strip:
00733           float dphi = phi_tower - track_phi;
00734           while(dphi >  M_PI)           {dphi -= 2.0 * M_PI;}
00735           while(dphi < -1.*M_PI)        {dphi += 2.0 * M_PI;}
00736 
00737           dtow.first = eta_tower - track_eta;
00738           dtow.second = dphi;
00739         }
00740         if(det==4){
00741           mSmdPGeom->getEtaPhi(id,eta_tower,phi_tower);
00742                                         
00743           //now calculate distance from center of strip:
00744           float dphi = phi_tower - track_phi;
00745           while(dphi >  M_PI)           {dphi -= 2.0 * M_PI;}
00746           while(dphi < -1.*M_PI)        {dphi += 2.0 * M_PI;}
00747 
00748           dtow.first = eta_tower - track_eta;
00749           dtow.second = dphi;
00750         }
00751         
00752         return dtow;
00753 }
00754 
00755 double StEmcOfflineCalibrationMaker::highestNeighbor(int id){
00756         double nSigma=0.;
00757         
00758         for(int deta=-1; deta<=1; deta++){
00759                 for(int dphi=-1; dphi<=1; dphi++){
00760                         if(deta==0 && dphi==0) continue;
00761                         int nextId = mEmcPosition->getNextTowerId(id,deta,dphi);
00762                         if(nextId<1 || nextId>4800) continue;
00763                         if((mADC[BTOW-1][nextId-1]-mPedestal[BTOW-1][nextId-1]) > nSigma*mPedRMS[BTOW-1][nextId-1]){
00764                                 nSigma = (mADC[BTOW-1][nextId-1] - mPedestal[BTOW-1][nextId-1])/mPedRMS[BTOW-1][nextId-1];
00765                         }
00766                 }
00767         }
00768         
00769         return nSigma;
00770 }
00771 
00772 void StEmcOfflineCalibrationMaker::addMinBiasTrigger(unsigned int trigId){ 
00773         mbTriggers.push_back(trigId);
00774 }
00775 
00776 void StEmcOfflineCalibrationMaker::addHighTowerTrigger(unsigned int trigId){ 
00777         htTriggers.push_back(trigId);
00778 }
00779 void StEmcOfflineCalibrationMaker::addHTTPTrigger(unsigned int trigId){
00780   httpTriggers.push_back(trigId);
00781 }
00782 void StEmcOfflineCalibrationMaker::addFastTrigger(unsigned int trigId){
00783   fastTriggers.push_back(trigId);
00784 }