StRoot: StSvtSeqAdjMaker/StSvtSeqAdjMaker.cxx Source File

00001  /***************************************************************************
00002  *
00003  * $Id: StSvtSeqAdjMaker.cxx,v 1.56 2007/04/28 17:57:10 perev Exp $
00004  *
00005  * Author: 
00006  ***************************************************************************
00007  *
00008  * Description: Sequence
00009  **************************************************************************
00010  *
00011  * 
00012  * Revision 1.29  2001/09/26 18:42:48  Takahashi
00013  * Added new bad anode list and switched ON the bad anode elimination
00014  *
00015  * $Log: StSvtSeqAdjMaker.cxx,v $
00016  * Revision 1.56  2007/04/28 17:57:10  perev
00017  * Redundant StChain.h removed
00018  *
00019  * Revision 1.55  2005/07/23 03:37:34  perev
00020  * IdTruth + Cleanup
00021  *
00022  * Revision 1.54  2004/04/03 01:17:25  caines
00023  * I think this time I should really have stopped valgrind complaing lets see
00024  *
00025  * Revision 1.53  2004/03/18 04:05:01  caines
00026  * Remove from global scope variables used in debug mode as they shouldnt be there and caused erratic behaviour, also initialise some variables that valgrind was complaining about - didnt really need to as they are sent back from function which initialises them properly always but doesnt hurt
00027  *
00028  * Revision 1.52  2004/02/11 00:42:17  caines
00029  * Put common mode back
00030  *
00031  * Revision 1.51  2004/01/27 02:44:16  perev
00032  * LeakOff
00033  *
00034  * Revision 1.50  2004/01/22 17:15:21  caines
00035  * Remove gain calibration files and common mode moise subtraction
00036  *
00037  * Revision 1.49  2003/09/18 18:16:44  caines
00038  * Initialise stuff for new compiler
00039  *
00040  * Revision 1.48  2003/07/31 19:01:14  caines
00041  * Make changes so the slow simulator can run
00042  *
00043  * Revision 1.47  2003/07/18 17:15:40  caines
00044  * Fix Pedoffset to be 20 not 10 change variables to int from floats to avoid casting problems, fix that when pedestal goes negative we dont
00045  *
00046  * Revision 1.46  2003/04/14 15:56:59  munhoz
00047  * changes for gain calibration file
00048  *
00049  * Revision 1.45  2003/01/28 20:21:17  munhoz
00050  * changing from Init() to InitRun()
00051  *
00052  * Revision 1.44  2003/01/21 01:27:44  jeromel
00053  * hfile->write(0 while NULL caused spurious crash.
00054  * Took the oportunity to add GetCVS()
00055  * Some maniaco-compulsive //! alignement fixes ...
00056  *
00057  * Revision 1.43  2002/09/26 20:46:06  caines
00058  * Make sure initRun works when no svt data present and correct code to allow for the ordering init and initrun are called in
00059  *
00060  * Revision 1.42  2002/09/20 19:35:19  caines
00061  * Change building of file name
00062  *
00063  * Revision 1.41  2002/09/19 16:17:48  caines
00064  * Add code to do Juns gain calibration
00065  *
00066  * Revision 1.40  2002/05/09 16:55:08  munhoz
00067  * add reading bad anodes from DB
00068  *
00069  * Revision 1.39  2002/04/21 20:28:01  caines
00070  * Put some of the print out only if in debug mode
00071  *
00072  * Revision 1.37  2002/01/11 22:49:15  caines
00073  * Fix sequence merging bugs-hopefully
00074  *
00075  * Revision 1.36  2001/12/13 03:08:11  caines
00076  * Can now subtract common mode noise via black anodes 239 and 2
00077  *
00078  * Revision 1.35  2001/10/19 23:31:30  caines
00079  * Correct problem that if anodes were missing didnt do average common mode noise calc
00080  *
00081  * Revision 1.32  2001/10/04 02:56:01  caines
00082  * Change default pedoffset and black anodes to current values
00083  *
00084  * Revision 1.31  2001/10/02 22:55:57  caines
00085  * Jun was coorect about B2L8D2H2
00086  *
00087  * Revision 1.30  2001/09/28 20:47:38  caines
00088  * Fix typo in bad anode listing
00089  *
00090  * Revision 1.29  2001/09/28 15:36:41  caines
00091  * Add in bad anode elimination
00092  *
00093  * Revision 1.28  2001/09/26 18:42:48  caines
00094  * Fix 2 anode subtraction routines
00095  *
00096  * Revision 1.27  2001/09/22 00:35:54  caines
00097  * Fixes now that AddData() is cleared everyevent
00098  *
00099  * Revision 1.26  2001/09/16 22:24:22  caines
00100  * Fix for when SVT isnt in every event
00101  *
00102  * Revision 1.25  2001/08/28 19:02:55  caines
00103  * Remove hardwired line that makes it always do common mode average (dumb!)
00104  *
00105  * Revision 1.24  2001/08/24 21:40:04  caines
00106  * Adjust rawHybridData not SeqHybridData in Common Mode noise subtraction
00107  *
00108  * Revision 1.23  2001/08/24 20:57:45  caines
00109  * Do common mode noise suppression from first two anodes
00110  *
00111  * Revision 1.22  2001/08/22 14:24:49  caines
00112  * Raise m_thresh_hi to 10
00113  *
00114  * Revision 1.21  2001/07/25 14:47:41  caines
00115  * Fix filling histogram only when debug is on
00116  *
00117  * Revision 1.20  2001/07/22 20:31:28  caines
00118  * Better tuning for real data. Common mode noise calc and sub. Avoid overlapping seq. Fill histograms only in debug
00119  *
00120  * Revision 1.19  2001/05/06 22:16:35  caines
00121  * Fix pointer problem for DAQ Ped files
00122  *
00123  * Revision 1.18  2001/04/30 22:27:40  caines
00124  * Fix adjusting for ZSP data
00125  *
00126  * Revision 1.17  2001/04/25 19:00:17  perev
00127  * HPcorrs
00128  *
00129  * Revision 1.16  2001/03/28 21:49:38  caines
00130  * Correct bug that allowed AdjustSeq1 code to fall off end
00131  *
00132  * Revision 1.15  2001/03/22 20:47:10  caines
00133  * Comment out some of the QA histograms
00134  *
00135  * Revision 1.14  2001/02/07 19:15:05  caines
00136  * Change char[3] to char[4] for full SVT running
00137  *
00138  * Revision 1.13  2000/11/30 20:45:56  caines
00139  * Dynamically calc prob values, use database
00140  *
00141  * Revision 1.12  2000/10/31 16:21:59  caines
00142  * Set up efaults for paramters in init
00143  *
00144  * Revision 1.11  2000/10/02 13:48:10  caines
00145  * Adjusting donw hybrid by hybrid
00146  *
00147  * Revision 1.10  2000/09/14 22:13:56  caines
00148  * Move histogram creation to better place
00149  *
00150  * Revision 1.9  2000/08/28 22:11:38  caines
00151  * Fixed check that data exists before using it
00152  *
00153  * Revision 1.8  2000/08/24 04:23:50  caines
00154  * Improved histograms
00155  *
00156  * Revision 1.7  2000/08/21 12:57:31  caines
00157  * Now opens and reads in ped using CalibMaker
00158  *
00159  * Revision 1.6  2000/08/09 02:05:08  caines
00160  * Only add pixels in the ASIC-like sequence adjusting
00161  *
00162  * Revision 1.5  2000/07/16 22:32:23  caines
00163  * Now also saves RAW data
00164  *
00165  * Revision 1.4  2000/07/13 14:50:13  caines
00166  * Fixed up problem when pointer went past 128
00167  *
00168  * Revision 1.3  2000/07/11 18:36:15  caines
00169  * Updates to save more of sequence for fitting
00170  *
00171  * Revision 1.2  2000/07/03 02:07:56  perev
00172  * StEvent: vector<TObject*>
00173  *
00174  * Revision 1.1  2000/06/15 20:04:54  caines
00175  * Initial versions of sequence adjusting codes
00176  *
00177  **************************************************************************/
00178 
00179 #include "TFile.h"
00180 #include "TH1.h"
00181 #include "TH2.h"
00182 #include "StSequence.hh"
00183 #include "StMCTruth.h"
00184 #include "StSvtClassLibrary/StSvtData.hh"
00185 #include "StSvtClassLibrary/StSvtHybridData.hh"
00186 #include "StSvtClassLibrary/StSvtHybridPed.hh"
00187 #include "StSvtClassLibrary/StSvtHybridCollection.hh"
00188 #include "StSvtInverseProducts.hh"
00189 #include "StSvtProbValues.hh"
00190 #include "StSvtPedSub.h"
00191 #include "StSvtClassLibrary/StSvtHybridBadAnodes.hh"
00192 #include "StSvtSeqAdjMaker.h"
00193 #include "StSvtCalibMaker/StSvtPedMaker.h"
00194 #include "StMessMgr.h"
00195 #include "StIOMaker/StIOMaker.h"
00196 
00197 #include "St_DataSetIter.h"
00198 #include "TObjectSet.h"
00199 
00200 
00201 
00202 //___________________________________________________________________________________________
00203 StSvtSeqAdjMaker::StSvtSeqAdjMaker(const char *name) : StMaker(name)
00204 {
00205 
00206   mSvtDataSet = NULL;
00207   mSvtAdjData = NULL;
00208   mSvtRawData = NULL;
00209   mSvtBadAnodes = NULL;
00210   mHybridRawData = NULL;
00211   mHybridAdjData = NULL;
00212   mSvtPedSub = NULL;
00213   mSvtPedColl = NULL;
00214   anolist = NULL;
00215   mInvProd = NULL;
00216   mProbValue = NULL;
00217 
00218   // Set up some defaults
00219 
00220   mPedFile = NULL;
00221   mPedOffSet = 20;
00222   m_thresh_lo = 3+mPedOffSet;
00223   m_thresh_hi = 5+mPedOffSet; 
00224   m_n_seq_lo  = 2;
00225   m_n_seq_hi  = 0;
00226   m_inv_prod_lo = 0;
00227  
00228 }
00229 
00230 //____________________________________________________________________________
00231 StSvtSeqAdjMaker::~StSvtSeqAdjMaker(){
00232 }
00233 //____________________________________________________________________________
00234 Int_t StSvtSeqAdjMaker::Init(){
00235 
00236 
00237   hfile = NULL;
00238   GetSvtRawData();
00239   
00240   SetSvtData();
00241   
00242   mInvProd = new StSvtInverseProducts();
00243   mProbValue = new StSvtProbValues();
00244   
00245   float sigma;
00246   // set probability values
00247   for (int barrel = 1;barrel <= mSvtRawData->getNumberOfBarrels();barrel++) {
00248     for (int ladder = 1;ladder <= mSvtRawData->getNumberOfLadders(barrel);ladder++) {
00249       for (int wafer = 1;wafer <= mSvtRawData->getNumberOfWafers(barrel);wafer++) {
00250         for (int hybrid = 1;hybrid <= mSvtRawData->getNumberOfHybrids();hybrid++) {
00251             int index = mSvtRawData->getHybridIndex(barrel,ladder,wafer,hybrid);
00252             if(index < 0) continue;
00253             if (mSvtPedColl && mSvtPedColl->at(index))
00254               sigma = ((StSvtHybridPed*)mSvtPedColl->at(index))->getRMS();
00255             else
00256               sigma = mProbValue->GetSigma();
00257             mProbValue->SetProbValue(sigma);    
00258             mInvProd->SetProbTable(mProbValue);
00259         }
00260       }
00261     }
00262   }
00263   
00264   return  StMaker::Init();  
00265 }
00266 //____________________________________________________________________________
00267 Int_t StSvtSeqAdjMaker::InitRun( int runnumber)
00268 {
00269 
00270   if(GetSvtRawData() ){
00271     
00272     St_DataSet *dataSet = NULL;
00273     dataSet = GetDataSet("StSvtConfig");
00274     
00275     if (!dataSet) return kStWarn;
00276 
00277     mSvtRawData = new StSvtData((StSvtConfig*)(dataSet->GetObject()));
00278     AddData(new TObjectSet("StSvtRawData",mSvtRawData));
00279   }
00280 
00281   mTotalNumberOfHybrids = mSvtRawData->getTotalNumberOfHybrids();
00282 
00283   if( hfile){
00284     hfile->Write();
00285     hfile->Close();
00286   }
00287   
00288   string ioMakerFileName;
00289   string FileName;
00290   string DirName;
00291   DirName=".";
00292 
00293   StIOMaker* mIOMaker = (StIOMaker*)GetMaker("inputStream");
00294   cout << "************************" << mIOMaker << endl;
00295   if (mIOMaker){
00296     ioMakerFileName = string(mIOMaker->GetFile()); 
00297     FileName = buildFileName( DirName+"/", baseName(ioMakerFileName),".SvtGainCal.root"); 
00298     // cout << "Heres my name: " << FileName << endl;
00299   }else
00300     {
00301       FileName=buildFileName( DirName+"/", "SeqAdj-Run"+GetRunNumber(),".SvtGainCal.root"); 
00302     }
00303   //hfile  = new TFile(FileName.c_str(),"RECREATE","Demo ROOT file");
00304   CreateHist(mTotalNumberOfHybrids);
00305 
00306   GetBadAnodes();
00307 
00308   GetSvtPedestals();
00309 
00310   gMessMgr->Info()<< " StSvtSeqAdjMaker-info:"<<endm;
00311   gMessMgr->Info() << "     PedOffSet = "<<mPedOffSet<<endm;
00312   gMessMgr->Info() << "    thresh_lo = "<<m_thresh_lo <<endm;
00313   gMessMgr->Info() << "     thresh_hi = "<<m_thresh_hi <<endm;
00314   gMessMgr->Info() << "     n_seq_lo  = "<<m_n_seq_lo <<endm;
00315   gMessMgr->Info() << "     n_seq_hi  = "<< m_n_seq_hi <<endm;
00316   gMessMgr->Info() << "    inv_prod_lo  = "<<m_inv_prod_lo <<endm;
00317 
00318   return  kStOK;  
00319 }
00320 
00321 //____________________________________________________________________________
00322 
00323 
00324 Int_t StSvtSeqAdjMaker::GetSvtRawData()
00325 {
00326   St_DataSet *dataSet;
00327   
00328   dataSet = GetDataSet("StSvtRawData");
00329   if( !dataSet) {
00330     gMessMgr->Warning() << " No Svt Raw data set" << endm;
00331     return kStWarn;
00332   }
00333 
00334   mSvtRawData = (StSvtData*)(dataSet->GetObject());
00335   if( !mSvtRawData) {
00336     gMessMgr->Warning() << " No Svt Raw data " << endm;
00337     return kStWarn;
00338   }
00339   return kStOK;
00340 }
00341 
00342 //____________________________________________________________________________
00343 
00344 
00345 Int_t StSvtSeqAdjMaker::GetSvtPedestals()
00346 {
00347   if (mPedFile) {
00348     StSvtPedMaker* sdaq = (StSvtPedMaker*)GetMaker("SvtPed");
00349 
00350     if( sdaq)
00351       sdaq->ReadFromFile( mPedFile);
00352   }
00353 
00354   St_DataSet *dataSet;
00355   
00356   dataSet = GetDataSet("StSvtPedestal");
00357 
00358   // get pedestals and set pedestal subtraction ONLY if there is a pedestal data set
00359   if (dataSet) {
00360     mSvtPedColl = (StSvtHybridCollection*)(dataSet->GetObject());
00361     mSvtPedSub = new StSvtPedSub(mSvtPedColl);
00362   
00363     //Set offset
00364     mSvtRawData->setPedOffset(mPedOffSet);
00365   }
00366 
00367   return kStOK;
00368 }
00369 
00370 //__________________________________________________________________________
00371 
00372 
00373 Int_t  StSvtSeqAdjMaker::SetSvtData()
00374 {
00375   mSvtDataSet = new St_ObjectSet("StSvtData");
00376   AddData(mSvtDataSet);  
00377 
00378   mSvtAdjData = new StSvtData(*mSvtRawData);
00379   //cout<<"mSvtAdjData  = "<<mSvtAdjData<<endl;
00380   mSvtDataSet->SetObject((TObject*)mSvtAdjData); 
00381   // assert(mSvtAdjData);
00382 
00383   return kStOK;
00384 }
00385 
00386 //__________________________________________________________________________
00387 
00388 Int_t StSvtSeqAdjMaker::SetMinAdcLevels( int MinAdc1,  int MinAbove1,
00389                                          int MinAdc2,   int MinAbove2, int PedOffset){
00390 
00391   m_thresh_lo = MinAdc1+PedOffset;
00392   m_thresh_hi = MinAdc2+PedOffset; 
00393   m_n_seq_lo  = MinAbove1;
00394   m_n_seq_hi  = MinAbove2;
00395 
00396   //if (PedOffset)
00397     mPedOffSet = PedOffset;
00398 
00399   return kStOK;
00400 }
00401 
00402 //____________________________________________________________________________
00403 
00404 Int_t StSvtSeqAdjMaker::SetPedestalFile(const char* pedFile)
00405 {
00406   if(Debug()) gMessMgr->Debug() << "opening file called " << pedFile << " " 
00407                                << GetName() << endm;  
00408 
00409   mPedFile = pedFile;
00410 
00411   return kStOK;
00412 }
00413 
00414 
00415 //____________________________________________________________________________
00416 
00417 Int_t StSvtSeqAdjMaker::SetLowInvProd(int LowInvProd)
00418 {
00419   m_inv_prod_lo = LowInvProd;
00420   
00421   return kStOK; 
00422 }
00423 
00424 //___________________________________________________________________________
00425 
00426 Int_t StSvtSeqAdjMaker::GetBadAnodes()
00427 {
00428 
00429   St_DataSet *dataSet;
00430   
00431   dataSet = GetDataSet("StSvtBadAnodes");
00432   if( !dataSet) {
00433     gMessMgr->Warning() << " No Svt Bad Anodes data set" << endm;
00434     return kStWarn;
00435   }
00436 
00437   mSvtBadAnodes = (StSvtHybridCollection*)(dataSet->GetObject());
00438   if( !mSvtBadAnodes) {
00439     gMessMgr->Warning() << " No Svt Bad Anodes data " << endm;
00440     return kStWarn;
00441   }
00442   return kStOK;
00443 }
00444 
00445 //_____________________________________________________________________________
00446 Int_t StSvtSeqAdjMaker::CreateHist(Int_t tNuOfHyb)
00447 {
00448    // Create Histograms
00449 
00450   
00451    //mInvProdSeqAdj= new TH1D*[tNuOfHyb];
00452    mRawAdc = new TH1F*[tNuOfHyb];
00453    mAdcAfter = new TH1F*[tNuOfHyb];
00454    //mTimeAn = new TH2F*[tNuOfHyb];
00455 
00456    char  RawTitle[25], AdcAfterTitle[25], TimeAnTitle[25];
00457    char  Index[4];
00458    char* RawTitle_cut;
00459    char* AdcAf_cut;
00460    char* TimeAnCh;
00461   for (int barrel = 1;barrel <= mSvtRawData->getNumberOfBarrels();barrel++) {
00462     for (int ladder = 1;ladder <= mSvtRawData->getNumberOfLadders(barrel);ladder++) {
00463       for (int wafer = 1;wafer <= mSvtRawData->getNumberOfWafers(barrel);wafer++) {
00464         for (int hybrid = 1;hybrid <= mSvtRawData->getNumberOfHybrids();hybrid++) {
00465            
00466             int index = mSvtRawData->getHybridIndex(barrel,ladder,wafer,hybrid);
00467             if(index < 0) continue;
00468             
00469             sprintf(RawTitle,"RawAdcIn");
00470             sprintf(AdcAfterTitle,"numOfhitsIn");
00471             sprintf(TimeAnTitle,"TimeAn");
00472             sprintf(Index,"%d", index);
00473             RawTitle_cut = strcat(RawTitle,Index);
00474             AdcAf_cut = strcat(AdcAfterTitle,Index);
00475             TimeAnCh = strcat(TimeAnTitle,Index);
00476 
00477             mRawAdc[index] = new TH1F(RawTitle_cut,"Raw ADC for each anode",241,0,241);
00478             mAdcAfter[index] = new TH1F(AdcAf_cut,"Number of hits in each anode",241,0,241);
00479             //mTimeAn[index] = new TH2F(TimeAnCh,"Time vs Anode",240,0,240,128,0,128);
00480         }
00481       }
00482     }
00483   }
00484   mOcupancyHisto = new TH1F("OcupancyHisto","Anode Occupancy in an event",129,0,128);
00485   EventOccupancy = new TH1F("EventOccup","Event Occupancy",200,0,100000);
00486 
00487   return kStOK;
00488 }
00489 
00490 //____________________________________________________________________________
00491 
00492 Int_t StSvtSeqAdjMaker::Make()
00493 {
00494 
00495   int Anode;
00496 
00497   StSvtHybridBadAnodes* BadAnode=NULL;
00498 
00499   //  GetBadAnodes();
00500   if ( GetSvtRawData() ) return kStWarn; // Return if SVT not there
00501 
00502   SetSvtData();
00503 
00504   // copy event information to adjusted data collection
00505 //VP  (*mSvtAdjData) = (*mSvtRawData);  // already made in SetSvtData()
00506 
00507   gMessMgr->Info() << "Working On Event: " << mSvtAdjData->getEventNumber() << endm;
00508 
00509   
00510   Evt_counts=0;
00511 
00512   for(int Barrel = 1;Barrel <= mSvtRawData->getNumberOfBarrels();Barrel++){
00513     for (int Ladder = 1;Ladder <= mSvtRawData->getNumberOfLadders(Barrel);Ladder++){
00514       for (int Wafer = 1;Wafer <= mSvtRawData->getNumberOfWafers(Barrel);Wafer++){
00515         for( int Hybrid = 1;Hybrid <=mSvtRawData->getNumberOfHybrids();Hybrid++){
00516 
00517           int index = mSvtRawData->getHybridIndex(Barrel,Ladder,Wafer,Hybrid);
00518 
00519           if( index < 0) continue;
00520 
00521           mHybridRawData = (StSvtHybridData *)mSvtRawData->at(index);
00522 
00523           // cout << "mHybridRawData = " << mHybridRawData << endl;
00524           if( !mHybridRawData) continue;
00525 
00526           // subtract pedestal from raw data and save pedestal subtracted data
00527           //if( mSvtPedSub){
00528           //  mSvtPedSub->SubtractPed(mHybridRawData, index, mPedOffSet);
00529           //}
00530 
00531           // retrieve bad anodes
00532           if (mSvtBadAnodes)
00533             BadAnode = (StSvtHybridBadAnodes*)mSvtBadAnodes->at(index);
00534 
00535           mHybridAdjData = (StSvtHybridData *)mSvtAdjData->at(index);
00536           mSvtAdjData->put_at(0,index);
00537           delete mHybridAdjData;
00538           mHybridAdjData = new StSvtHybridData(Barrel,Ladder,Wafer,Hybrid); 
00539           mHybridAdjData->setTimeZero(mHybridRawData->getTimeZero());
00540           mHybridAdjData->setSCAZero(mHybridRawData->getSCAZero());
00541 
00542           doCommon = 1;
00543           mNAnodes = FindBlackAnodes();
00544 
00545           if( doCommon || !mNAnodes ){
00546             if (Debug())   gMessMgr->Debug() << "Doing Common mode average for index " << index << endl;
00547 
00548             // Zero Common Mode Noise arrays
00549             for(int Timebin=0; Timebin<128; Timebin++){
00550               mCommonModeNoise[Timebin]=0;
00551               mCommonModeNoiseAn[Timebin]=0;    
00552             } 
00553           }
00554           for( int iAnode= 0; iAnode< mHybridRawData->getAnodeList(anolist); iAnode++)
00555            {
00556              Anode = anolist[iAnode];  
00557 
00558           // here anode is real anode number (1-240)
00559           //if (Debug())MakeHistogramsAdc(mHybridRawData,index,Anode,1);
00560              if( doCommon)  CommonModeNoiseCalc(iAnode);
00561            }
00562          
00563           if( doCommon){
00564             int TimeLast, TimeAv, TimeSum, TimeAvSav;
00565             TimeLast = 0;
00566             TimeSum = 0;
00567             TimeAv=0;
00568             TimeAvSav = 0;
00569             for( int TimeBin=0; TimeBin<128; TimeBin++){
00570               if(mCommonModeNoiseAn[TimeBin] > 30)
00571                 mCommonModeNoise[TimeBin] /= mCommonModeNoiseAn[TimeBin];
00572               else  mCommonModeNoise[TimeBin] =0;
00573               
00574               if( Debug()){
00575                 if( index < 4 && mCommonModeNoiseAn[TimeBin] > 20){
00576                   if( TimeLast < TimeBin-3 && TimeSum > 0){
00577                     
00578                     TimeAv /= TimeSum;
00579                     TimeLast = TimeBin;
00580                     TimeAvSav = TimeAv;
00581                     TimeAv = 0;
00582                     TimeSum=0;
00583                   }
00584                   else{
00585                     TimeAv +=TimeBin;
00586                     TimeSum++;
00587                     TimeLast = TimeBin;
00588                   }
00589                 }
00590               } // End of if If(Debug)
00591             }
00592           } // End of If(Common)
00593         
00594           for( int iAnode= 0; iAnode<mHybridRawData->getAnodeList(anolist); iAnode++)
00595             {
00596               Anode = anolist[iAnode];
00597 
00598               // Skip Bad anodes
00599               //
00600               if( BadAnode && BadAnode->isBadAnode(Anode)) continue;
00601 
00602               if( doCommon) CommonModeNoiseSub(iAnode);
00603               else SubtractFirstAnode(iAnode);
00604               //if (Debug() && !doCommon)MakeHistogramsAdc(mHybridRawData,index,Anode,1);
00605               //Perform Asic like zero suppression
00606               AdjustSequences1(iAnode, Anode);
00607               
00608               //Perform E896 type zero-suppresion (look for non-noise like signals
00609               if(m_inv_prod_lo){
00610                 mInvProd->FindInvProducts(mHybridAdjData,iAnode,mPedOffSet);
00611                 // if(Debug())MakeHistogramsProb(index,Anode);
00612                 AdjustSequences2(iAnode, Anode);
00613 
00614               }
00615               // Here anode is index to real anode number (1-240)
00616               MakeHistogramsAdc(mHybridRawData,index,Anode,2);
00617             }
00618           
00619 
00620           mHybridAdjData->setAnodeList();
00621           mSvtAdjData->put_at(mHybridAdjData,index);
00622         }
00623 
00624         mInvProd->ResetBuffer();
00625       }
00626     }
00627   }
00628   cout << endl;
00629   cout << " Event Occupancy = " << Evt_counts << endl;
00630 
00631   EventOccupancy->Fill((float)Evt_counts);
00632   Evt_counts=0;
00633 
00634   return kStOK;
00635 
00636 }
00637 
00638 //_____________________________________________________________________________
00639 
00640 Int_t StSvtSeqAdjMaker::AdjustSequences1(int iAnode, int Anode){
00641 
00642   //Perform ASIC like zero suppression. Need > m_n_seq_lo pixels with a count
00643   // > than m_thresh_lo and  > m_n_seq_hi pixels with a count
00644   // > than m_thresh_hi  
00645 
00646   int  nSeqOrig, nSeqNow, nSeqTruth, length, count1, count2;
00647   int startTimeBin,  status;
00648   StSequence* Sequence;
00649   StMCTruth*   SeqTruth;
00650   unsigned char* adc;
00651   //int ExtraBefore = 1;
00652   // int ExtraAfter = 3;
00653   int ExtraBefore = 0;
00654   int ExtraAfter = 0;
00655   int firstTimeBin;
00656 
00657   //Anode is the index into the anolist array
00658    status= mHybridRawData->getListSequences(iAnode,nSeqOrig,Sequence);
00659    status= mHybridRawData->getListTruth     (iAnode,nSeqTruth,SeqTruth );
00660   //status= mHybridRawData->getSequences(Anode,nSeqOrig,Sequence);
00661 
00662   nSeqNow = 0;
00663   StSequence tempSeq1[128];
00664   StMCTruth  tempTru1[128];
00665     for( int nSeq=0; nSeq< nSeqOrig ; nSeq++){
00666     
00667     adc=Sequence[nSeq].firstAdc;
00668     length = Sequence[nSeq].length;
00669     startTimeBin=Sequence[nSeq].startTimeBin;
00670     
00671     int j =0;
00672     while( j<length){
00673       count1=0;
00674       count2=0;
00675       
00676       while( (int)adc[j] > m_thresh_lo && j<length){
00677         count1++;
00678 
00679         if( (int)adc[j] > m_thresh_hi){
00680           count2++;
00681         }
00682 
00683         j++;
00684       }
00685 
00686       if( count2 > m_n_seq_hi && count1 > m_n_seq_lo){
00687 
00688         firstTimeBin = startTimeBin + j - count1 - ExtraBefore;
00689 
00690         tempSeq1[nSeqNow].firstAdc=&adc[j- count1 - ExtraBefore];
00691         tempSeq1[nSeqNow].startTimeBin = firstTimeBin;
00692         tempSeq1[nSeqNow].length = 0;
00693         if((startTimeBin + j - count1 - ExtraBefore)  < 0){
00694           tempSeq1[nSeqNow].startTimeBin=0;
00695           tempSeq1[nSeqNow].firstAdc=&adc[0];
00696           // Make length temp. negative to account for the fact we went
00697           //past the start of the sequence get adjusted properly in a bit
00698           tempSeq1[nSeqNow].length = startTimeBin + j - count1 - ExtraBefore;
00699         }
00700 
00701         if((startTimeBin + j - count1 - ExtraBefore)  < startTimeBin){
00702           tempSeq1[nSeqNow].startTimeBin=startTimeBin;
00703           tempSeq1[nSeqNow].firstAdc=&adc[0];
00704           // Make length temp. negative to account for the fact we went
00705           //past the start of the sequence get adjusted properly in a bit
00706           tempSeq1[nSeqNow].length = startTimeBin + j - count1 - ExtraBefore
00707             -startTimeBin;
00708         }
00709 
00710         // Make length proper size even if it went negative
00711         tempSeq1[nSeqNow].length += count1+ ExtraAfter+ ExtraBefore;
00712 
00713         // Check dont go past end of sequence if do adjust length correctly
00714         if( tempSeq1[nSeqNow].length + tempSeq1[nSeqNow].startTimeBin  > 128) 
00715           tempSeq1[nSeqNow].length=128-tempSeq1[nSeqNow].startTimeBin +1;
00716         if ( tempSeq1[nSeqNow].startTimeBin+ tempSeq1[nSeqNow].length >
00717              startTimeBin +length){
00718           tempSeq1[nSeqNow].length =  (startTimeBin +length) -
00719             tempSeq1[nSeqNow].startTimeBin ;
00720         }
00721         if (SeqTruth) tempTru1[nSeqNow]=SeqTruth[nSeq];
00722         nSeqNow++;
00723       }
00724       
00725       
00726       j++;
00727     }
00728   }
00729 
00730   mNumOfSeq=0;
00731   if( nSeqNow > 0){
00732   if (SeqTruth) SeqTruth=tempTru1;
00733   mNumOfSeq = MergeSequences(tempSeq1,nSeqNow,SeqTruth);
00734   }
00735 
00736   mHybridAdjData->setListSequences(iAnode, Anode,mNumOfSeq, tempSeq1);
00737   if (SeqTruth)
00738   mHybridAdjData->setListTruth    (iAnode, Anode,mNumOfSeq, tempTru1);
00739   
00740   //if (nSeqNow)
00741   //  cout << "For Anode=" << Anode << " Number sequnces was=" << nSeqOrig << " Number now=" << nSeqNow << endl;
00742   
00743   //  data << endl;
00744   
00745   return kStOK;
00746 }
00747 
00748 //_____________________________________________________________________________
00749 
00750 Int_t StSvtSeqAdjMaker::AdjustSequences2(int iAnode, int Anode){
00751 
00752   //Perform E896 like zero suppression. Find pixels that have consecutive ADC 
00753   // counts that do not have the shape of noise
00754 
00755   int nSeqBefore, nSeqNow, count1, nSeqTruth;
00756   int startTimeBin, length, status;
00757   StSequence* Sequence;
00758   StMCTruth*  SeqTruth;
00759   unsigned char* adc;
00760   //int ExtraBefore=1;
00761   //int ExtraAfter=3;
00762  int ExtraBefore=0;
00763   int ExtraAfter=0;
00764 
00765   int firstTimeBin;
00766   
00767   double tempBuffer = 0;
00768   
00769   nSeqNow = 0;
00770   StSequence tempSeq1[128];
00771   StMCTruth  tempTru1[128];
00772   status = mHybridAdjData->getListSequences(iAnode,nSeqBefore,Sequence);
00773   status = mHybridRawData->getListTruth    (iAnode,nSeqTruth ,SeqTruth );
00774   
00775   for(int Seq = 0; Seq < nSeqBefore; Seq++) 
00776     {
00777       startTimeBin =Sequence[Seq].startTimeBin; 
00778       length = Sequence[Seq].length;
00779       adc = Sequence[Seq].firstAdc;
00780       
00781       int j=0;
00782       while( j < length)
00783         { 
00784           count1 = 0;
00785           tempBuffer = mInvProd->GetBuffer(startTimeBin + j);
00786           
00787           while(tempBuffer > m_inv_prod_lo && j < length)
00788             {
00789               ++count1;
00790               ++j;
00791               tempBuffer = mInvProd->GetBuffer(startTimeBin + j);
00792             }
00793           if(count1 > 0  && (tempBuffer < m_inv_prod_lo || j == length))
00794             {
00795               firstTimeBin = startTimeBin + j - count1 - ExtraBefore;
00796               
00797               tempSeq1[nSeqNow].firstAdc=&adc[j- count1 - ExtraBefore];
00798               tempSeq1[nSeqNow].startTimeBin = firstTimeBin;
00799               tempSeq1[nSeqNow].length = 0;
00800               if((startTimeBin + j - count1 - ExtraBefore)  < 0){
00801                 tempSeq1[nSeqNow].startTimeBin=0;
00802                 tempSeq1[nSeqNow].firstAdc=&adc[0];
00803                 // Make length temp. negative to account for the fact we went
00804                 //past the start of the sequence get adjusted properly in a bit
00805                 tempSeq1[nSeqNow].length = startTimeBin + j - count1 - ExtraBefore;
00806               }
00807               
00808               if((startTimeBin + j - count1 - ExtraBefore)  < startTimeBin){
00809                 tempSeq1[nSeqNow].startTimeBin=startTimeBin;
00810                 tempSeq1[nSeqNow].firstAdc=&adc[0];
00811                 // Make length temp. negative to account for the fact we went
00812                 //past the start of the sequence get adjusted properly in a bit
00813                 tempSeq1[nSeqNow].length = startTimeBin + j - 
00814                   count1 - ExtraBefore -startTimeBin;
00815               }
00816               
00817               // Make length proper size even if it went negative
00818               tempSeq1[nSeqNow].length += count1+ ExtraAfter+ ExtraBefore;
00819               
00820               // Check dont go past end of sequence if do adjust length correctly
00821               if( tempSeq1[nSeqNow].length + tempSeq1[nSeqNow].startTimeBin  > 128) 
00822                 tempSeq1[nSeqNow].length=128-tempSeq1[nSeqNow].startTimeBin +1;
00823               if ( tempSeq1[nSeqNow].startTimeBin+ tempSeq1[nSeqNow].length >
00824                    startTimeBin +length){
00825                 tempSeq1[nSeqNow].length =  (startTimeBin +length) -
00826                   tempSeq1[nSeqNow].startTimeBin ;
00827               }
00828               if(SeqTruth) tempTru1[nSeqNow] = SeqTruth[Seq];
00829               nSeqNow++;
00830             
00831             }
00832           j++;
00833         }
00834           
00835     } // Sequence loop
00836   
00837   
00838   if( nSeqBefore >0){
00839     if (SeqTruth) SeqTruth = tempTru1;
00840     mNumOfSeq = MergeSequences(tempSeq1, nSeqNow,SeqTruth);
00841     mHybridAdjData->setListSequences(iAnode, Anode,mNumOfSeq, tempSeq1);
00842     if (SeqTruth)
00843     mHybridAdjData->setListTruth    (iAnode, Anode,mNumOfSeq, tempTru1);
00844   }
00845   //cout << "For Anode=" << Anode << " Number of sequnces was=" << nSeqBefore << " Number now=" << nSeqNow << endl;
00846   
00847   return kStOK;
00848 
00849 }
00850 //_____________________________________________________________________________
00851 
00852 Int_t StSvtSeqAdjMaker::MergeSequences( StSequence* seq, int nSeq, StMCTruth* tru){
00853 
00854   // Check and see if the end of one sequence overlaps the start of the next
00855   // if it does merge sequences
00856 
00857   int nSeqNow = 0;
00858   int EndTime;
00859   StMCPivotTruth pivo(1);
00860   if (tru) pivo.Add(tru[0]);
00861   for( int i=1; i<nSeq; i++){
00862     
00863     if( (seq[nSeqNow].startTimeBin + seq[nSeqNow].length) 
00864         >= seq[i].startTimeBin){
00865       EndTime = seq[i].startTimeBin + seq[i].length;
00866       seq[nSeqNow].length = EndTime - seq[nSeqNow].startTimeBin;
00867       if (tru) { pivo.Add(tru[i]); tru[nSeqNow] = pivo.Get();}
00868       
00869     }
00870     else{
00871       nSeqNow++;
00872       seq[nSeqNow] = seq[i];
00873       if (tru){ pivo.Reset(); pivo.Add(tru[i]);} 
00874     }
00875 
00876   }
00877 
00878   return nSeqNow+1;
00879 }
00880 //_____________________________________________________________________________
00881 
00882 void StSvtSeqAdjMaker::CommonModeNoiseCalc(int iAnode){
00883 
00884   // Calc common mode noise
00885 
00886   int  nSeqOrig, length;
00887   int startTimeBin,  status;
00888   StSequence* Sequence;
00889   unsigned char* adc;
00890 
00891   //Anode is the index into the anolist array
00892   
00893   status= mHybridRawData->getListSequences(iAnode,nSeqOrig,Sequence);
00894 
00895   for( int nSeq=0; nSeq< nSeqOrig ; nSeq++){
00896   
00897     adc=Sequence[nSeq].firstAdc;
00898     length = Sequence[nSeq].length;
00899     startTimeBin=Sequence[nSeq].startTimeBin;
00900 
00901     int j =0;
00902     while( j<length){
00903       
00904       mCommonModeNoise[j+startTimeBin] += (int)adc[j] - mPedOffSet;
00905       mCommonModeNoiseAn[j+startTimeBin]++;
00906       j++;
00907  
00908 
00909     }
00910   }
00911   return;
00912 }
00913 //_____________________________________________________________________________
00914 
00915 void StSvtSeqAdjMaker::CommonModeNoiseSub(int iAnode){
00916 
00917   // Calc common mode noise
00918 
00919   int  nSeqOrig, length;
00920   int startTimeBin,  status;
00921   StSequence* Sequence;
00922   unsigned char* adc;
00923 
00924   //Anode is the index into the anolist array
00925   
00926   status= mHybridRawData->getListSequences(iAnode,nSeqOrig,Sequence);
00927 
00928   for( int nSeq=0; nSeq< nSeqOrig ; nSeq++){
00929   
00930     adc=Sequence[nSeq].firstAdc;
00931     length = Sequence[nSeq].length;
00932     startTimeBin=Sequence[nSeq].startTimeBin;
00933 
00934     int j =0;
00935     while( j<length){
00936       if( (int) adc[j]-mCommonModeNoise[j+startTimeBin] < 0) adc[j]=0;
00937       else{
00938         adc[j] -=mCommonModeNoise[j+startTimeBin];
00939       }
00940       j++;
00941     }
00942   }
00943   return;
00944 }
00945 
00946 //_____________________________________________________________________________
00947 
00948 void StSvtSeqAdjMaker::SubtractFirstAnode(int iAnode){
00949 
00950   // Calc common mode noise
00951 
00952   int  nSeq, nSeqOrig, length;
00953   int startTimeBin,  status, j;
00954   int adcMean;
00955   StSequence* Sequence;
00956   unsigned char* adc;
00957 
00958 
00959   status= mHybridRawData->getListSequences(iAnode,nSeqOrig,Sequence);
00960 
00961   adcMean = 0;
00962   for( nSeq=0; nSeq< nSeqOrig ; nSeq++){
00963   
00964     adc=Sequence[nSeq].firstAdc;
00965     length = Sequence[nSeq].length;
00966     startTimeBin=Sequence[nSeq].startTimeBin; 
00967     j =0;
00968     while( j<length){
00969 
00970      if( mNAnodes) adcMean = adcCommon[j+startTimeBin]/mNAnodes;
00971 
00972      //cout << "iAnode = " << iAnode << "ADC value is: " << (int) adc[j] << " for time = " << startTimeBin + j << " and ADCMean = "<< adcMean << " mNAnodes = " << mNAnodes <<endl;
00973       if( (int)adc[j]- adcMean < 0) adc[j]=0;
00974       else if( adcMean < 0){
00975         adc[j]  += (unsigned char)(-1*adcMean);
00976       }
00977       else{
00978         adc[j] -= (unsigned char)adcMean;
00979       }
00980       //cout << "Final ADC= " << (int)adc[j] << endl;
00981       j++;
00982     }
00983   }
00984   return;
00985 }
00986 //____________________________________________________________________________
00987 
00988 int StSvtSeqAdjMaker::FindBlackAnodes(){
00989   // Decide whether to do common mode noise from average timebucket 
00990   //value or via first two black anodes.Have many options andoes 1&240,
00991   //1&2 or 2&239. IN some case one of the black anodes is bad so only 
00992   //subtract one.
00993   int status, length;
00994   int i, j, startTimeBin;
00995   int nSequence = 0;
00996   StSequence* Seq = NULL;
00997   unsigned char *adc;
00998   
00999   doCommon = 0;
01000   mNAnodes=2;
01001   int adcAv=0;
01002   
01003   status= mHybridRawData->getSequences(240,nSequence,Seq);
01004   length = 0;
01005   for(j=0; j<128; j++) adcCommon[j]=0;
01006   for( i=0; i<nSequence; i++)  length += Seq[i].length;
01007   if( length < 126){
01008     // IF anode 240 isnt black try 239
01009     status= mHybridRawData->getSequences(239,nSequence,Seq);
01010     length = 0;
01011     
01012     for( i=0; i<nSequence; i++)  length += Seq[i].length;
01013     if( length < 126){
01014       //if that one isnt black either try 2
01015       status= mHybridRawData->getSequences(2,nSequence,Seq);
01016       length = 0;
01017       
01018       for( i=0; i<nSequence; i++)  length += Seq[i].length;
01019       if( length < 126){
01020         // If that isnt black too do common mode noise subtraction
01021         doCommon =1;
01022       }
01023     }
01024   }
01025   
01026   // If one of those anodes was black fill adc array
01027   if( !doCommon){
01028     
01029     // Fill in array of adc values of first anode
01030     for( int nSeq=0; nSeq< nSequence ; nSeq++){
01031       
01032       adc=Seq[nSeq].firstAdc;
01033       length = Seq[nSeq].length;
01034       startTimeBin=Seq[nSeq].startTimeBin;
01035       j=0;
01036       while( j<length){
01037         adcCommon[startTimeBin+j] = adc[j]-mPedOffSet;
01038         adcAv += (int) adc[j];
01039         j++;
01040       }
01041     }
01042     // If this black anode was bad dont count it in subtraction
01043     if( adcAv <100 || adcAv> 25000) {
01044       mNAnodes--;
01045       for(j=0; j<128; j++) adcCommon[j]=0;
01046     }
01047   }
01048   
01049   // now find second black anode
01050   
01051   status= mHybridRawData->getSequences(2,nSequence,Seq);
01052   length = 0;
01053   
01054   for( i=0; i<nSequence; i++)  length += Seq[i].length;
01055   if( length < 126){
01056     status= mHybridRawData->getSequences(1,nSequence,Seq);
01057     length = 0;   
01058     
01059     for( i=0; i<nSequence; i++)  length += Seq[i].length;
01060     if( length < 126){
01061       doCommon=1;
01062     }
01063     
01064   }
01065   if( !doCommon){
01066     // Fill in array of adc values of second anode
01067     adcAv=0;
01068     for( int nSeq=0; nSeq< nSequence ; nSeq++){
01069       
01070       adc=Seq[nSeq].firstAdc;
01071       length = Seq[nSeq].length;   
01072       startTimeBin=Seq[nSeq].startTimeBin;
01073       j=0;
01074       while( j<length){
01075         adcCommon[startTimeBin+j] += adc[j]-mPedOffSet;
01076         adcAv += (int)  adc[j];
01077         j++;
01078       }
01079     }
01080 
01081     if( adcAv <100 || adcAv> 25000) {
01082       mNAnodes--;
01083       // If black anode is bad subtract back of wrong adc values
01084       // from average will subtract from all other anodes
01085       if( mNAnodes){
01086         for( int nSeq=0; nSeq< nSequence ; nSeq++){
01087           
01088           adc=Seq[nSeq].firstAdc;
01089           length = Seq[nSeq].length;
01090           startTimeBin=Seq[nSeq].startTimeBin;
01091           j=0;
01092           while( j<length){
01093             adcCommon[startTimeBin+j] -= adc[j]-mPedOffSet;
01094             adcAv += (int)  adc[j];
01095             j++;
01096           }
01097         }
01098       }
01099     }
01100   }
01101   return mNAnodes;
01102 }
01103 
01104 //_____________________________________________________________________________
01105 void StSvtSeqAdjMaker::MakeHistogramsAdc(StSvtHybridData* hybridData, int index, int Anode, int Count){
01106   
01107   // Count is 1 for before CM-Noise subtraction and 2 for after CMN subtraction.
01108   int mSequence;
01109   int len,status;
01110   unsigned char* adc;
01111 
01112   StSequence* svtSequence;
01113   mSequence = 0;
01114   svtSequence = NULL;
01115 
01116   // Anode goes between 1-240
01117   
01118   status = hybridData->getSequences(Anode,mSequence,svtSequence);
01119   
01120   int numOfPixels=0;
01121   for(int mSeq = 0; mSeq < mSequence; mSeq++) 
01122     { 
01123       adc = svtSequence[mSeq].firstAdc;
01124       len = svtSequence[mSeq].length;
01125       if (len>12) continue;
01126       for(int j = 0 ; j < len; j++){
01127         mRawAdc[index]->Fill(Anode,(int)adc[j]);
01128         //mTimeAn[index]->Fill(Anode,j+svtSequence[mSeq].startTimeBin,(int)adc[j]);
01129         mAdcAfter[index]->Fill(Anode);
01130         //cout << "Index= " << index << " Anode " << Anode << " " << (int)adc[j] - mPedOffSet << " count= " << Evt_counts <<endl;
01131         if (adc[j] >0) numOfPixels++;
01132         if (adc[j] >0) Evt_counts++;
01133       }
01134     }
01135   mOcupancyHisto->Fill(numOfPixels);
01136 }
01137 
01138 //_____________________________________________________________________________
01139 Int_t StSvtSeqAdjMaker::Reset(){
01140 
01141   //delete mInvProd;
01142   //delete mProbValue;
01143 
01144   mSvtDataSet = NULL;
01145   mSvtAdjData = NULL;
01146   mSvtRawData = NULL;
01147   mHybridRawData = NULL;
01148   mHybridAdjData = NULL;
01149   mSvtPedSub = NULL;
01150   mSvtPedColl = NULL;
01151   mSvtBadAnodes = NULL;
01152   mInvProd = NULL;
01153   mProbValue = NULL;
01154 
01155   m_ConstSet->Delete();
01156   
01157   return kStOK;
01158 }
01159 //_____________________________________________________________________________
01160 string StSvtSeqAdjMaker::baseName(string s){
01161 
01162   string name(s);
01163   size_t pos;
01164   pos = name.find_last_of("/");
01165   if (pos!=string::npos ) name.erase(0, pos );
01166   pos = name.find_first_of(".");
01167   if (pos!=string::npos ) name.erase(pos,name.length()-pos );
01168   return name;
01169 } 
01170 //_____________________________________________________________________________
01171 string StSvtSeqAdjMaker::buildFileName(string dir, string fileName, string extention){
01172 
01173   fileName = dir + fileName + extention;
01174   while (fileName.find("//")!=string::npos) {
01175     int pos = fileName.find("//");
01176     fileName.erase(pos,1);
01177   }
01178   return fileName;
01179 }
01180 //_____________________________________________________________________________
01181 Int_t StSvtSeqAdjMaker::Finish(){
01182   
01183   
01184   if (Debug()) gMessMgr->Debug() << "In StSvtSeqAdjMaker::Finish() "
01185                                  << GetName() << endm; 
01186 
01187   if ( hfile ){
01188     hfile->Write();
01189     hfile->Close();
01190   }
01191   
01192   return kStOK;
01193 }
01194 
01195 //_____________________________________________________________________________
01196 ClassImp(StSvtSeqAdjMaker)