StRoot: StTpcHitMaker/StTpcRTSHitMaker.cxx Source File

00001 /***************************************************************************
00002  *
00003  * $Id: StTpcRTSHitMaker.cxx,v 1.30 2012/05/17 20:05:59 fisyak Exp $
00004  *
00005  * Author: Valeri Fine, BNL Feb 2007
00006  ***************************************************************************
00007  *
00008  * Description:  Make clusters from StTpcRawData and fill the StEvent      */
00009 #include <assert.h>
00010 #include <sys/types.h>
00011 #include <stdio.h>
00012 #include "StTpcHitMaker.h"
00013 #include "StTpcRTSHitMaker.h"
00014 
00015 #include "TString.h"
00016 
00017 #include "StTpcRawData.h"
00018 #include "StEvent/StTpcRawData.h"
00019 #include "StEvent/StTpcHit.h"
00020 #include "StEvent/StEvent.h"
00021 #include "StEvent/StTpcHitCollection.h"
00022 #include "StThreeVectorF.hh"
00023 #include "StTpcDb/StTpcDb.h"
00024 #include "StDbUtilities/StCoordinates.hh"
00025 #include "StDetectorDbMaker/St_tss_tssparC.h"
00026 #include "StDetectorDbMaker/St_tpcPadGainT0C.h"
00027 #include "StDetectorDbMaker/St_tpcAnodeHVavgC.h"
00028 #include "StDetectorDbMaker/St_tpcMaxHitsC.h"
00029 #include "StDetectorDbMaker/StDetectorDbTpcRDOMasks.h"
00030 #include "StDetectorDbMaker/St_tpcPadPlanesC.h"
00031 #include "StMessMgr.h" 
00032 
00033 #ifndef NEW_DAQ_READER
00034 #  include "RTS/include/rtsLog.h"
00035 #  include "RTS/src/RTS_READER/rts_reader.h"
00036 #  include "RTS/src/RTS_READER/daq_dta.h"
00037 #  include "RTS/src/DAQ_TPX/daq_tpx.h"
00038 #else /* NEW_DAQ_READER */
00039 #  include "StDAQMaker/StDAQReader.h"
00040 #  include "StRtsTable.h"
00041 #  include "DAQ_TPX/daq_tpx.h"
00042 #  include "DAQ_READER/daq_dta.h"
00043 #  include "DAQ_READER/daqReader.h"
00044 #endif /* NEW_DAQ_READER */
00045 ClassImp(StTpcRTSHitMaker); 
00046 #define __DEBUG__
00047 #ifdef __DEBUG__
00048 #define PrPP(A,B) if (Debug()%10 > 1) {LOG_INFO << "StTpcRTSHitMaker::" << (#A) << "\t" << (#B) << " = \t" << (B) << endm;}
00049 #else
00050 #define PrPP(A,B)
00051 #endif
00052 //________________________________________________________________________________
00053 StTpcRTSHitMaker::~StTpcRTSHitMaker() {
00054   SafeDelete(fTpx);
00055   if (mTpx_RowLen) delete [] mTpx_RowLen;
00056 }
00057 //________________________________________________________________________________
00058 Int_t StTpcRTSHitMaker::Init() {
00059   memset(maxHits,0,sizeof(maxHits));
00060   maxBin0Hits = 0;
00061   bin0Hits = 0;
00062   return kStOK;
00063 }
00064 //________________________________________________________________________________
00065 Int_t StTpcRTSHitMaker::InitRun(Int_t runnumber) {
00066   SetAttr("minSector",1);
00067   SetAttr("maxSector",24);
00068   SetAttr("minRow",1);
00069   Int_t NoRowsInner = St_tpcPadPlanesC::instance()->innerPadRows();
00070   Int_t NoRowsOuter = St_tpcPadPlanesC::instance()->outerPadRows();
00071   NoRows = NoRowsInner + NoRowsOuter;
00072   if (NoRows != 45) {
00073     // Fill no. of pad per row 
00074     mTpx_RowLen = new UChar_t[NoRows+1];
00075     mTpx_RowLen[0] = 0;
00076     for (Int_t i = 1; i <= NoRows; i++) {
00077       mTpx_RowLen[i] = St_tpcPadPlanesC::instance()->padsPerRow(i);
00078     }
00079   }
00080   SetAttr("maxRow",NoRows);
00081   SafeDelete(fTpx);
00082   fTpx = new daq_tpx() ; 
00083   if (GetDate() >= 20091215) fTpx->fcf_run_compatibility = 10 ;
00084   if (GetDate() <= 20090101) fminCharge = 40;
00085   // do gains example; one loads them from database but I don't know how...
00086   if (NoRows <= 45) { // hack for now take Tonko's defaults for iTpx
00087     daq_dta *dta  = fTpx->put("gain");
00088     
00089     // Prepare scaled hit maxima
00090     
00091     // No hit maxima if these DB params are 0
00092     Int_t maxHitsPerSector = St_tpcMaxHitsC::instance()->maxSectorHits();
00093     Int_t maxBinZeroHits = St_tpcMaxHitsC::instance()->maxBinZeroHits();
00094     Int_t livePads = 0;
00095     Int_t totalPads = 0;
00096     Float_t liveFrac = 1;
00097     for(Int_t sector=1;sector<=24;sector++) {
00098       Int_t liveSecPads = 0;
00099       Int_t totalSecPads = 0;
00100       for(Int_t row=1;row<=NoRows;row++) {
00101         Int_t numPadsAtRow = St_tpcPadPlanesC::instance()->padsPerRow(row);
00102         daq_det_gain *gain = (daq_det_gain *) dta->request(183);        // max pad+1            
00103         assert(gain);
00104         gain[0].gain = 0.0;     // kill pad0 just in case..
00105         gain[0].t0   = 0.0;
00106         for(Int_t pad = 1; pad <= numPadsAtRow; pad++) {
00107           if (m_Mode == 2) {
00108             if (St_tpcPadGainT0C::instance()->Gain(sector,row,pad) > 0) 
00109               gain[pad].gain = 1.;
00110             else 
00111               gain[pad].gain = .0;
00112             gain[pad].t0   = 0.;
00113           } else {
00114             if (St_tpcPadGainT0C::instance()->Gain(sector,row,pad) <= 0) continue;
00115             gain[pad].gain = St_tpcPadGainT0C::instance()->Gain(sector,row,pad);
00116             gain[pad].t0   = St_tpcPadGainT0C::instance()->T0(sector,row,pad);
00117           }
00118         }
00119         dta->finalize(183,sector,row);
00120         if (maxHitsPerSector > 0 || maxBinZeroHits > 0) {
00121           totalSecPads += numPadsAtRow;
00122           if (StDetectorDbTpcRDOMasks::instance()->isOn(sector,
00123             StDetectorDbTpcRDOMasks::instance()->rdoForPadrow(row)) &&
00124               St_tpcAnodeHVavgC::instance()->livePadrow(sector,row))
00125             liveSecPads += numPadsAtRow;
00126         }
00127       }
00128       livePads += liveSecPads;
00129       totalPads += totalSecPads;
00130       if (maxHitsPerSector > 0) {
00131         liveFrac = TMath::Max((Float_t) 0.1,
00132                               ((Float_t) liveSecPads) / ((Float_t) totalSecPads));
00133         maxHits[sector-1] = (Int_t) (liveFrac * maxHitsPerSector);
00134         if (Debug()) {LOG_INFO << "maxHits in sector " << sector
00135                                << " = " << maxHits[sector-1] << endm;}
00136       } else {
00137         maxHits[sector-1] = 0;
00138         if (Debug()) {LOG_INFO << "No maxHits in sector " << sector << endm;}
00139       }
00140     }
00141     if (maxBinZeroHits > 0) {
00142       liveFrac = TMath::Max((Float_t) 0.1,
00143                             ((Float_t) livePads) / ((Float_t) totalPads));
00144       maxBin0Hits = (Int_t) (liveFrac * maxBinZeroHits);
00145       if (Debug()) {LOG_INFO << "maxBinZeroHits " << maxBin0Hits << endm;}
00146     } else {
00147       maxBin0Hits = 0;
00148       if (Debug()) {LOG_INFO << "No maxBinZeroHits" << endm;}
00149     }
00150   }
00151   /*
00152     InitRun will setup the internal representations of gain 
00153     and other necessary structures but if the gains have not
00154     been previously loaded as shown in the example above they
00155     will be set to 1.0!
00156   */
00157   fTpx->InitRun(runnumber);
00158   return kStOK;
00159 }
00160 //________________________________________________________________________________
00161 Int_t StTpcRTSHitMaker::Make() {
00162   static  Short_t ADCs[__MaxNumberOfTimeBins__];
00163   static UShort_t IDTs[__MaxNumberOfTimeBins__];
00164   StEvent*   rEvent      = (StEvent*)    GetInputDS("StEvent");
00165   if (! rEvent) {
00166     LOG_WARN << "There is no StEvent" << endm;
00167     return kStWarn;
00168   }
00169   StTpcHitCollection *hitCollection = rEvent->tpcHitCollection();
00170   TDataSet*  tpcRawEvent =               GetInputDS("Event");
00171   if (! tpcRawEvent) {
00172     LOG_WARN << "There is not Tpc Raw Event" << endm;
00173     return kStWarn;
00174   }
00175   StTpcRawData *tpcRawData = (StTpcRawData *) tpcRawEvent->GetObject();
00176   if (! tpcRawData) {
00177     LOG_WARN << "There is not Tpc Raw Data" << endm;
00178     return kStWarn;
00179   }
00180   // create (or reuse) the adc_sim bank...
00181   // add a bunch of adc data for a specific sector:row:pad
00182   Int_t minSector = IAttr("minSector");
00183   Int_t maxSector = IAttr("maxSector");
00184   Int_t minRow    = IAttr("minRow");
00185   Int_t maxRow    = IAttr("maxRow");
00186 
00187   bin0Hits = 0;
00188   daq_dta *dta = 0;
00189   for (Int_t sec = minSector; sec <= maxSector; sec++) {
00190     StTpcDigitalSector *digitalSector = tpcRawData->GetSector(sec);
00191     if (! digitalSector) continue;
00192     UShort_t Id = 0;
00193     if (NoRows != 45) dta = fTpx->put("adc_sim",0,NoRows+1,0,mTpx_RowLen); // used for any kind of data; transparent pointer
00194     else              dta = fTpx->put("adc_sim");
00195     Int_t hitsAdded = 0;
00196     Int_t nup = 0;
00197     Int_t NoAdcs = 0;
00198     for (Int_t row = minRow; row <= maxRow; row++) {
00199       if (! St_tpcPadGainT0C::instance()->livePadrow(sec,row)) continue;
00200       Int_t Npads = digitalSector->numberOfPadsInRow(row);
00201       if (! Npads) continue;
00202       for(Int_t pad = 1; pad <= Npads; pad++) {
00203         UInt_t ntimebins = digitalSector->numberOfTimeBins(row,pad);
00204         if (! ntimebins) continue;
00205         // allocate space for at least 512 pixels (timebins)
00206         daq_sim_adc_tb *d = (daq_sim_adc_tb *) dta->request(__MaxNumberOfTimeBins__);
00207         // add adc data for this specific sector:row:pad
00208         memset (ADCs, 0, sizeof(ADCs));
00209         memset (IDTs, 0, sizeof(IDTs));
00210         digitalSector->getTimeAdc(row,pad,ADCs,IDTs);
00211         UInt_t l = 0;
00212         for (UInt_t k = 0; k < __MaxNumberOfTimeBins__; k++) {
00213           if (ADCs[k]) {
00214             d[l].adc = ADCs[k];
00215             d[l].tb  = k;
00216             d[l].track_id = IDTs[k];
00217             l++;
00218           }
00219         }
00220         if (l > 0) {
00221           dta->finalize(l,sec,row,pad);
00222           NoAdcs += l;
00223         }
00224       }
00225     }      
00226     if (! NoAdcs) continue;
00227     if (Debug() > 0) {
00228       // verify data!
00229       dta = fTpx->get("adc_sim");
00230       while(dta && dta->iterate()) {
00231         LOG_INFO << Form("*** sec %2d, row %2d, pad %3d: %3d pixels",dta->sec,dta->row,dta->pad,dta->ncontent) << endm;
00232         for(UInt_t i=0;i<dta->ncontent;i++) {
00233           if (Debug() > 1 || dta->sim_adc[i].track_id) {
00234             LOG_INFO << Form("    %2d: adc %4d, tb %3d: track %4d",i,
00235                              dta->sim_adc[i].adc,
00236                              dta->sim_adc[i].tb,
00237                              dta->sim_adc[i].track_id
00238                              ) << endm;
00239           }
00240         }
00241       }
00242     }
00243     static StTpcCoordinateTransform transform(gStTpcDb);
00244     static StThreeVectorF hard_coded_errors;
00245     //      fTpx->put("cld_sim");       // clean up clusters
00246     dta = fTpx->get("cld_sim"); // rerun the cluster finder on the simulated data...
00247     Double_t ADC2GeV = 0;
00248     Int_t rowOld = -1;
00249     static Int_t iBreak = 0;
00250     while(dta && dta->iterate()) {
00251       if (Debug() > 0) {
00252         LOG_INFO << Form("CLD sec %2d: row %2d: %d clusters",dta->sec, dta->row, dta->ncontent) << endm;
00253       }
00254       for(UInt_t i=0;i<dta->ncontent;i++) {
00255         if (! dta->sim_cld[i].cld.pad) continue;
00256         if (dta->sim_cld[i].cld.tb >= __MaxNumberOfTimeBins__) continue;
00257         if (dta->sim_cld[i].cld.charge < fminCharge) continue;
00258         if (Debug() > 0) {
00259           if (Debug() > 1 || ( dta->sim_cld[i].cld.p2 - dta->sim_cld[i].cld.p1 <= 1 )) {
00260             LOG_INFO << Form("    pad %f[%d:%d], tb %f[%d:%d], cha %d, fla 0x%X, Id %d, Q %d ",
00261                              dta->sim_cld[i].cld.pad,
00262                              dta->sim_cld[i].cld.p1,
00263                              dta->sim_cld[i].cld.p2,
00264                              dta->sim_cld[i].cld.tb,
00265                              dta->sim_cld[i].cld.t1,
00266                              dta->sim_cld[i].cld.t2,
00267                              dta->sim_cld[i].cld.charge,
00268                              dta->sim_cld[i].cld.flags,
00269                              dta->sim_cld[i].track_id,
00270                              dta->sim_cld[i].quality
00271                              ) << endm;
00272             iBreak++;
00273           }
00274         }
00275         if (! hitCollection )  {
00276           hitCollection = new StTpcHitCollection();
00277           rEvent->setTpcHitCollection(hitCollection);
00278         }
00279         StTpcPadCoordinate Pad(dta->sec, dta->row, dta->sim_cld[i].cld.pad, dta->sim_cld[i].cld.tb); PrPP(Make,Pad);
00280         static StTpcLocalSectorCoordinate LS;
00281         static StTpcLocalCoordinate L;
00282         transform(Pad,LS,kFALSE,kTRUE); PrPP(Make,LS); // don't useT0, useTau                  
00283         transform(LS,L);                                                                             PrPP(Make,L);
00284         if (dta->row != rowOld) {
00285           rowOld = dta->row;
00286           Double_t gain = (dta->row<=13) ? St_tss_tssparC::instance()->gain_in() : St_tss_tssparC::instance()->gain_out();
00287           Double_t wire_coupling = (dta->row<=13) ? 
00288             St_tss_tssparC::instance()->wire_coupling_in() : 
00289             St_tss_tssparC::instance()->wire_coupling_out();
00290           ADC2GeV = ((Double_t) St_tss_tssparC::instance()->ave_ion_pot() * 
00291                      (Double_t) St_tss_tssparC::instance()->scale())/(gain*wire_coupling) ;
00292         }
00293         UInt_t hw = 1;   // detid_tpc
00294         hw += dta->sec << 4;     // (row/100 << 4);   // sector
00295         hw += dta->row << 9;     // (row%100 << 9);   // row
00296         
00297         Int_t npads = TMath::Abs(dta->sim_cld[i].cld.p2 - dta->sim_cld[i].cld.p1) + 1;
00298         hw += (npads   << 15);  // npads
00299         
00300         Int_t ntmbk = TMath::Abs(dta->sim_cld[i].cld.t2 - dta->sim_cld[i].cld.t1) + 1;
00301         hw += (ntmbk << 22);  // ntmbks...
00302         Double_t q = ADC2GeV*dta->sim_cld[i].cld.charge;
00303         Id++;
00304         StTpcHit *hit = StTpcHitMaker::StTpcHitFlag(L.position(),hard_coded_errors,hw,q
00305                                                     , (UChar_t ) 0  // counter 
00306                                                     , (UShort_t) dta->sim_cld[i].track_id  // idTruth=0
00307                                                     , (UShort_t) dta->sim_cld[i].quality   // quality=0,
00308                                                     , Id                                   // id =0,
00309                                                     , dta->sim_cld[i].cld.p1 //  mnpad
00310                                                     , dta->sim_cld[i].cld.p2 //  mxpad
00311                                                     , dta->sim_cld[i].cld.t1 //  mntmbk
00312                                                     , dta->sim_cld[i].cld.t2 //  mxtmbk
00313                                                     , dta->sim_cld[i].cld.pad
00314                                                     , dta->sim_cld[i].cld.tb 
00315                                                     , dta->sim_cld[i].cld.charge
00316                                                     , dta->sim_cld[i].cld.flags);
00317         /*tpxFCF.h
00318           #define FCF_ROW_EDGE            16      // 0x10 touched end of row
00319           #define FCF_BROKEN_EDGE         32      // 0x20 touches one of the mezzanine edges
00320           #define FCF_DEAD_EDGE           64      // 0x40 touches a dead pad 
00321         */
00322         assert(dta->sim_cld[i].cld.pad >  0 && dta->sim_cld[i].cld.pad <= 182 && 
00323                dta->sim_cld[i].cld.tb  >= 0 && dta->sim_cld[i].cld.tb  <  512);
00324         hitsAdded++;
00325         if (hit->minTmbk() == 0) bin0Hits++;
00326         hitCollection->addHit(hit);
00327       }
00328     }
00329     // Set IdTruth
00330     for (Int_t row = minRow; row <= maxRow; row++) {
00331       Int_t Npads = digitalSector->numberOfPadsInRow(row);
00332       for(Int_t pad = 1; pad <= Npads; pad++) {
00333         UInt_t ntimebins = digitalSector->numberOfTimeBins(row,pad);
00334         if (! ntimebins) continue;
00335         digitalSector->getTimeAdc(row,pad,ADCs,IDTs);
00336         // Update pixels if any (for data)
00337         dta = fTpx->get("adc_sim",sec);
00338         Int_t Updated = 0;
00339         while(dta && dta->iterate()) {
00340           Int_t secC  = dta->sec;
00341           Int_t rowC  = dta->row;
00342           Int_t padC  = dta->pad;
00343           if (secC != sec || rowC != row || padC != pad) continue;
00344           for(UInt_t i=0;i<dta->ncontent;i++) {
00345             Int_t tb = dta->sim_adc[i].tb;
00346             if (ADCs[tb] != dta->sim_adc[i].adc ||
00347                 IDTs[tb] != dta->sim_adc[i].track_id) {
00348               if (Debug() > 1) {
00349                 LOG_INFO << Form("tb %3d    adc %4d => %4d, track %4d => %4d",tb,
00350                                  ADCs[tb],dta->sim_adc[i].adc,
00351                                  IDTs[tb],dta->sim_adc[i].track_id) << endm;
00352               }
00353               ADCs[tb] = dta->sim_adc[i].adc;
00354               IDTs[tb] = dta->sim_adc[i].track_id;
00355               Updated++;
00356             }
00357           }
00358         }
00359         if (Updated) {
00360           digitalSector->putTimeAdc(row,pad,ADCs,IDTs);
00361           nup += Updated;
00362         }
00363       }
00364       if (nup && Debug() > 1) {
00365         LOG_INFO << "Update total " << nup << " pixels from Sector / row = " << sec << " / " << row <<   endm;
00366       }
00367     }
00368     if (maxHits[sec-1] && hitsAdded > maxHits[sec-1]) {
00369       LOG_ERROR << "Too many hits (" << hitsAdded << ") in one sector ("
00370                 << sec << "). Skipping event." << endm;
00371       return kStSkip;
00372     }
00373   }
00374   if (maxBin0Hits && bin0Hits > maxBin0Hits) {
00375     LOG_ERROR << "Too many hits (" << bin0Hits
00376               << ") starting at time bin 0. Skipping event." << endm;
00377     return kStSkip;
00378   }
00379   StTpcHitMaker::AfterBurner(hitCollection);
00380   return kStOK;
00381 }