StRoot: StGenericVertexMaker/StGenericVertexMaker.cxx Source File

00001 
00002 // Revision 1.1.1.1  2001/01/31 14:00:07  balewski
00003 // First release
00004 //
00005 //
00007 //                                                                      //
00008 //   Maker to run minuit based vertex finder                            //
00009 //
00010 //                                                                      //
00012 
00013 #include <strings.h>
00014 #include <math.h>
00015 
00016 #include "StGenericVertexMaker.h"
00017 #include "St_DataSetIter.h"
00018 #include "StEventTypes.h"
00019 #include "TH2.h"
00020 #include "TNtuple.h"
00021 #include "StMessMgr.h"
00022 
00023 #include "StGenericVertexFinder.h"
00024 #include "StppLMVVertexFinder.h"
00025 #include "StFixedVertexFinder.h"
00026 
00027 #include "StTreeMaker/StTreeMaker.h"
00028 
00029 // Vertex finder implemtations
00030 #include "Minuit/StMinuitVertexFinder.h"
00031 #include "StiPPVertex/StPPVertexFinder.h"
00032 
00033 
00034 #include "tables/St_g2t_vertex_Table.h" // tmp for Dz(vertex)
00035 #include "tables/St_vertexSeed_Table.h" //
00036 
00037 // for Helix model
00038 #include "StarCallf77.h"
00039 extern "C" {void type_of_call F77_NAME(gufld,GUFLD)(float *x, float *b);}
00040 #define gufld F77_NAME(gufld,GUFLD)
00041 #include "SystemOfUnits.h"
00042 #ifndef ST_NO_NAMESPACES
00043 using namespace units;
00044 #endif
00045 
00046 
00047 ClassImp(StGenericVertexMaker)
00048 //___________________________________________________________
00049 StGenericVertexMaker::StGenericVertexMaker(const char *name):StMaker(name)
00050 {
00051   useITTF=kTRUE;
00052   useBeamline = kFALSE;
00053   calibBeamline = kFALSE;
00054   useCTB = kFALSE;
00055   usePCT = kFALSE;
00056   useBTOF = kFALSE;
00057   eval = kFALSE;
00058   nEvTotal=nEvGood=0;
00059   externalFindUse=kTRUE; 
00060   mEvalNtuple = 0;
00061   mEvent = 0;
00062   primV = 0;
00063   theFinder = 0;
00064   minTracks = 0;
00065 }
00066 //_____________________________________________________________________________
00067 StGenericVertexMaker::~StGenericVertexMaker()
00068 {
00069 
00070   SafeDelete(theFinder);
00071 
00072 }
00073 
00089 Int_t StGenericVertexMaker::Init()
00090 {
00091   // setup params
00092   useITTF       = IAttr("ITTF");
00093   useBeamline   = IAttr("beamLine");
00094   calibBeamline = IAttr("calibBeamline");
00095   useCTB        = IAttr("CTB");
00096   usePCT        = IAttr("PCT");
00097   useBTOF       = IAttr("BTOF");
00098   eval          = IAttr("eval");
00099   minTracks     = IAttr("minTracks");
00100 
00101   Bool_t isMinuit=kFALSE;
00102 
00103   if ( IAttr("VFMinuit") || IAttr("VFMinuit2") || IAttr("VFMinuit3")){ // 3 versions of Minuit for ranking modes
00104     theFinder= new StMinuitVertexFinder();
00105     if (IAttr("VFMinuit") ) ((StMinuitVertexFinder*) theFinder)->useOldBEMCRank();
00106     if (IAttr("VFMinuit3") ) ((StMinuitVertexFinder*) theFinder)->lowerSplitVtxRank();
00107     if (minTracks > 0) ((StMinuitVertexFinder*) theFinder)->SetMinimumTracks(minTracks);
00108     isMinuit=kTRUE;
00109 
00110   } else if ( IAttr("VFppLMV")){
00111     theFinder= new StppLMVVertexFinder();
00112     theFinder->SetMode(0);                 // this mode is an internal to ppLMV option switch
00113 
00114   } else if ( IAttr("VFppLMV5")){
00115     theFinder= new StppLMVVertexFinder();
00116     theFinder->SetMode(1);                 // this mode is an internal to ppLMV option switch
00117 
00118   } else if ( IAttr("VFPPV") ||  IAttr("VFPPVnoCTB")){ // 2 version of PPV w/ & w/o CTB
00119     LOG_INFO << "StGenericVertexMaker::Init: uses PPVertex finder"<<  endm;
00120     theFinder= new StPPVertexFinder();
00121     if ( IAttr("VFPPVnoCTB")) ((StPPVertexFinder*) theFinder)->useCTB(kFALSE);  
00122     if(GetMaker("emcY2")) {//very dirty, but detects if it is M-C or real data
00123       ((StPPVertexFinder*) theFinder)->setMC(kTRUE);
00124     }
00125 
00126   } else if ( IAttr("VFFV") || IAttr("VFMCE")) {
00127       theFinder = new StFixedVertexFinder();
00128       if (IAttr("VFMCE")){
00129         LOG_INFO << "StGenericVertexMaker::Init: fixed vertex using MC vertex" << endm;
00130         theFinder->SetMode(1);
00131       } else {
00132         LOG_INFO << "StGenericVertexMaker::Init: fixed vertex 'finder' selected" << endm;
00133       }
00134 
00135   } else {
00136     // Later, this would NEVER make multiple possible vertex
00137     // finder unlike for option 0x1 .
00138     theFinder= new StMinuitVertexFinder();
00139     isMinuit=kTRUE;
00140 
00141   }
00142 
00143   theFinder->UsePCT(usePCT);
00144   theFinder->UseBTOF(useBTOF);
00145 
00146   if (calibBeamline) theFinder->CalibBeamLine();
00147 
00148   if(isMinuit) { // this is ugly, one should abort at 'else' above, Jan
00149     if (useITTF)  ((StMinuitVertexFinder*)theFinder)->DoUseITTF();
00150     if (useCTB) ((StMinuitVertexFinder*)theFinder)->CTBforSeed();
00151   } else {
00152     assert(!eval); // current implementation support only Minuit Vertex finder, JB 
00153   }
00154  
00155   //    theFinder->UseVertexConstraint(-0.265,0.4088,-0.00135,0.0004333,0.0001);
00156   //theFinder->UseVertexConstraint(0.0,0.0,0.0,0.0,0.0001);
00157   if (eval) mEvalNtuple = new TNtuple("results","results","thX:thY:thZ:thStat:goodGlob:evX:evY:evZ:evStat:nPrim:nCTB:geantX:geantY:geantZ");
00158 
00159   theFinder->Init();
00160   return StMaker::Init();
00161 }
00162 
00163 //_____________________________________________________________________________
00164 void StGenericVertexMaker::Clear(const char* opt){
00165   LOG_INFO <<" StGenericVertexMaker::Clear()"<<endm; 
00166   theFinder->Clear();
00167 }
00168 
00169 
00170 //_____________________________________________________________________________
00171 
00180 Int_t StGenericVertexMaker::InitRun(int runnumber){
00181   theFinder->InitRun(runnumber);
00182   if (useBeamline) {
00183      double x0 = 0.;
00184      double y0 = 0.;
00185      double dxdz = 0.;
00186      double dydz = 0.;
00187 
00188      // Get Current Beam Line Constraint from database
00189      TDataSet* dbDataSet = this->GetDataBase("Calibrations/rhic");
00190 
00191      if (dbDataSet) {
00192        vertexSeed_st* vSeed = ((St_vertexSeed*) (dbDataSet->FindObject("vertexSeed")))->GetTable();
00193        if(vSeed==0){
00194          LOG_ERROR << "StGenericVertexMaker -- No 'vertexSeed' table in Database for beamline, makse no sens to proceed, Jan" << endm;
00195          assert(1==2);
00196        }
00197        
00198        
00199        x0 = vSeed->x0;
00200        y0 = vSeed->y0;
00201        dxdz = vSeed->dxdz;
00202        dydz = vSeed->dydz;
00203      } else {
00204        LOG_ERROR << "StGenericVertexMaker -- No 'Calibrations/rhic' Database for beamline, makse no sens to proceed, Jan" << endm;
00205        assert(1==2);
00206      }
00207      LOG_INFO << "BeamLine Constraint: " << endm;
00208      LOG_INFO << "x(z) = " << x0 << " + " << dxdz << " * z" << endm;
00209      LOG_INFO << "y(z) = " << y0 << " + " << dydz << " * z" << endm;
00210      theFinder->UseVertexConstraint(x0,y0,dxdz,dydz,0.0001);
00211   }
00212   return StMaker::InitRun(runnumber);
00213 }
00214 
00215 
00216 //_____________________________________________________________________________
00217 Int_t StGenericVertexMaker::Finish()
00218 {
00219 
00220   LOG_INFO << "StGenericVertexMaker::Finish " <<GetName() <<endm;
00221   LOG_INFO << " Total events: " << nEvTotal << endm;
00222   LOG_INFO << " Good events:  " << nEvGood  << endm;
00223 
00224 
00225   //LSB TODO Leave this for now. Should really be using STAR/ROOT I/O scheme?
00226   if (eval) {
00227    TFile out("MinuitVertexEval.root","RECREATE");
00228    mEvalNtuple->Write();
00229    out.Close();
00230   }
00231   
00232   if(theFinder) theFinder->Finish();  
00233   return  kStOK;
00234 }
00235 
00236 //_____________________________________________________________________________
00237 //_____________________________________________________________________________
00238 //_____________________________________________________________________________
00239 Bool_t StGenericVertexMaker::DoFit(){
00240   StThreeVectorD myvertex;
00241 
00242   StEvent *event = (StEvent *) GetInputDS("StEvent");
00243   assert(event);
00244 
00245   if (theFinder->fit(event)) {
00246     theFinder->printInfo();
00247   }  else {
00248     LOG_ERROR << "StGenericVertexMaker::DoFit: vertex fit failed, no vertex." << endm;
00249     return kFALSE;
00250   }
00251 
00252   return kTRUE;
00253 
00254 }
00255 //_____________________________________________________________________________
00256 //_____________________________________________________________________________
00257 //_____________________________________________________________________________
00258 
00259 
00260 //_____________________________________________________________________________
00261 Int_t StGenericVertexMaker::Make()
00262 {
00263   nEvTotal++;
00264   primV  = NULL;
00265   mEvent = NULL;
00266   mEvent = (StEvent *)GetInputDS("StEvent");
00267   LOG_DEBUG << "StGenericVertexMaker::Make: StEvent pointer " << mEvent << endm;
00268   LOG_DEBUG << "StGenericVertexMaker::Make: external find use " << externalFindUse << endm;
00269 
00270   if(!externalFindUse){
00271     DoFit();
00272   }
00273 
00274 //   //For testing purposes
00275 //   theFinder->DoNotUseITTF();
00276 //   if (theFinder->fit(mEvent)) theFinder->printInfo();
00277 //   if (theFinder->status()!=-1) theFinder->FillStEvent(mEvent);
00278 
00279 //   theFinder->DoUseITTF();
00280 //   if (theFinder->fit(mEvent)) theFinder->printInfo();
00281 //   if (theFinder->status()!=-1) {
00282 //     theFinder->FillStEvent(mEvent);
00283 //     nEvGood++;
00284 //   }
00285 
00286 
00287   if (eval)MakeEvalNtuple();
00288 
00289   if(!externalFindUse){
00291     if (theFinder->size()>0){
00292       theFinder->FillStEvent(mEvent);
00293       nEvGood++;
00294     }
00295   }
00296   return kStOK;
00297 }
00298 
00299 //-----------------------------------------------------------------------------
00300 
00301 void StGenericVertexMaker::MakeEvalNtuple(){ // only for Minuit vertex finder
00302 
00303   // get geant vertex
00304   St_DataSet *gds=GetDataSet("geant");
00305   St_g2t_vertex  *g2t_ver=0;
00306   g2t_vertex_st *gver=0;
00307   if(gds)  g2t_ver=( St_g2t_vertex *)gds->Find("g2t_vertex");
00308   if(g2t_ver)gver=g2t_ver->GetTable();
00309 
00310   double gx = -999.;
00311   double gy = -999.;
00312   double gz = -999.;
00313 
00314   if(gver) {
00315     gx=gver->ge_x[0];
00316     gy=gver->ge_x[1];
00317     gz=gver->ge_x[2];
00318   }
00319 
00320 
00321   //  G E T     P R I M     V E R T E X
00322   primV=mEvent->primaryVertex();
00323   int nCtb= ((StMinuitVertexFinder*)theFinder)->NCtbMatches(); 
00324   int stat= ((StMinuitVertexFinder*)theFinder)->statusMin();
00325   
00326   if(!primV) {
00327     LOG_INFO <<"primaryVertex()=NULL"<<endm;
00328     // why would one access x,y,z of the vertex if it is not found, Jan ???
00329     float x=999,y=999,z=999;
00330     mEvalNtuple->Fill(x,y,z,stat,mEvent->summary()->numberOfGoodTracks(),-999.,-999.,-999.,-999.,-999.,nCtb,gx,gy,gz);
00331     }
00332   else {
00333      LOG_INFO << Form("primaryVertex()= %f, %f %f, nTracks=%d\n",primV->position().x(),primV->position().y(),primV->position().z(),primV->numberOfDaughters())<<endm;
00334      mEvalNtuple->Fill(primV->position().x(),primV->position().y(),primV->position().z(),stat               ,mEvent->summary()->numberOfGoodTracks(),primV->position().x(),primV->position().y(),primV->position().z(),primV->flag(),primV->numberOfDaughters(), nCtb,gx,gy,gz);
00335   }
00336 }
00337 
00338 //____________________________________________________________________________