StRoot: StSvtCalibMaker/StSvtDriftVelocityMaker.cxx Source File

00001 /***************************************************************************
00002  *
00003  * $Id: StSvtDriftVelocityMaker.cxx,v 1.6 2004/01/26 23:10:49 perev Exp $
00004  *
00005  * Author: Marcelo Munhoz
00006  ***************************************************************************
00007  *
00008  * Description: SVT Drift Velocity calculation Maker
00009  *
00010  ***************************************************************************
00011  *
00012  * $Log: StSvtDriftVelocityMaker.cxx,v $
00013  * Revision 1.6  2004/01/26 23:10:49  perev
00014  * Leak off
00015  *
00016  * Revision 1.5  2003/09/02 17:59:05  perev
00017  * gcc 3.2 updates + WarnOff
00018  *
00019  * Revision 1.4  2002/03/04 17:06:48  willson
00020  * Laser spot positions recorded
00021  *
00022  * Revision 1.3  2002/02/06 00:10:45  willson
00023  * File entries, variable names.
00024  *
00025  * Revision 1.2  2002/01/18 20:53:47  willson
00026  * Some bug fixes by Helen
00027  *
00028  * Revision 1.1  2002/01/18 19:57:42  willson
00029  * Drift Velocity Calculation v.1
00030  *
00031  *
00032  **************************************************************************/
00033 
00034 #include <assert.h>
00035 #include "StSvtDriftVelocityMaker.h"
00036 
00037 #include "StMessMgr.h"
00038 
00039 // #include "StSvtDb/StSvtDbMaker.h"
00040 
00041 
00042 #include "StEventTypes.h"
00043 #include "StEvent/StPrimaryVertex.h"
00044 #include "StEvent/StEvent.h"
00045 #include "StarClassLibrary/StSequence.hh"
00046 #include "StSvtClassLibrary/StSvtHybridData.hh"
00047 #include "StSvtClassLibrary/StSvtData.hh"
00048 #include "StSvtClassLibrary/StSvtConfig.hh"
00049 #include "StSvtClassLibrary/StSvtHybridPed.hh"
00050 #include "StSvtClassLibrary/StSvtHybridDriftVelocity.hh"
00051 #include "StSvtClassLibrary/StSvtHybridCollection.hh"
00052 #include "StSvtClassLibrary/StSvtHybridPixels.hh"
00053 #include "StSvtClusterMaker/StSvtAnalysedHybridClusters.hh"
00054 #include "Stiostream.h"
00055 #include <fstream>
00056 
00057 #define N_INJECTOR_LINES 4
00058 
00059 ClassImp(StSvtDriftVelocityMaker)
00060 
00061 //_____________________________________________________________________________
00062   StSvtDriftVelocityMaker::StSvtDriftVelocityMaker(const char *name):StMaker(name)
00063 {
00064   int n = (char*)&mDebug - (char*)&mSvtData + sizeof(mDebug);
00065   memset(&mSvtData,0,n);
00066 
00067   mNumTimeBins = 128;
00068   mMaximumTB = 128;
00069   mMinimumTB = 0;
00070   mDebug = false;
00071   mFraction = 0.5;
00072   mMoveForward = true;
00073   mT0Guess = 0.5;
00074 }
00075 
00076 //_____________________________________________________________________________
00077 StSvtDriftVelocityMaker::~StSvtDriftVelocityMaker()
00078 {
00079   delete mSvtDriftVeloc; mSvtDriftVeloc=0;
00080   for (int i=0;i<mNHybridDriftVelocityHisto;i++) {
00081     delete mHybridDriftVelocityHisto[i]; 
00082     delete mHybridDriftVelocity2DHisto[i];}
00083   delete [] mHybridDriftVelocityHisto;
00084   delete [] mHybridDriftVelocity2DHisto;
00085 
00086   delete mGlobalDriftVelocityHisto; 
00087   delete mGlobalDriftVelocity2DHisto; 
00088   delete mCalculatedDriftVelocity; 
00089   delete mLaserSpotDistL07B3_1; 
00090   delete mLaserSpotDistL15B3_1; 
00091   delete mLaserSpotDistL15B3_2; 
00092 
00093 }
00094 
00095 //_____________________________________________________________________________
00096 Int_t StSvtDriftVelocityMaker::Init()
00097 {
00098   //cout << "StSvtDriftVelocityMaker::Init" << endl;
00099 
00100   SetSvtRawData();
00101 
00102   SetSvtDriftVelocity();
00103 
00104   int indexHybrid;
00105   StSvtHybridDriftVelocity* hybridDriftVeloc;
00106   TH1D* hybridHisto;
00107   TH2D* hybrid2DHisto;
00108   char CharString1[100];
00109   char CharString2[100];
00110   mNHybridDriftVelocityHisto  = mSvtRawData->getTotalNumberOfHybrids();
00111   mHybridDriftVelocityHisto = new TH1D*[mNHybridDriftVelocityHisto];
00112   
00113   if (mDebug)
00114     mHybridDriftVelocity2DHisto = new TH2D*[mNHybridDriftVelocityHisto];
00115 
00116   // Loop over barrels, ladders, wafers and hybrids
00117   for (int barrel = 1;barrel <= mSvtRawData->getNumberOfBarrels();barrel++) {
00118     for (int ladder = 1;ladder <= mSvtRawData->getNumberOfLadders(barrel);ladder++) {
00119       for (int wafer = 1;wafer <= mSvtRawData->getNumberOfWafers(barrel);wafer++) {
00120         for (int hybrid = 1;hybrid <= mSvtRawData->getNumberOfHybrids();hybrid++) {
00121 
00122           indexHybrid = mSvtRawData->getHybridIndex(barrel, ladder, wafer, hybrid);
00123 
00124           if (indexHybrid < 0) continue;
00125 
00126           hybridDriftVeloc = (StSvtHybridDriftVelocity*)mSvtRawData->at(indexHybrid);
00127           if (!hybridDriftVeloc)
00128             hybridDriftVeloc = new StSvtHybridDriftVelocity(barrel, ladder, wafer, hybrid);
00129             
00130           mSvtRawData->put_at(hybridDriftVeloc,indexHybrid);
00131 
00132           sprintf(CharString1, "Histo%d", indexHybrid);
00133           sprintf(CharString2, "2DHisto%d", indexHybrid);
00134 
00135           hybridHisto = new TH1D(CharString1, CharString1, mNumTimeBins, mMinimumTB, mMaximumTB);
00136           mHybridDriftVelocityHisto[indexHybrid] = hybridHisto;
00137 
00138           if (mDebug) {
00139             hybrid2DHisto = new TH2D(CharString1, CharString1, 128, 0, 128, 240, 0, 240);
00140             mHybridDriftVelocity2DHisto[indexHybrid] = hybrid2DHisto;
00141           }
00142 
00143         }  // end of loop over hybrids
00144       }  // end of loop over wafers
00145     }  // end of loop over ladders
00146   }  // end of loop over barrels
00147   
00148   if (mDebug) {
00149     mGlobalDriftVelocityHisto = new TH1D("Global1", "Global1", 640, 0, 128); 
00150     mGlobalDriftVelocity2DHisto = new TH2D("Global2", "Global2", 640, 0, 128, 1200, 0, 240); 
00151   }
00152     
00153   mCalculatedDriftVelocity = new TH1D("Final", "Final", 432, 0, 431);
00154   mLaserSpotDistL07B3_1 = new TH1D("L07B3 Laser Distribution 1", "L07B3 Laser Distribution 1", 1280, 0, 128); 
00155   mLaserSpotDistL15B3_1 = new TH1D("L15B3 Laser Distribution 1", "L15B3 Laser Distribution 1", 1280, 0, 128); 
00156   mLaserSpotDistL15B3_2 = new TH1D("L15B3 Laser Distribution 2", "L15B3 Laser Distribution 2", 1280, 0, 128); 
00157 
00158   return StMaker::Init();
00159 }
00160 //______________________________________________
00161 Int_t StSvtDriftVelocityMaker::SetSvtRawData()
00162 {
00163   St_DataSet *dataSet;
00164   
00165   dataSet = GetDataSet("StSvtRawData");
00166   if( !dataSet) {
00167     gMessMgr->Warning() << " No Svt Raw data set" << endm;
00168     dataSet = new St_ObjectSet("DriftVelocity");
00169     mSvtRawData = new StSvtData("FULL");
00170     return kStWarn;
00171   }
00172     
00173 
00174   mSvtRawData = (StSvtData*)(dataSet->GetObject());
00175   if( !mSvtRawData) {
00176     gMessMgr->Warning() << " No Svt Raw data " << endm;
00177     return kStWarn;
00178   }
00179   
00180   return kStOK;
00181 }
00182 //_____________________________________________________________________________
00183 Int_t StSvtDriftVelocityMaker::SetSvtData()
00184 {
00185   St_DataSet *dataSet;
00186 
00187   dataSet = GetDataSet("StSvtAnalResults");
00188   if (dataSet) mRaw = true; else
00189     mRaw = false;
00190 
00191   if (mRaw) {
00192     mSvtData = (StSvtHybridCollection*)(dataSet->GetObject());
00193   }
00194 
00195   return kStOK;
00196 }
00197 
00198 //_____________________________________________________________________________
00199 Int_t StSvtDriftVelocityMaker::SetSvtDriftVelocity()
00200 {
00201   assert(!mSvtDriftVeloc);
00202   if (mSvtRawData)
00203     mSvtDriftVeloc = new StSvtHybridCollection(mSvtRawData->getConfiguration());
00204   else
00205     mSvtDriftVeloc = new StSvtHybridCollection("FULL");
00206 
00207     return kStOK;  
00208 }
00209 
00210 //_____________________________________________________________________________
00211 Int_t StSvtDriftVelocityMaker::Make()
00212 {
00213   cout << "StSvtDriftVelocityMaker::Make" << endl;
00214   
00215   SetSvtData();
00216 
00217   if (mRaw)
00218     FillHistogramsRaw();
00219   else
00220     FillHistogramsStEvent();
00221 
00222   cout << "Current Event  Number is " << ((StSvtData*)mSvtRawData)->getEventNumber() << endl;
00223   cout << "Current Run Number is " << ((StSvtData*)mSvtRawData)->getRunNumber() << endl;
00224 
00225 
00226   //WriteToDb("2000-02-03 07:31:00");
00227 
00228   return kStOK;
00229 }
00230 
00231 //_____________________________________________________________________________
00232 Int_t StSvtDriftVelocityMaker::FillHistogramsRaw()
00233 {
00234 
00235   //cout << "StSvtDriftVelocityMaker::FillHistogramsRaw" << endl;
00236   StSvtAnalysedHybridClusters* hybridData;
00237   //StSvtHybridDriftVelocity* hybridDriftVeloc;
00238 
00239   int indexHybrid, index, nHits;
00240 
00241   mEventCounter++;  // increase counter
00242 
00243   // Loop over barrels, ladders, wafers and hybrids
00244   for (int barrel = 1;barrel <= mSvtData->getNumberOfBarrels();barrel++) {
00245     for (int ladder = 1;ladder <= mSvtData->getNumberOfLadders(barrel);ladder++) {
00246       for (int wafer = 1;wafer <= mSvtData->getNumberOfWafers(barrel);wafer++) {
00247         for (int hybrid = 1;hybrid <= mSvtData->getNumberOfHybrids();hybrid++) {
00248 
00249           indexHybrid = mSvtData->getHybridIndex(barrel, ladder, wafer, hybrid);
00250           if (indexHybrid < 0) continue;
00251 
00252           // Get hybrid data
00253           hybridData = (StSvtAnalysedHybridClusters*)(mSvtData->at(indexHybrid));
00254 
00255           if (!hybridData) continue;
00256 
00257           nHits = hybridData->numOfHits();
00258           mHitCounter += nHits;
00259 
00260           for (index=0; index<nHits; index++) {
00261             
00262             if (accept(hybridData->svtHit())) {
00263               mHybridDriftVelocityHisto[indexHybrid]->Fill(hybridData->WaferPosition()[index].x());
00264               
00265               if (indexHybrid==293) {
00266                 if (hybridData->WaferPosition()[index].y()>196.5 && hybridData->WaferPosition()[index].y()<200.5 && hybridData->WaferPosition()[index].x()>60)
00267                   mLaserSpotDistL07B3_1->Fill(hybridData->WaferPosition()[index].x());
00268               }
00269 
00270               if (indexHybrid==416) {
00271                 if (hybridData->WaferPosition()[index].y()>198.5 && hybridData->WaferPosition()[index].y()<200.5 && hybridData->WaferPosition()[index].x()>60)
00272                   mLaserSpotDistL15B3_1->Fill(hybridData->WaferPosition()[index].x());
00273                 if (hybridData->WaferPosition()[index].y()>196.5 && hybridData->WaferPosition()[index].y()<198.5 && hybridData->WaferPosition()[index].x()>60)
00274                   mLaserSpotDistL15B3_2->Fill(hybridData->WaferPosition()[index].x());
00275               }
00276 
00277               if (mDebug) {
00278                 mHybridDriftVelocity2DHisto[indexHybrid]->Fill(hybridData->WaferPosition()[index].x(), hybridData->WaferPosition()[index].y());
00279                 mGlobalDriftVelocityHisto->Fill(hybridData->WaferPosition()[index].x());
00280                 mGlobalDriftVelocity2DHisto->Fill(hybridData->WaferPosition()[index].x(), hybridData->WaferPosition()[index].y());
00281               }
00282 
00283             }
00284           }       
00285         }  // end of loop over hybrids
00286       }  // end of loop over wafers
00287     }  // end of loop over ladders
00288   }  // end of loop over barrels
00289       
00290   return kStOK;
00291 }
00292 
00293 
00294 //_____________________________________________________________________________
00295 Int_t StSvtDriftVelocityMaker::FillHistogramsStEvent()
00296 {
00297   //cout << "StSvtDriftVelocityMaker::FillHistogramsStEvent" << endl;
00298   //StSvtHybridData* hybridData;
00299   //StSvtHybridDriftVelocity* hybridDriftVeloc;
00300   
00301   //
00302   //    Get pointer to StEvent
00303   //
00304   StEvent* event;
00305   event = (StEvent *) GetInputDS("StEvent");
00306 
00307   cout << "StEvent = " << event << endl;
00308 
00309   if (!event) return kStOK;        // if no event, we're done
00310   
00311   StSvtConfig* fSvtConfig = new StSvtConfig();
00312   fSvtConfig->setConfiguration("FULL");
00313   //
00314   //  See if this event survives the event filter.
00315   //  If not we stop here right away.
00316   //
00317   if (!accept(event)) return kStOK;  
00318   
00319   mEventCounter++;  // increase counter
00320  
00321   // Get pointers to svt hit collection 
00322   
00323   StSvtHitCollection* rSvtHitColl = event->svtHitCollection();
00324   
00325   cout << "Number of hits = " << rSvtHitColl->numberOfHits() << "." << endl;
00326   mHitCounter += rSvtHitColl->numberOfHits();
00327 
00328   StSvtHit* sCurrentHit;
00329 
00330   for (unsigned int br=0; br<rSvtHitColl->numberOfBarrels(); br++) {
00331     for (unsigned int ld=0; ld<rSvtHitColl->barrel(br)->numberOfLadders(); ld++) {
00332       for (unsigned int w=0; w<rSvtHitColl->barrel(br)->ladder(ld)->numberOfWafers(); w++) {
00333         StSPtrVecSvtHit& hits = rSvtHitColl->barrel(br)->ladder(ld)->wafer(w)->hits();
00334         for (StSvtHitIterator i = hits.begin(); i != hits.end(); i++) {
00335           sCurrentHit = (*i);
00336           if( accept(sCurrentHit)) {
00337             mHybridDriftVelocityHisto[sCurrentHit->index()]->Fill(sCurrentHit->timebucket());
00338             if (mDebug)
00339               mGlobalDriftVelocityHisto->Fill(sCurrentHit->timebucket());
00340           }
00341         }
00342       }
00343     } 
00344   }
00345       
00346   return kStOK;
00347 
00348 }
00349 
00350 
00351 //_____________________________________________________________________________
00352 Int_t StSvtDriftVelocityMaker::CalcDriftVelocity()
00353 {
00354 
00355   double Aver=0, Num=0;
00356   char filename[100];
00357   Int_t Offset= (mMaximumTB - mMinimumTB)/10;
00358   Int_t Start = Offset + mMinimumTB;
00359   Int_t End   = mMaximumTB - Offset;
00360   Int_t i, j;
00361   double DriftVelocity, TotalDriftTime;
00362   ofstream file;
00363 
00364   sprintf(filename, "Run%d_DV.out", mSvtRawData->getRunNumber());
00365 
00366   file.open(filename);
00367 
00368   for (i=0; i<mSvtDriftVeloc->getTotalNumberOfHybrids(); i++) {
00369     
00370     for (j=Start; j<End; j++) {
00371       Aver += mHybridDriftVelocityHisto[i]->GetBinContent(j);
00372       Num++;
00373     }
00374     
00375     Aver /= Num;
00376 
00377     Aver *= mFraction;
00378     
00379     if (mMoveForward) {
00380       for (j=mMaximumTB/2; j<mMaximumTB+1; j++) {
00381         if (mHybridDriftVelocityHisto[i]->GetBinContent(j) < Aver) 
00382           break;
00383       }
00384     }
00385 
00386     else {
00387       for (j=mMaximumTB-2; j>mMaximumTB/2; j--) {
00388         if (mHybridDriftVelocityHisto[i]->GetBinContent(j) > Aver) 
00389           break;
00390       }
00391     }
00392     
00393       
00394 
00395     if (mDebug)
00396       cout << j << endl;
00397 
00398     TotalDriftTime = (j-mMinimumTB-mT0Guess)*(128.0/mNumTimeBins)*0.00000004;
00399 
00400     DriftVelocity = 3.0/TotalDriftTime;
00401 
00402     if (mDebug)
00403       mCalculatedDriftVelocity->Fill(i, DriftVelocity);
00404   
00405     file << i << "\t" << DriftVelocity << "\t\t" << j*128.0/mNumTimeBins << "\n";
00406     
00407     Aver=Num=0;
00408 
00409   }
00410 
00411   file << "Laser Spot Time Bucket Positions:\n";
00412 
00413   file << "L07B3_1" << "\t" << mLaserSpotDistL07B3_1->GetMean() << "\n";
00414   file << "L15B3_1" << "\t" << mLaserSpotDistL15B3_1->GetMean() << "\n";
00415   file << "L15B3_2" << "\t" << mLaserSpotDistL15B3_2->GetMean() << "\n";
00416 
00417   file.close();
00418 
00419   return kStOK;
00420 }
00421 
00422 //_____________________________________________________________________________
00423 Int_t StSvtDriftVelocityMaker::GetInjectorLine(float peak)
00424 {
00425   return 0;
00426 }
00427 
00428 //_____________________________________________________________________________
00429 Int_t  StSvtDriftVelocityMaker::GetInjectorLine(float* peak)
00430 {
00431   return 0;
00432 }
00433 
00434 //_____________________________________________________________________________
00435 Float_t  StSvtDriftVelocityMaker::GetClosestToLine(float peak1, float peak2)
00436 {
00437   return 0;
00438 }
00439 
00440 //_____________________________________________________________________________
00441 Float_t  StSvtDriftVelocityMaker::GetTimeZero(int anode)
00442 {
00443   return 0;
00444 }
00445 
00446 //_____________________________________________________________________________
00447 Float_t  StSvtDriftVelocityMaker::GetDistanceInjectorLine(int line)
00448 {
00449   return 0;
00450 }
00451 
00452 //_____________________________________________________________________________
00453 Float_t  StSvtDriftVelocityMaker::FitVelocity(int nInjectorsFired, float* peak, float t0)
00454 {
00455   return 0;
00456 }
00457 
00458 //_____________________________________________________________________________
00459 Float_t  StSvtDriftVelocityMaker::CalcV1(int anode, float velocity)
00460 {
00461   return 0;
00462 }
00463 
00464 //_____________________________________________________________________________
00465 Float_t  StSvtDriftVelocityMaker::CalcV2(int anode, float velocity)
00466 {
00467   return 0;
00468 }
00469 
00470 //_____________________________________________________________________________
00471 Float_t  StSvtDriftVelocityMaker::CalcV3(int anode, float velocity)
00472 {
00473   return 0;
00474 }
00475 
00476 //_____________________________________________________________________________
00477 Int_t  StSvtDriftVelocityMaker::FillAllAnodes()
00478 {
00479   return 0;
00480 }
00481 
00482 //_____________________________________________________________________________
00483 Int_t  StSvtDriftVelocityMaker::WriteToDb(Text_t *timestamp)
00484 {
00485   cout << "StSvtDriftVelocityMaker::WriteToDb" << endl;
00486 
00487   if (!GetMaker("SvtDb")) {
00488     gMessMgr->Warning() << " SVT database was not instantiated. Impossible to write! " << endm;
00489     return kStWarn;
00490   }
00491 
00492   //if (timestamp)
00493   //  ((StSvtDbMaker*)GetMaker("SvtDb"))->setTimeStamp(timestamp);
00494 
00495   //((StSvtDbMaker*)GetMaker("SvtDb"))->writeSvtDriftVelocity(mSvtDriftVeloc);
00496 
00497   return kStOK;
00498 }
00499 
00500 //_____________________________________________________________________________
00501 Int_t StSvtDriftVelocityMaker::Finish()
00502 {
00503   //cout << "StSvtDriftVelocityMaker::Finish" << endl;
00504 
00505   CalcDriftVelocity();
00506 
00507   Int_t indexHybrid;
00508 
00509   if (!mDebug) {
00510     // Loop over barrels, ladders, wafers and hybrids
00511     for (int barrel = 1;barrel <= mSvtDriftVeloc->getNumberOfBarrels();barrel++) {
00512       for (int ladder = 1;ladder <= mSvtDriftVeloc->getNumberOfLadders(barrel);ladder++) {
00513         for (int wafer = 1;wafer <= mSvtDriftVeloc->getNumberOfWafers(barrel);wafer++) {
00514           for (int hybrid = 1;hybrid <= mSvtDriftVeloc->getNumberOfHybrids();hybrid++) {
00515             
00516             indexHybrid = mSvtDriftVeloc->getHybridIndex(barrel, ladder, wafer, hybrid);
00517             
00518             if (indexHybrid < 0) continue;
00519             
00520             delete mHybridDriftVelocityHisto[indexHybrid];
00521             
00522           }  // end of loop over hybrids
00523         }  // end of loop over wafers
00524       }  // end of loop over ladders
00525     }  // end of loop over barrels
00526     
00527     delete [] mHybridDriftVelocityHisto;
00528   }
00529   
00530   return kStOK;
00531 }
00532 
00533 bool StSvtDriftVelocityMaker::accept(StEvent* event)
00534 {
00535     //
00536     //  This is a kind of very simple event filter.
00537     //  We select only events with a valid event vertex,
00538     //  i.e. event->primaryVertex() returns a non-zero pointer.
00539     // and  with the primary vertex under the svt
00540 
00541   if( !(event->primaryVertex())) return 0;
00542   //cout << event->primaryVertex()->position().z() << endl;GetGlobalHisto()
00543   //if( fabs(event->primaryVertex()->position().z()) < 20) return 1;
00544   
00545   return 1;
00546 }
00547 
00548 bool StSvtDriftVelocityMaker::accept(StSvtHit* hit)
00549 {
00550   //Accept hits with flag < 4 and peak ADC>15
00551 
00552 
00553   if( hit->flag() < 100) return 1;
00554   return 0;
00555 }