StRoot: StMcEventMaker/StMcEventMaker.cxx Source File

00001 
00271 #include "Stiostream.h"
00272 #include <stdlib.h>
00273 #include <string>
00274 #include <algorithm>
00275 #ifndef ST_NO_NAMESPACES
00276 using std::string;
00277 using std::sort;
00278 using std::find;
00279 #endif
00280 
00281 #include "TStyle.h"
00282 #include "TCanvas.h"
00283 #include "StMcEventMaker.h"
00284 
00285 #include "PhysicalConstants.h"
00286 #include "SystemOfUnits.h"
00287 #include "StGlobals.hh"
00288 #include "StMessMgr.h"
00289 #include "StMemoryInfo.hh"
00290 #include "StTimer.hh"
00291 
00292 #include "StThreeVectorF.hh"
00293 #include "StParticleDefinition.hh"
00294 
00295 #include "St_DataSet.h"
00296 #include "St_DataSetIter.h"
00297 
00298 #include "tables/St_g2t_event_Table.h"
00299 #include "tables/St_g2t_ftp_hit_Table.h"
00300 #include "tables/St_g2t_rch_hit_Table.h"
00301 #include "tables/St_g2t_ctf_hit_Table.h"
00302 #include "tables/St_g2t_mtd_hit_Table.h"
00303 #include "tables/St_g2t_svt_hit_Table.h"
00304 #include "tables/St_g2t_ssd_hit_Table.h"
00305 #include "tables/St_g2t_tpc_hit_Table.h"
00306 #include "tables/St_g2t_emc_hit_Table.h"
00307 #include "tables/St_g2t_pix_hit_Table.h"
00308 #include "tables/St_g2t_ist_hit_Table.h"
00309 #include "tables/St_g2t_fgt_hit_Table.h"
00310 #include "tables/St_g2t_etr_hit_Table.h"
00311 #include "tables/St_g2t_track_Table.h"
00312 #include "tables/St_g2t_vertex_Table.h"
00313 #include "tables/St_particle_Table.h"
00314 
00315 #include "StMcEventTypes.hh"
00316 
00317 #include "StEmcUtil/geometry/StEmcGeom.h" // For Barrel Emc
00318 
00319 //#include "StEEmcUtil/EEmcGeom/EEmcGeomDefs.h" // For Endcap Emc
00320 #include "StEEmcUtil/EEmcMC/EEmcMCData.h"     // geant unpacker, EEmcHit
00321 
00322 static double vertexCut = .0000025; // 25 nm (lifetime of the pi0)
00323 struct vertexFlag {
00324               StMcVertex* vtx;
00325               int primaryFlag; };
00326 
00327 static const char rcsid[] = "$Id: StMcEventMaker.cxx,v 1.76 2012/03/22 01:20:55 perev Exp $";
00328 ClassImp(StMcEventMaker)
00329 #define AddHit2Track(G2Type,DET) \
00330   Int_t iTrkId = ( G2Type ## HitTable[ihit].track_p) - 1;       \
00331   if (iTrkId >= 0 && iTrkId < NTracks ) {                       \
00332     th->setParentTrack(ttemp[iTrkId]);                          \
00333     ttemp[iTrkId]->add ## DET ## Hit(th);                       \
00334   }
00335 #define AddHits(G2Type,det,DET)                                 \
00336   if (doUse ## DET) {                                                   \
00337     if (g2t_ ## G2Type ## _hitTablePointer) {                           \
00338       StMc ## DET ## Hit* th = 0;                                       \
00339       Int_t  NHits = g2t_ ## G2Type ## _hitTablePointer->GetNRows();    \
00340       for(Int_t ihit=0; ihit<NHits; ihit++) {                           \
00341         th = new StMc ## DET ## Hit(& G2Type ## HitTable[ihit]);        \
00342         mCurrentMcEvent-> det ## HitCollection()->addHit(th);           \
00343         AddHit2Track(G2Type,DET);                                       \
00344       }                                                                 \
00345       if (Debug())  cout << "Filled " << mCurrentMcEvent-> det ## HitCollection()->numberOfHits() << #DET << " Hits" << endl; \
00346     } else if (Debug()) cout << "No " << #DET << "/" << #det << " Hits in this event" << endl; \
00347   }
00348 //_____________________________________________________________________________
00349 
00350     
00351 //_____________________________________________________________________________
00352 StMcEventMaker::StMcEventMaker(const char*name, const char * title) :
00353     StMaker(name,title),
00354     doPrintEventInfo (kFALSE),  
00355     doPrintMemoryInfo(kFALSE),  
00356     doPrintCpuInfo   (kFALSE), 
00357     doUseTpc         (kTRUE),
00358     doUseSvt         (kTRUE),
00359     doUseSsd         (kTRUE),
00360     doUseFtpc        (kTRUE),
00361     doUseRich        (kTRUE),
00362     doUseBemc        (kTRUE),
00363     doUseBsmd        (kTRUE),
00364     doUseCtb         (kTRUE),
00365     doUseTofp        (kFALSE),
00366     doUseTof         (kTRUE),
00367     doUseMtd         (kTRUE),
00368     doUseEemc        (kTRUE),
00369     doUseFpd         (kTRUE),
00370     doUseFsc         (kTRUE),
00371     doUsePixel       (kTRUE),
00372     doUseIst         (kTRUE),
00373     doUseFgt         (kTRUE),
00374     doUseEtr         (kTRUE),
00375     ttemp(),
00376     ttempParticle(),
00377     mCurrentMcEvent(0)    
00378 {
00379     // StMcEventMaker - constructor
00380     // - set all pointers defined in the header file to zero
00381 
00382 }
00383 //_____________________________________________________________________________
00384 StMcEventMaker::~StMcEventMaker()
00385 {
00386     // StMcEventMaker - destructor
00387     if (Debug()>=2) cout << "Inside ReaderMaker Destructor" << endl;
00388     //    SafeDelete(mCurrentMcEvent);  //
00389 
00390 }
00391 
00392 
00393 
00394 //_____________________________________________________________________________
00395 void StMcEventMaker::Clear(const char*)
00396 {
00397     if (doPrintMemoryInfo) StMemoryInfo::instance()->snapshot();    
00398     // StMcEventMaker - Clear,
00399     if (mCurrentMcEvent) {
00400       //        delete mCurrentMcEvent;
00401         mCurrentMcEvent = 0;
00402     }
00403     if (doPrintMemoryInfo) {
00404         StMemoryInfo::instance()->snapshot();
00405         StMemoryInfo::instance()->print();
00406     }
00407     StMaker::Clear();
00408 }
00409 
00410 
00411 //_____________________________________________________________________________
00412 Int_t StMcEventMaker::Finish()
00413 {
00414     // StMcEventMaker - Finish, Draw histograms if SetDraw true
00415 
00416     // Right now I'm not doing any histograms, later on, I would need to uncomment
00417     // the next line, and add a DrawHists() method.  Look in St_QA_Maker.cxx
00418     //if (drawinit)  DrawHists();
00419   return StMaker::Finish();
00420 }
00421 
00422 
00423 //_____________________________________________________________________________
00424 Int_t StMcEventMaker::Init()
00425 {
00426     if (Debug()>=1) cout << "This is StMcEventMaker::Init() - by Ming" << endl;
00427     return StMaker::Init();
00428 }
00429 
00430 //_____________________________________________________________________________
00431 Int_t StMcEventMaker::Make()
00432 {
00433     // StMcEventMaker - Make; fill StMcEvent objects
00434     StTimer timer;
00435     if (doPrintCpuInfo) timer.start();
00436     if (doPrintMemoryInfo) StMemoryInfo::instance()->snapshot();
00437     
00438     if (Debug()>=1) cout << "Inside StMcEventMaker::Make()" << endl;
00439     // We're supposed to get the dataset from the chain. I don't know how yet. I think it is:
00440     
00441     const Char_t* geaTmp[3]={"geant","event/geant/Event","bfcTree/geantBranch"};
00442     St_DataSet* dsGeant = 0;
00443     for(UInt_t i=0; i<3; i++){ 
00444       dsGeant = GetDataSet(geaTmp[i]);
00445       if(!dsGeant || !dsGeant->GetList()) {
00446         gMessMgr->Warning() << "Could not find dataset " << geaTmp[i] << endm;
00447         dsGeant = GetDataSet("event/geant/Event");
00448       } else {
00449           if (Debug()) gMessMgr->Info() << "Get Geant info from  dataset " << geaTmp[i] << endm;
00450         break;
00451       }
00452     }
00453     if(!dsGeant) return kStWarn;;
00454     /* 
00455     if(!dsGeant || !dsGeant->GetList()) {
00456         gMessMgr->Warning() << "Could not find dataset geant" << endm;
00457         dsGeant = GetDataSet("event/geant/Event");
00458         if(!dsGeant || !dsGeant->GetList()) {  // Try direct output from Geant
00459           gMessMgr->Warning() << "Could not find dataset event/geant/Event" << endm;
00460           dsGeant = GetDataSet("bfcTree/geantBranch");
00461           if(!dsGeant || !dsGeant->GetList()) {  // Try output from bfc forGeant
00462             gMessMgr->Warning() << "Could not find dataset bfcTree/geantBranch" << endm;
00463             return kStWarn;
00464           } else cout<< "Find dataset bfcTree/geantBranch "<<endl;
00465         }
00466     }
00467     */
00468     // This is done only for one file, though.  I haven't put functionality for
00469     // multiple file handling.  If needed, it should start in the StMcEventReadMacro
00470     // and try to follow the doEvents.C macro to read in multiple files.
00471     // Then, we would need to mimic the nextRootFile() methods and so on.
00472     
00473   
00474     // Now we have the DataSet, but for some reason, we need the Iterator to navigate
00475     
00476     St_DataSetIter geantDstI(dsGeant);
00477       
00478     // Now the Iterator is set up, and this allows us to access the tables
00479     // This is done like so:
00480     // TableClass *instanceOfTableClassPointer = cast to TableClassPointer instanceOfDataSetIter("actual name of table in data set");
00481     
00482     St_g2t_event   *g2t_eventTablePointer   =  (St_g2t_event   *) geantDstI("g2t_event");
00483     St_g2t_vertex  *g2t_vertexTablePointer  =  (St_g2t_vertex  *) geantDstI("g2t_vertex");
00484     St_g2t_track   *g2t_trackTablePointer   =  (St_g2t_track   *) geantDstI("g2t_track");
00485     St_g2t_tpc_hit *g2t_tpc_hitTablePointer =  (St_g2t_tpc_hit *) geantDstI("g2t_tpc_hit");
00486     St_g2t_svt_hit *g2t_svt_hitTablePointer =  (St_g2t_svt_hit *) geantDstI("g2t_svt_hit");
00487     St_g2t_ssd_hit *g2t_ssd_hitTablePointer =  (St_g2t_ssd_hit *) geantDstI("g2t_ssd_hit");
00488     St_g2t_ftp_hit *g2t_ftp_hitTablePointer =  (St_g2t_ftp_hit *) geantDstI("g2t_ftp_hit");
00489     St_g2t_rch_hit *g2t_rch_hitTablePointer =  (St_g2t_rch_hit *) geantDstI("g2t_rch_hit");
00490     St_g2t_emc_hit *g2t_emc_hitTablePointer =  (St_g2t_emc_hit *) geantDstI("g2t_emc_hit");
00491     St_g2t_emc_hit *g2t_smd_hitTablePointer =  (St_g2t_emc_hit *) geantDstI("g2t_smd_hit");
00492     St_g2t_ctf_hit *g2t_ctb_hitTablePointer =  (St_g2t_ctf_hit *) geantDstI("g2t_ctb_hit"); // Added CTB Hits
00493     St_g2t_ctf_hit *g2t_tof_hitTablePointer =  (St_g2t_ctf_hit *) geantDstI("g2t_tof_hit");
00494     St_g2t_ctf_hit *g2t_tfr_hitTablePointer =  (St_g2t_ctf_hit *) geantDstI("g2t_tfr_hit");
00495     St_g2t_mtd_hit *g2t_mtd_hitTablePointer =  (St_g2t_mtd_hit *) geantDstI("g2t_mtd_hit");
00496     St_g2t_emc_hit *g2t_eem_hitTablePointer =  (St_g2t_emc_hit *) geantDstI("g2t_eem_hit");
00497     St_g2t_emc_hit *g2t_esm_hitTablePointer =  (St_g2t_emc_hit *) geantDstI("g2t_esm_hit");
00498     St_g2t_emc_hit *g2t_fpd_hitTablePointer =  (St_g2t_emc_hit *) geantDstI("g2t_fpd_hit"); // Added FPD Hits
00499     St_g2t_emc_hit *g2t_fsc_hitTablePointer =  (St_g2t_emc_hit *) geantDstI("g2t_fsc_hit"); // Added FSC Hits
00500     St_g2t_pix_hit *g2t_pix_hitTablePointer =  (St_g2t_pix_hit *) geantDstI("g2t_pix_hit");
00501     St_g2t_ist_hit *g2t_ist_hitTablePointer =  (St_g2t_ist_hit *) geantDstI("g2t_ist_hit");
00502     St_g2t_fgt_hit *g2t_fgt_hitTablePointer =  (St_g2t_fgt_hit *) geantDstI("g2t_fgt_hit");
00503     St_g2t_etr_hit *g2t_etr_hitTablePointer =  (St_g2t_etr_hit *) geantDstI("g2t_etr_hit");
00504     St_particle    *particleTablePointer    =  (St_particle    *) geantDstI("particle");
00505 
00506     // For backwards compatibility, look for the rch and particle tables also in the dstBranch
00507     TDataSet *dst = GetDataSet("dst");
00508     if (dst) {
00509       St_DataSetIter dstDstI(dst);
00510       if (!particleTablePointer)
00511         particleTablePointer    = (St_particle    *) dstDstI("particle");
00512       if (!g2t_rch_hitTablePointer)
00513         g2t_rch_hitTablePointer = (St_g2t_rch_hit *) dstDstI("g2t_rch_hit");
00514     }
00515     // Now we check if we have the pointer, if we do, then we can access the tables!
00516   
00517     if (g2t_vertexTablePointer && g2t_trackTablePointer){
00518         
00519         //
00520         // g2t_event Table
00521         //
00522         g2t_event_st   *eventTable = 0;
00523         
00524         // Check Pointer
00525         if (g2t_eventTablePointer)
00526             eventTable  = g2t_eventTablePointer->GetTable();
00527         else
00528             if (Debug()) cerr << "Table g2t_event Not found in Dataset " << geantDstI.Pwd()->GetName() << endl;
00529 
00530         // Check Table
00531         if(!g2t_eventTablePointer->GetNRows()) {
00532             cout << "Event Table is EMPTY!! " << endl;
00533             PR(g2t_eventTablePointer->GetNRows());
00534             PR(eventTable)
00535         }
00536         //
00537         // Vertex, Track and table don't have problems normally.
00538         //
00539         g2t_vertex_st  *vertexTable = g2t_vertexTablePointer->GetTable();
00540         g2t_track_st   *trackTable  = g2t_trackTablePointer->GetTable();
00541 
00542         //
00543         // Tpc Hit Table (might not be there for photon events
00544         //
00545         g2t_tpc_hit_st *tpcHitTable = 0;
00546         if (g2t_tpc_hitTablePointer)
00547             tpcHitTable = g2t_tpc_hitTablePointer->GetTable();
00548         else
00549             if (Debug()) cerr << "Table g2t_tpc_hit Not found in Dataset " << geantDstI.Pwd()->GetName() << endl;
00550         //
00551         // Ftpc Hit Table
00552         //
00553         g2t_ftp_hit_st *ftpHitTable = 0;
00554         if (g2t_ftp_hitTablePointer)
00555             ftpHitTable = g2t_ftp_hitTablePointer->GetTable();
00556         else 
00557             if (Debug()) cerr << "Table g2t_ftp_hit Not found in Dataset " << geantDstI.Pwd()->GetName() << endl;
00558 
00559         //
00560         // Svt Hit Table
00561         //
00562         g2t_svt_hit_st *svtHitTable = 0;
00563         if (g2t_svt_hitTablePointer)
00564             svtHitTable = g2t_svt_hitTablePointer->GetTable();
00565         else
00566             if (Debug()) cerr << "Table g2t_svt_hit Not found in Dataset " << geantDstI.Pwd()->GetName() << endl;
00567         
00568         //
00569         // Ssd Hit Table
00570         //
00571         g2t_ssd_hit_st *ssdHitTable = 0;
00572         if (g2t_ssd_hitTablePointer)
00573             ssdHitTable = g2t_ssd_hitTablePointer->GetTable();
00574         else
00575             if (Debug()) cerr << "Table g2t_ssd_hit Not found in Dataset " << geantDstI.Pwd()->GetName() << endl;
00576 
00577         //
00578         // Rich Hit Table
00579         //
00580         g2t_rch_hit_st *rchHitTable = 0;
00581         if (g2t_rch_hitTablePointer)
00582             rchHitTable = g2t_rch_hitTablePointer->GetTable();
00583         else
00584             if (Debug()) cerr << "Table g2t_rch_hit Not found in Dataset " << geantDstI.Pwd()->GetName() << endl;
00585         //
00586         // Ctb Hit Table
00587         //
00588         g2t_ctf_hit_st *ctbHitTable = 0;
00589         if (g2t_ctb_hitTablePointer)
00590             ctbHitTable = g2t_ctb_hitTablePointer->GetTable();
00591         else
00592             if (Debug()) cerr << "Table g2t_rch_hit Not found in Dataset " << geantDstI.Pwd()->GetName() << endl;
00593 
00594         // TOF Hit Tables
00595         //
00596         g2t_ctf_hit_st *tofHitTable = 0;
00597         if (g2t_tof_hitTablePointer)
00598             tofHitTable = g2t_tof_hitTablePointer->GetTable();
00599         else
00600             if (Debug()) cerr << "Table g2t_tof_hit Not found in Dataset " << geantDstI.Pwd()->GetName() << endl;
00601         
00602         g2t_ctf_hit_st *tfrHitTable = 0;
00603         if (g2t_tfr_hitTablePointer)
00604             tfrHitTable = g2t_tfr_hitTablePointer->GetTable();
00605         else
00606             if (Debug()) cerr << "Table g2t_tfr_hit Not found in Dataset " << geantDstI.Pwd()->GetName() << endl;
00607 
00608         // MTD Hit Tables
00609         //
00610         g2t_mtd_hit_st *mtdHitTable = 0;
00611         if (g2t_mtd_hitTablePointer)
00612             mtdHitTable = g2t_mtd_hitTablePointer->GetTable();
00613         else
00614             if (Debug()) cerr << "Table g2t_mtd_hit Not found in Dataset " << geantDstI.Pwd()->GetName() << endl;
00615 
00616         //
00617         // BEMC and BPRS Hit Table
00618         //
00619         g2t_emc_hit_st *emcHitTable = 0;
00620         if (g2t_emc_hitTablePointer)
00621             emcHitTable = g2t_emc_hitTablePointer->GetTable();
00622         else
00623             if (Debug()) cerr << "Table g2t_emc_hit Not found in Dataset " << geantDstI.Pwd()->GetName() << endl;
00624 
00625         //
00626         // BSMDE and BSMDP Hit Table
00627         //
00628         g2t_emc_hit_st *smdHitTable = 0;
00629         if (g2t_smd_hitTablePointer)
00630             smdHitTable = g2t_smd_hitTablePointer->GetTable();
00631         else
00632             if (Debug()) cerr << "Table g2t_smd_hit Not found in Dataset " << geantDstI.Pwd()->GetName() << endl;
00633 
00634         //
00635         // EEMC and EPRS Hit Table
00636         //
00637         g2t_emc_hit_st *eemHitTable = 0;
00638         if (g2t_eem_hitTablePointer)
00639             eemHitTable = g2t_eem_hitTablePointer->GetTable();
00640         else
00641             if (Debug()) cerr << "Table g2t_eem_hit Not found in Dataset " << geantDstI.Pwd()->GetName() << endl;
00642 
00643         //
00644         // ESMDU and ESMDV Hit Table
00645         //
00646         g2t_emc_hit_st *esmHitTable = 0;
00647         if (g2t_esm_hitTablePointer)
00648             esmHitTable = g2t_esm_hitTablePointer->GetTable();
00649         else
00650             if (Debug()) cerr << "Table g2t_esm_hit Not found in Dataset " << geantDstI.Pwd()->GetName() << endl;
00651 
00652         //
00653         // FPD Hit Table
00654         //
00655         g2t_emc_hit_st* fpdHitTable = 0;
00656         if (g2t_fpd_hitTablePointer)
00657           fpdHitTable = g2t_fpd_hitTablePointer->GetTable();
00658         else
00659           if (Debug()) cout << "Table g2t_fpd_hit Not found in Dataset " << geantDstI.Pwd()->GetName() << endl;
00660 
00661         //
00662         // FSC Hit Table
00663         //
00664         g2t_emc_hit_st* fscHitTable = 0;
00665         if (g2t_fsc_hitTablePointer)
00666           fscHitTable = g2t_fsc_hitTablePointer->GetTable();
00667         else
00668           if (Debug()) cout << "Table g2t_fsc_hit Not found in Dataset " << geantDstI.Pwd()->GetName() << endl;
00669 
00670         //      
00671         // Pixel Hit Table
00672         //
00673         g2t_pix_hit_st *pixHitTable=0;
00674         if (g2t_pix_hitTablePointer)
00675             pixHitTable = g2t_pix_hitTablePointer->GetTable();
00676         else 
00677             if (Debug()) cerr << "Table g2t_pix_hit Not found in Dataset " << geantDstI.Pwd()->GetName() << endl;
00678 
00679         //      
00680         // Ist Hit Table
00681         //
00682         g2t_ist_hit_st *istHitTable=0;
00683         if (g2t_ist_hitTablePointer)
00684             istHitTable = g2t_ist_hitTablePointer->GetTable();
00685         else 
00686             if (Debug()) cerr << "Table g2t_ist_hit Not found in Dataset " << geantDstI.Pwd()->GetName() << endl;
00687         //      
00688         // Fgt Hit Table
00689         //
00690         g2t_fgt_hit_st *fgtHitTable=0;
00691         if (g2t_fgt_hitTablePointer)
00692             fgtHitTable = g2t_fgt_hitTablePointer->GetTable();
00693         else 
00694             if (Debug()) cerr << "Table g2t_fgt_hit Not found in Dataset " << geantDstI.Pwd()->GetName() << endl;
00695 
00696         //      
00697         // Etr Hit Table
00698         //
00699         g2t_etr_hit_st *etrHitTable=0;
00700         if (g2t_etr_hitTablePointer)
00701             etrHitTable = g2t_etr_hitTablePointer->GetTable();
00702             if (Debug()) cerr << "Table g2t_etr_hit found in Dataset " << geantDstI.Pwd()->GetName() << endl;
00703         else 
00704             if (Debug()) cerr << "Table g2t_etr_hit Not found in Dataset " << geantDstI.Pwd()->GetName() << endl;
00705 
00706         //
00707         // particle Table
00708         //
00709         particle_st *particleTable = 0;
00710         if (particleTablePointer)
00711             particleTable = particleTablePointer->GetTable();
00712         else
00713             if (Debug()) cerr << "Table particle Not found in Dataset " << geantDstI.Pwd()->GetName() << endl;
00714        
00715        // Before filling StMcEvent, we can check whether we can actually
00716        // access the tables.
00717        
00718 //        cout << "Event Table Examples: " << endl;
00719 //        cout << "eg_label: " << eventTable->eg_label << endl; 
00720 //        cout << "n_event : " << eventTable->n_event << endl;
00721 //        cout << "n_run   : " << eventTable->n_run << endl;
00722 //        cout << "n_part_prot_west: " <<  eventTable->n_part_prot_west << endl;
00723 //        cout << "n_part_neut_west: " <<  eventTable->n_part_neut_west << endl;
00724 //        cout << "n_part_prot_east: " <<  eventTable->n_part_prot_east << endl;
00725 //        cout << "n_part_neut_east: " <<  eventTable->n_part_neut_east << endl;
00726 
00727 //        cout << "Vertex Table Examples:" << endl;
00728 //        cout << "ge_x[0] :" << vertexTable[0].ge_x[0] << endl;
00729 //        cout << "ge_x[1] :" << vertexTable[0].ge_x[1] << endl;
00730 //        cout << "ge_x[2] :" << vertexTable[0].ge_x[2] << endl;
00731 
00732 //        cout << "Track Table Examples:" << endl;
00733 //        cout << "p[0] :" << trackTable[0].p[0] << endl;
00734 //        cout << "p[1] :" << trackTable[0].p[1] << endl;
00735 //        cout << "p[2] :" << trackTable[0].p[2] << endl;
00736           
00737 //        cout << "Tpc Hit Table Examples:" << endl;
00738 //        cout << "p[0] :" << tpcHitTable[0].p[0] << endl;
00739 //        cout << "p[1] :" << tpcHitTable[0].p[1] << endl;
00740 //        cout << "p[2] :" << tpcHitTable[0].p[2] << endl;
00741 
00742 //        cout << "Svt Hit Table Examples:" << endl;
00743 //        cout << "p[0] :" << svtHitTable[0].p[0] << endl;
00744 //        cout << "p[1] :" << svtHitTable[0].p[1] << endl;
00745 //        cout << "p[2] :" << svtHitTable[0].p[2] << endl;
00746 
00747 //        cout << "Ftpc Hit Table Examples:" << endl;
00748 //        cout << "p[0] :" << ftpHitTable[0].p[0] << endl;
00749 //        cout << "p[1] :" << ftpHitTable[0].p[1] << endl;
00750 //        cout << "p[2] :" << ftpHitTable[0].p[2] << endl;
00751           
00752 //        cout << "Rich Hit Table Examples:" << endl;
00753 //        cout << "p[0] :" << rchHitTable[0].p[0] << endl;
00754 //        cout << "p[1] :" << rchHitTable[0].p[1] << endl;
00755 //        cout << "p[2] :" << rchHitTable[0].p[2] << endl;
00756           
00757         // Ok, now we have the g2t tables for this event, now we can create the
00758         // StMcEvent with them.
00759         
00760         // Now Here goes the scheme Mike Lisa and I cooked up
00761         
00762         //______________________________________________________________________
00763         // Step 1 - fill the StMcEvent, already created in the header file
00764         
00765         if (eventTable)
00766             mCurrentMcEvent = new StMcEvent(eventTable);
00767         else {
00768             gMessMgr->Warning() << "StMcEventMaker::Make(): g2t_event Table not found.  Using default constructor." << endm;
00769             mCurrentMcEvent = new StMcEvent;  
00770         }
00771         
00772         if (mCurrentMcEvent) {
00773             if (Debug()) cout << "****  Created new StMcEvent, Event No. " << mCurrentMcEvent->eventNumber() << endl;
00774         }
00775         else {
00776             gMessMgr->Warning() << "Could not create StMcEvent! Exit from StMcEventMaker." << endm;
00777             return kStWarn;
00778         }
00779         AddData(mCurrentMcEvent);
00780         //______________________________________________________________________
00781         // Step 2 - Fill Vertices - we do not fill parent/daughters until Step 3
00782         
00783         long NVertices = g2t_vertexTablePointer->GetNRows();
00784           
00785 #ifndef ST_NO_TEMPLATE_DEF_ARGS
00786         vector<vertexFlag> vtemp(NVertices); // Temporary array for Step 3
00787 #else
00788         vector<vertexFlag, allocator<vertexFlag> > vtemp(NVertices);
00789 #endif
00790        
00791         if (Debug()>=1) cout << "Preparing to process and fill VERTEX information ....." << endl;
00792         StMcVertex* v = 0;
00793         
00794         
00795         // First vertex is PRIMARY
00796         
00797         v = new StMcVertex(&(vertexTable[0]));
00798         mCurrentMcEvent->vertices().push_back(v);
00799         vtemp[vertexTable[0].id - 1].vtx = v;
00800         vtemp[vertexTable[0].id - 1].primaryFlag = 1;
00801         mCurrentMcEvent->setPrimaryVertex(v);
00802         
00803         StThreeVectorF primVertexPos = v->position();
00804         
00805         StThreeVectorF testVertexPos;
00806         int nThrownVertices = 0;
00807         for (long ivtx=1; ivtx<NVertices; ivtx++)
00808             {    
00809                 v = new StMcVertex(&(vertexTable[ivtx]));
00810                 
00811                 // get the position of the current vertex
00812                 testVertexPos = v->position();
00813                 
00814                 if (vertexTable[ivtx].eg_label == 0 ||
00815                     abs(primVertexPos - testVertexPos) > vertexCut ) {
00816                     // GEANT vertex or (generator vertex that satisfies cut) ...
00817                     mCurrentMcEvent->vertices().push_back(v);  // adds vertex v to master collection
00818                     vtemp[vertexTable[ivtx].id - 1].primaryFlag = 0;
00819                     vtemp[vertexTable[ivtx].id - 1].vtx = v;
00820                 }
00821                 else { // These "vertices" are the same as primary, so flag them
00822                     nThrownVertices++;
00823                     vtemp[vertexTable[ivtx].id - 1].primaryFlag = 1;
00824                     
00825                     delete v;
00826                     vtemp[vertexTable[ivtx].id - 1].vtx = mCurrentMcEvent->primaryVertex();
00827                 }
00828             }
00829         if (nThrownVertices) 
00830             gMessMgr->Warning() << "StMcEventMaker::Make(): Throwing " << nThrownVertices
00831                                 << " that are the same as the primary vertex." << endm;
00832        
00833         //______________________________________________________________________
00834         // Step 3 - Fill Tracks - we do not fill associated hits until Step 4
00835         
00836         long NTracks = g2t_trackTablePointer->GetNRows();
00837         size_t usedTracksG2t = 0;
00838         long NGeneratorTracks = (particleTablePointer) ? particleTablePointer->GetNRows() : 0;
00839         size_t usedTracksEvGen = 0;
00840         ttemp.resize(NTracks);
00841         ttempParticle.resize(NGeneratorTracks);
00842         if (Debug()>=1) cout << "Preparing to process and fill TRACK information ....." << endl;
00843         StMcTrack* egTrk = 0;
00844         size_t nParticlesInBothTables = 0;
00845         if (Debug()>=1) cout << "Event Generator Tracks..." << endl;
00846 
00847         // The filling of the event generator track assumes that if we encounter a particle that
00848         // has a mother, the mother should ALREADY HAVE BEEN CREATED, i.e. that the mother indices
00849         // of the particles in the particle table is never more than the current index.
00850 
00851         long motherIndex = -1;  // Set it to some unused number. 
00852         {for (long gtrk=0; gtrk<NGeneratorTracks; gtrk++) {
00853             //      if (particleTable[gtrk].isthep > 3) continue; // skip comment fields
00854             egTrk = new StMcTrack(&(particleTable[gtrk]));
00855             egTrk->setEventGenLabel(gtrk+1);
00856             egTrk->setPrimary(0);
00857             ttempParticle[gtrk] = egTrk;
00858             mCurrentMcEvent->tracks().push_back(egTrk); // adds track egTrk to master collection 
00859             usedTracksEvGen++;
00860             
00861         }} // Generator Tracks Loop
00862         
00863         StMcTrack* t = 0;
00864         long iStartVtxId = 0;
00865         long iStopVtxId = 0;
00866         long iItrmdVtxId = 0;
00867         
00868         int nThrownTracks = 0;
00869         if (Debug()>=1) cout << "g2t Tracks..." << endl;
00870         
00871         {for (long itrk=0; itrk<NTracks; itrk++) {
00872             iStopVtxId = (trackTable[itrk].stop_vertex_p) - 1;
00873                                 
00874             if (iStopVtxId >= 0) {
00875                 if (vtemp[iStopVtxId].primaryFlag == 1) {
00876                     
00877                     nThrownTracks++;
00878                     continue; // This skips until the next itrk
00879                     
00880                 }
00881             }    
00882 //              Fix(hack)  garbage (VP)
00883             if (trackTable[itrk].e<trackTable[itrk].ptot) 
00884                 trackTable[itrk].ptot=trackTable[itrk].e;
00885 //              Fix(hack) tracks with garbage geand id(VP)
00886             if (trackTable[itrk].ge_pid<0)      trackTable[itrk].ge_pid=0;
00887             if (trackTable[itrk].ge_pid>0xffff) trackTable[itrk].ge_pid=0;
00888 
00889             t = new StMcTrack(&(trackTable[itrk]));
00890             usedTracksG2t++;
00891             ttemp[ trackTable[itrk].id - 1 ] = t; // This we do for all accepted tracks
00892             
00893             
00894             mCurrentMcEvent->tracks().push_back(t); // adds track t to master collection
00895             
00896             // point track to its stop vertex,
00897             // and tell stop vertex that this is its parent
00898             if (iStopVtxId >= 0) {
00899                 t->setStopVertex(vtemp[iStopVtxId].vtx);
00900                 vtemp[iStopVtxId].vtx->setParent(t);
00901             }
00902             //cout << "Established relation btw track & STOP vertex" << endl;
00903             
00904             // point track to its parent vertex,
00905             // and tell parent vertex that this is its daughter
00906             
00907             iStartVtxId = trackTable[itrk].start_vertex_p - 1;
00908             
00909             v = vtemp[iStartVtxId].vtx; // This should already know the right vertex.
00910             
00911             t->setStartVertex(v);
00912             t->setPrimary(t->startVertex() == mCurrentMcEvent->primaryVertex() ? kTRUE : kFALSE);
00913             v->addDaughter(t);
00914             //cout << "Established relation btw track & START vertex" << endl;
00915             
00916             // now fill track's intermediate vertex collection,
00917             // and tell those vertices their parents
00918             
00919             iItrmdVtxId = (trackTable[itrk].itrmd_vertex_p) - 1;
00920             
00921             if (iItrmdVtxId >= 0) {
00922                 t->intermediateVertices().push_back(vtemp[iItrmdVtxId].vtx);
00923                 
00924                 vtemp[iItrmdVtxId].vtx->setParent(t);
00925                 
00926                 // follow the "linked list" of g2t_vertex table to get the rest
00927                 
00928                 while(vertexTable[iItrmdVtxId].next_itrmd_p != 0) {
00929                     iItrmdVtxId = (vertexTable[iItrmdVtxId].next_itrmd_p) - 1;
00930                     t->intermediateVertices().push_back(vtemp[iItrmdVtxId].vtx);
00931                     vtemp[iItrmdVtxId].vtx->setParent(t);
00932                 }
00933                 
00934                 
00935             } // Intermediate vertices
00936 
00937             // Look in the particle table
00938             // for particles from event generator
00939             long iEventGeneratorLabel = (trackTable[itrk].eg_label) - 1;
00940             if (iEventGeneratorLabel >=0 ) {
00941 
00942                 // Now make sure that this track is really from the table,
00943                 // When embedding, the particle got an eg_label = 99999 even
00944                 // though there was only one entry in the particle table.
00945                 if (iEventGeneratorLabel < NGeneratorTracks) {
00946                     // Track should already be loaded from the particle table
00947                     // i.e. t & ttempParticle[iEventGeneratorLabel] are the same tracks,
00948                     // obtained from different tables.
00949                     // We should rather keep the one in the g2t table, but we
00950                     // need to keep the information of its parentage.
00951                     nParticlesInBothTables++;
00952                     if (Debug()>=2) {
00953                         if (particleTable[iEventGeneratorLabel].jmohep[0] && !(ttempParticle[particleTable[iEventGeneratorLabel].jmohep[0]])) {
00954                             cout << "There should be a parent and there isn't one!\n";
00955                             PR(iEventGeneratorLabel);
00956                             PR(particleTable[iEventGeneratorLabel].jmohep[0]);
00957                             PR(ttempParticle[particleTable[iEventGeneratorLabel].jmohep[0]]);
00958                             
00959                         }
00960                     }
00961                     if (particleTable[iEventGeneratorLabel].jmohep[0])
00962                         t->setParent(ttempParticle[particleTable[iEventGeneratorLabel].jmohep[0]]);
00963                     StMcTrackIterator trkToErase = find (mCurrentMcEvent->tracks().begin(),
00964                                                          mCurrentMcEvent->tracks().end(),
00965                                                          ttempParticle[iEventGeneratorLabel]);
00966                     delete *trkToErase; // if we delete using the iterator, deleting should be done before erasing!
00967                     mCurrentMcEvent->tracks().erase(trkToErase);
00968                     
00969                     ttempParticle[iEventGeneratorLabel] = t;
00970                 }
00971                 
00972             }
00973             else {
00974                 // Track from GEANT, use next_parent_p to get to the parent
00975                 // track.  Use the same scheme as for the particle table.
00976                 motherIndex = trackTable[itrk].next_parent_p;
00977                 if ((motherIndex > 0) && (motherIndex < NTracks)) {
00978                     if (motherIndex > itrk) { 
00979                         if (Debug()) {
00980                             gMessMgr->Warning()
00981                                 << "Wrong ordering!  Track " << itrk+1 << "from particle table: "
00982                                 << "Can't assign mother track " << motherIndex
00983                                 << "because it has not been created yet!" << endm;
00984                         }
00985                     }
00986                     else {t->setParent(ttemp[motherIndex-1]);}
00987             } }
00988             
00989         }} // Track loop
00990         {for (long gtrk=0; gtrk<NGeneratorTracks; gtrk++) {
00991             // Find Mother...
00992             motherIndex = particleTable[gtrk].jmohep[0];
00993             if ((motherIndex > 0) && (motherIndex < NGeneratorTracks)) {
00994                 if (motherIndex > gtrk) {
00995                     if (Debug()) {
00996                         gMessMgr->Warning()
00997                             << "Wrong ordering!  Track " << gtrk+1 << " from particle table: "
00998                             << "Can't assign mother track " << motherIndex
00999                             << " to track with index " << gtrk << endm;
01000                     }
01001                 }
01002                 else ttempParticle[gtrk]->setParent(ttempParticle[motherIndex-1]);
01003                 if (Debug()>=2) {
01004                   cout << "Particle table (generator) track " << gtrk << ", key " << ttempParticle[gtrk]->key() << endl;
01005                   cout << "Mother Index-1 (off-by-1)        " << motherIndex-1 << endl;
01006                   cout << "PDG ID of generator track        " << particleTable[gtrk].idhep << endl;
01007                   cout << "PDG ID of mother track           " << particleTable[motherIndex-1].idhep << endl;
01008                   if (ttempParticle[gtrk]->particleDefinition())
01009                     cout << "particle                         " << ttempParticle[gtrk]->particleDefinition()->name() << endl;
01010                   if (ttempParticle[gtrk]->parent() && ttempParticle[gtrk]->parent()->particleDefinition())
01011                     cout << "parent                           " <<  ttempParticle[gtrk]->parent()->particleDefinition()->name() << endl;
01012                   
01013                   if (motherIndex && !(ttempParticle[gtrk]->parent())) {
01014                       cout << "Error in assigning parent to particle table!\n There should be a parent and there isn't one!\n";
01015                       PR(particleTable[gtrk].jmohep[0]);
01016                       PR(ttempParticle[gtrk]->parent());
01017                   }
01018                 }
01019             }
01020         }}
01021         if (nThrownTracks)
01022             gMessMgr->Warning() << "StMcEventMaker::Make(): Throwing " << nThrownTracks
01023                                 << " whose stop vertex is the same as primary vertex." << endm;
01024         if (Debug()>=2) {
01025             // Check the whole ttempParticle for entries with problems
01026             for (long gtrk=0; gtrk<NGeneratorTracks; gtrk++)
01027                 if (ttempParticle[gtrk]->parent() && ttempParticle[gtrk]->parent() != ttempParticle[particleTable[gtrk].jmohep[0]-1]) {
01028                     cout << "The indexing got screwed up!" << endl;
01029                     PR(ttempParticle[gtrk]->eventGenLabel());
01030                     PR(ttempParticle[gtrk]->key());
01031                     PR(ttempParticle[gtrk]->parent());
01032                     PR(ttempParticle[particleTable[gtrk].jmohep[0]-1]);
01033                     if (ttempParticle[gtrk]->particleDefinition()) cout << "particle " << ttempParticle[gtrk]->particleDefinition()->name() << endl;
01034                     if (ttempParticle[gtrk]->parent() && ttempParticle[gtrk]->parent()->particleDefinition()) cout << "parent   " <<  ttempParticle[gtrk]->parent()->particleDefinition()->name() << endl;
01035 
01036                 }
01037             cout << "Used   tracks from g2t_track table: " << usedTracksG2t << endl;
01038             cout << "Avail. tracks from g2t_track table: " << NTracks       << endl;
01039             cout << "Used   tracks from particle  table: " << usedTracksEvGen  << endl;
01040             cout << "Avail. tracks from particle  table: " << NGeneratorTracks << endl;
01041             cout << "Tracks appearing in both tables   : " << nParticlesInBothTables << endl;
01042             cout << "Total tracks in StMcEvent         : " << mCurrentMcEvent->tracks().size() << endl;
01043         }
01044         vtemp.clear();
01045         ttempParticle.clear();
01046         
01047         //______________________________________________________________________
01048         // Step 4 - Fill Hits
01049                 
01050         if (Debug()) cout << "Preparing to process and fill HIT information ....." << endl;
01051         
01052 
01053         //
01054         // TPC Hits
01055         //
01056         if (doUseTpc) {
01057             if (g2t_tpc_hitTablePointer) {
01058                 StMcTpcHit* th = 0;
01059                 long NHits = g2t_tpc_hitTablePointer->GetNRows();
01060                 long nBadVolId = 0;
01061                 long nPseudoPadrow = 0;
01062                 long ihit;
01063                 for(ihit=0; ihit<NHits; ihit++) {
01064                     if (tpcHitTable[ihit].volume_id < 101 || tpcHitTable[ihit].volume_id > 2445) {
01065                         if (tpcHitTable[ihit].volume_id <= 202445 &&
01066                             tpcHitTable[ihit].volume_id > 2445) nPseudoPadrow++; 
01067                         else nBadVolId++;
01068                         continue;
01069                     }
01070                     //              g2t_tpc_hitTablePointer->Print(ihit,1);               
01071                     th = new StMcTpcHit(&tpcHitTable[ihit]);
01072                     //              cout << "McTpcHit\t" << ihit << *th << endl;
01073                     if(!mCurrentMcEvent->tpcHitCollection()->addHit(th)) {// adds hit th to collection
01074                         nBadVolId++;
01075                         delete th;
01076                         th = 0;
01077                         continue;
01078                     }
01079                     // point hit to its parent and add it to collection
01080                     // of the appropriate track
01081                     AddHit2Track(tpc,Tpc);
01082                 } // hit loop
01083                 if (Debug()) {
01084                     cout << "Filled " << mCurrentMcEvent->tpcHitCollection()->numberOfHits() << " TPC Hits" << endl;
01085                     cout << "Found " << nPseudoPadrow << " Hits in Pseudo-Padrows." << endl;
01086                     if (nBadVolId) {gMessMgr->Warning() << "StMcEventMaker::Make(): cannot store " << nBadVolId
01087                                                         << " TPC hits, wrong Volume Id." << endm;}
01088                 }
01089                 // Sort the hits
01090                 for (unsigned int iSector=0;
01091                       iSector<mCurrentMcEvent->tpcHitCollection()->numberOfSectors(); iSector++)
01092                     for (unsigned int iPadrow=0;
01093                           iPadrow<mCurrentMcEvent->tpcHitCollection()->sector(iSector)->numberOfPadrows();
01094                           iPadrow++) {
01095                         StSPtrVecMcTpcHit& tpcHits = mCurrentMcEvent->tpcHitCollection()->sector(iSector)->padrow(iPadrow)->hits();
01096                         if (Debug() > 2) {
01097                           Int_t nhits = tpcHits.size();
01098                           cout << "Tpc hits before sort" << endl;
01099                           for (int i = 0; i < nhits; i++) {
01100                             cout << *(tpcHits[i]) << endl;
01101                           }
01102                         }
01103                         sort (tpcHits.begin(), tpcHits.end(), compMcHit() );
01104                         if (Debug() > 2) {
01105                           Int_t nhits = tpcHits.size();
01106                           cout << "Tpc hits after sort" << endl;
01107                           for (int i = 0; i < nhits; i++) {
01108                             cout << *(tpcHits[i]) << endl;
01109                           }
01110                         }
01111                         
01112                     }
01113             } // pointer exists
01114             else {
01115                 if (Debug()) cout << "No TPC Hits in this event" << endl;
01116             }
01117         } // doUseTpc
01118         
01119         //
01120         // SVT Hits
01121         //
01122         if (doUseSvt) {
01123         if (g2t_svt_hitTablePointer) {
01124             StMcSvtHit* th = 0;
01125             long NHits = g2t_svt_hitTablePointer->GetNRows();
01126             long nBadVolId = 0;
01127             long ihit;
01128             for(ihit=0; ihit<NHits; ihit++) {
01129                 if (svtHitTable[ihit].volume_id < 1101 || svtHitTable[ihit].volume_id > 9000) {
01130                     nBadVolId++;
01131                     continue;
01132                 }
01133                 th = new StMcSvtHit(&svtHitTable[ihit]);
01134                 if (!mCurrentMcEvent->svtHitCollection()->addHit(th)) {// adds hit th to collection
01135                     nBadVolId++;
01136                     delete th; // If the hit couldn't be assigned, delete it.
01137                     th = 0;
01138                     continue;
01139                 }
01140 
01141                 // point hit to its parent and add it to collection
01142                 // of the appropriate track
01143                 AddHit2Track(svt,Svt);
01144             }
01145             if (Debug()) {
01146                 cout << "Filled " << mCurrentMcEvent->svtHitCollection()->numberOfHits() << " SVT Hits" << endl;
01147                 cout << "Nhits, " << NHits << " rows from table" << endl;
01148                 if (nBadVolId)
01149                     gMessMgr->Warning() << "StMcEventMaker::Make(): cannot store " << nBadVolId
01150                                         << " SVT hits, wrong Volume Id." << endm;
01151             }
01152             // Sort the hits
01153             for (unsigned int iBarrel=0;
01154                  iBarrel<mCurrentMcEvent->svtHitCollection()->numberOfBarrels(); iBarrel++)
01155                 for (unsigned int iLadder=0;
01156                      iLadder<mCurrentMcEvent->svtHitCollection()->barrel(iBarrel)->numberOfLadders();
01157                      iLadder++)
01158                     for (unsigned int iWafer=0;
01159                          iWafer<mCurrentMcEvent->svtHitCollection()->barrel(iBarrel)->ladder(iLadder)->numberOfWafers();
01160                          iWafer++) {
01161                         StSPtrVecMcSvtHit& svtHits = mCurrentMcEvent->svtHitCollection()->barrel(iBarrel)->ladder(iLadder)->wafer(iWafer)->hits();
01162                         sort (svtHits.begin(), svtHits.end(), compMcHit() );
01163                 
01164             }
01165 
01166         }
01167         else {
01168             if (Debug()) cout << "No SVT Hits in this event" << endl;
01169         }
01170         } // do use svt
01171         
01172         //
01173         // SSD Hits
01174         //
01175         if (doUseSsd) {
01176           if (g2t_ssd_hitTablePointer) {
01177             StMcSsdHit* th = 0;
01178             long NHits = g2t_ssd_hitTablePointer->GetNRows();
01179             long nBadVolId = 0;
01180             long ihit;
01181             for(ihit=0; ihit<NHits; ihit++) {
01182                 if (ssdHitTable[ihit].volume_id < 7000 || ssdHitTable[ihit].volume_id > 9000) {
01183                     nBadVolId++;
01184                     continue;
01185                 }
01186                 th = new StMcSsdHit(&ssdHitTable[ihit]);
01187                 if (!mCurrentMcEvent->ssdHitCollection()->addHit(th)) {// adds hit th to collection
01188                     nBadVolId++;
01189                     delete th; // If the hit couldn't be assigned, delete it.
01190                     th = 0;
01191                     continue;
01192                 }
01193                 
01194                 // point hit to its parent and add it to collection
01195                 // of the appropriate track
01196                 AddHit2Track(ssd,Ssd);
01197             }
01198             if (Debug()) {
01199                 cout << "Filled " << mCurrentMcEvent->ssdHitCollection()->numberOfHits() << " SSD Hits" << endl;
01200                 cout << "Nhits, " << NHits << " rows from table" << endl;
01201                 if (nBadVolId)
01202                     gMessMgr->Warning() << "StMcEventMaker::Make(): cannot store " << nBadVolId
01203                                         << " SSD hits, wrong Volume Id." << endm;
01204             }
01205 
01206         }
01207         else {
01208             if (Debug()) cout << "No SSD Hits in this event" << endl;
01209         }
01210         } // do use ssd
01211 
01212 
01213         //
01214         // FTPC Hits
01215         //
01216         if (doUseFtpc) {
01217         if (g2t_ftp_hitTablePointer) {
01218             StMcFtpcHit* th = 0;
01219             long  NHits = g2t_ftp_hitTablePointer->GetNRows();
01220             long  nBadVolId = 0;
01221             long ihit;
01222             for(ihit=0; ihit<NHits; ihit++) {
01223                 // old volume id scheme, 101 - 209
01224                 // changed in oct 2003
01225                 // planes are included, valid volume id's are from 1000 - 2906
01226                 // keep the checks from 101 to 2906 for backward compatibility.
01227                 if (ftpHitTable[ihit].volume_id < 101 || ftpHitTable[ihit].volume_id > 2906) {
01228                     nBadVolId++;
01229                     continue;
01230                 }
01231 
01232                 th = new StMcFtpcHit(&ftpHitTable[ihit]);
01233 
01234                 if (!mCurrentMcEvent->ftpcHitCollection()->addHit(th)){ // adds hit th to collection
01235                     nBadVolId++;
01236                     delete th;
01237                     th = 0;
01238                     continue;
01239                 }
01240                 // point hit to its parent and add it to collection
01241                 // of the appropriate track
01242                 AddHit2Track(ftp,Ftpc);
01243             }
01244             if (Debug()) {
01245                 cout << "Filled " << mCurrentMcEvent->ftpcHitCollection()->numberOfHits() << " FTPC Hits" << endl;
01246                 if (nBadVolId)
01247                     gMessMgr->Warning() << "StMcEventMaker::Make(): cannot store " << nBadVolId
01248                                         << " FTPC hits, wrong Volume Id." << endm;
01249             }
01250             // Sort the hits
01251             for (unsigned int iPlane=0;
01252                  iPlane<mCurrentMcEvent->ftpcHitCollection()->numberOfPlanes(); iPlane++) {
01253                 StSPtrVecMcFtpcHit& ftpcHits = mCurrentMcEvent->ftpcHitCollection()->plane(iPlane)->hits();
01254                 sort (ftpcHits.begin(), ftpcHits.end(), compMcFtpcHit() );
01255                 
01256             }
01257         }
01258         else {
01259             if (Debug()) cout << "No FTPC Hits in this event" << endl;
01260         }
01261         }// do use ftpc
01262         AddHits(rch,rich,Rich);
01263         AddHits(ctb,ctb,Ctb);
01264         AddHits(tof,tof,Tof);
01265         AddHits(tfr,tof,Tof);
01266         AddHits(mtd,mtd,Mtd);
01267         AddHits(pix,pixel,Pixel);
01268         AddHits(ist,ist,Ist);
01269         AddHits(fgt,fgt,Fgt);
01270         AddHits(etr,etr,Etr);
01271 
01272         // BEMC and BPRS Hits
01273         if (doUseBemc) fillBemc(g2t_emc_hitTablePointer);
01274 
01275         // BSMDE and BSMDP Hits
01276         if (doUseBsmd) fillBsmd(g2t_smd_hitTablePointer);
01277         
01278         // EEMC and EPRS Hits, ESMDU & ESMDV Hits
01279         if (doUseEemc) fillEemc(g2t_eem_hitTablePointer,g2t_esm_hitTablePointer);
01280 
01281         // FPD Hits
01282         if (doUseFpd) fillFpd(g2t_fpd_hitTablePointer);
01283 
01284         // FSC Hits
01285         if (doUseFsc) fillFsc(g2t_fsc_hitTablePointer);
01286 
01287         ttemp.clear();
01288         
01289         //_______________________________________________________________
01290         // At this point StMcEvent should be loaded.
01291                 
01292     }
01293     
01294     
01295     if (!g2t_vertexTablePointer)
01296         gMessMgr->Warning() << "StMcEventMaker:   g2t_vertex  Table Not found " << endm;
01297     if (!g2t_trackTablePointer)
01298         gMessMgr->Warning() << "StMcEventMaker:   g2t_track   Table Not found " << endm;
01299     if (!g2t_tpc_hitTablePointer)
01300         gMessMgr->Info()    << "StMcEventMaker:   g2t_tpc_hit Table Not found " << endm;
01301     if (!particleTablePointer)
01302         gMessMgr->Info()    << "StMcEventMaker:   particle    Table Not found " << endm;
01303     
01304     if (doPrintEventInfo) printEventInfo();
01305     if (doPrintMemoryInfo) {
01306         StMemoryInfo::instance()->snapshot();
01307         StMemoryInfo::instance()->print();
01308     }
01309     if (doPrintCpuInfo) {
01310         timer.stop();
01311         cout << "CPU time for StMcEventMaker::Make(): "
01312              << timer.elapsedTime() << " sec\n" << endl;
01313     }
01314   
01315   return kStOK;
01316 
01317 }
01318 
01319 //_____________________________________________________________________________
01320 void StMcEventMaker::fillBemc(St_g2t_emc_hit* g2t_emc_hitTablePointer)
01321 {
01322     if (g2t_emc_hitTablePointer == 0) {
01323         if (Debug()) cout << "No BEMC and BPRS Hits in this event" << endl;
01324         return;
01325     }
01326 
01327     g2t_emc_hit_st* emcHitTable = g2t_emc_hitTablePointer->GetTable();
01328     StEmcGeom *geomBemc = StEmcGeom::getEmcGeom(1);
01329     int module, eta, sub, detector; 
01330     float de;
01331     StMcTrack *tr; 
01332     StMcCalorimeterHit *emchBemc, *emchBprs;
01333         
01334     StMcEmcHitCollection *bemcColl=mCurrentMcEvent->bemcHitCollection();
01335     StMcEmcHitCollection  *bprsColl=mCurrentMcEvent->bprsHitCollection();
01336     long NHits = g2t_emc_hitTablePointer->GetNRows();
01337                 
01338     for(long ihit=0; ihit<NHits; ihit++,emcHitTable++) { 
01339 
01340         geomBemc->getVolIdBemc(emcHitTable->volume_id, module,eta,sub,detector); // Must check ??
01341 //      cout << "hit       " << ihit << endl;
01342 //      cout << "volume id " << emcHitTable->volume_id << endl;
01343 //      cout << "module    " << module   << endl;
01344 //      cout << "eta       " << eta      << endl;
01345 //      cout << "sub       " << sub      << endl;
01346 //      cout << "detector  " << detector << endl;
01347         tr   = ttemp[emcHitTable->track_p - 1];
01348         de   = emcHitTable->de;
01349                     
01350         if (detector == 1 || detector == 2) {
01351             emchBemc = new StMcCalorimeterHit(module,eta,sub,de, tr);
01352             emchBprs = 0;                                                      // For safety
01353             StMcEmcHitCollection::EAddHit bemcNew = bemcColl->addHit(emchBemc);
01354 
01355             if (bemcNew == StMcEmcHitCollection::kNew){ 
01356                 tr->addBemcHit(emchBemc);
01357                 if(detector == 2) emchBprs = new StMcCalorimeterHit(module,eta,sub,de, tr);
01358             }
01359             else if(bemcNew == StMcEmcHitCollection::kAdd){ 
01360                 emchBprs = emchBemc;
01361             }
01362             else if(bemcNew == StMcEmcHitCollection::kErr){ 
01363                 delete emchBemc;
01364                 emchBemc = 0;
01365                 gMessMgr->Warning()<<"<E> Bad hit in Bemc collection " << endm;
01366             }
01367             else { // Unused branch
01368                 delete emchBemc;
01369                 emchBemc = 0;
01370                 emchBprs = 0;
01371                 gMessMgr->Warning()<<"<E> Funny return value! EAddHit = " << static_cast<int>(bemcNew) <<endm;
01372             }
01373 
01374             if (detector == 2 && emchBprs) {
01375                 StMcEmcHitCollection::EAddHit bprsNew = bprsColl->addHit(emchBprs);
01376                 if (bprsNew == StMcEmcHitCollection::kNew) {
01377                     tr->addBprsHit(emchBprs);
01378                 }
01379                 else { 
01380                     delete emchBprs;
01381                     emchBprs = 0;
01382                 }
01383             }
01384             else if(emchBprs) {
01385                  delete emchBprs;
01386                  emchBprs = 0;
01387             }
01388         }
01389         else {
01390             gMessMgr->Warning() << "<E> Bad EMC detector number " << detector << " volume_id "
01391                                 << emcHitTable->volume_id << " detector "<< detector << endm;
01392         }
01393     }
01394     if (Debug()) {
01395         cout << "Filled " << mCurrentMcEvent->bemcHitCollection()->numberOfHits() << " BEMC Hits" << endl;
01396         cout << "Filled " << mCurrentMcEvent->bprsHitCollection()->numberOfHits() << " BPRS Hits" << endl;
01397     }
01398 }
01399 
01400 //_____________________________________________________________________________
01401 void StMcEventMaker::fillBsmd(St_g2t_emc_hit* g2t_smd_hitTablePointer)
01402 {
01403     if (g2t_smd_hitTablePointer == 0) {
01404         if (Debug()) cout << "No BSMDE and BSMDP Hits in this event" << endl;
01405         return;
01406     }
01407 
01408     g2t_emc_hit_st* smdHitTable = g2t_smd_hitTablePointer->GetTable();
01409     StEmcGeom *geomBsmd = StEmcGeom::getEmcGeom(3);
01410     int module, eta, sub, detector; 
01411     float de;
01412     StMcTrack *tr; 
01413     StMcCalorimeterHit *emchBsmde, *emchBsmdp;
01414 
01415     StMcEmcHitCollection *bsmdeColl=mCurrentMcEvent->bsmdeHitCollection();
01416     StMcEmcHitCollection *bsmdpColl=mCurrentMcEvent->bsmdpHitCollection();
01417     long NHits = g2t_smd_hitTablePointer->GetNRows();
01418     for(long ihit=0; ihit<NHits; ihit++,smdHitTable++) { 
01419         geomBsmd->getVolIdBsmd(smdHitTable->volume_id, module,eta,sub,detector); // Must check ??
01420         tr   = ttemp[smdHitTable->track_p - 1];
01421         de   = smdHitTable->de;
01422 
01423         if (detector == 3) {
01424             emchBsmde = new StMcCalorimeterHit(module,eta,sub,de, tr);
01425             StMcEmcHitCollection::EAddHit bsmdeNew = bsmdeColl->addHit(emchBsmde);
01426             if (bsmdeNew == StMcEmcHitCollection::kNew) {
01427                 tr->addBsmdeHit(emchBsmde);
01428             }
01429             else {
01430                 delete emchBsmde;
01431                 emchBsmde = 0;
01432             } 
01433         }
01434         else if (detector == 4) {
01435             emchBsmdp = new StMcCalorimeterHit(module,eta,sub,de, tr);
01436             StMcEmcHitCollection::EAddHit bsmdpNew = bsmdpColl->addHit(emchBsmdp);
01437             if (bsmdpNew == StMcEmcHitCollection::kNew) {
01438                 tr->addBsmdpHit(emchBsmdp);
01439             }
01440             else {
01441                 delete emchBsmdp;
01442                 emchBsmdp = 0;
01443             } 
01444         }
01445         else {
01446             gMessMgr->Warning() << "<E> Bad SMD detector number " << detector << " volume_id "
01447                                 << smdHitTable->volume_id << " detector "<< detector << endm;
01448         }
01449     }
01450     if (Debug()) {
01451         cout << "Filled " << mCurrentMcEvent->bsmdeHitCollection()->numberOfHits() << " BSMDE Hits" << endl;
01452         cout << "Filled " << mCurrentMcEvent->bsmdpHitCollection()->numberOfHits() << " BSMDP Hits" << endl;
01453     }
01454 }
01455 
01456 //_____________________________________________________________________________
01457 void StMcEventMaker::fillEemc(St_g2t_emc_hit* g2t_tile, St_g2t_emc_hit* g2t_smd){
01458     if (Debug()) {
01459        gMessMgr->Info() << GetName() << "fillEemc() called" << endm;
01460     }
01461     EEmcMCData mEemcGeant;
01462     mEemcGeant.unpackGeantHits(g2t_tile, g2t_smd);
01463     if (Debug()>1) {
01464         mEemcGeant.print();
01465     }
01466     
01467     StMcEmcHitCollection *eemcColl=mCurrentMcEvent->eemcHitCollection();
01468     StMcEmcHitCollection *eprsColl=mCurrentMcEvent->eprsHitCollection();
01469     StMcEmcHitCollection *esmduColl=mCurrentMcEvent->esmduHitCollection();
01470     StMcEmcHitCollection *esmdvColl=mCurrentMcEvent->esmdvHitCollection();
01471     
01472     if (!eemcColl || !eprsColl || !esmduColl || !esmdvColl) {
01473         gMessMgr->Warning()<<GetName() <<"::fillEemc(), sth wrong with StMcEEmcCollection,\n    skip EEMC GEANT hits"<<endm;
01474         return;
01475     }
01476     
01477     
01478     int nHit;
01479     const EEmcMCHit *h  = mEemcGeant.getGeantHits(nHit);
01480     for(Int_t i=0; i<nHit; i++,h++) {
01481         int detId=h->detector;
01482         assert(h->track_p>0); // tmp, to catch bugs,JB
01483         StMcTrack *tr = ttemp[h->track_p - 1];
01484         int Beta=0,Bsub=0,Bmodule=0; // barrel indexes
01485         /* barrel indexes are used to lable eemc hits 
01486            B-names are misleading but preserved for consistency
01487            The mapping must follow the exactly the scheme as at
01488            /StRoot/StEEmcSimulatorMaker/StEEmcFastMaker::mEE2ST()
01489            Otherwise internal indexing in EmcCollection can overwrite
01490            adrresses and/or association will be broken. Jan Balewski
01491         */
01492         Bmodule=h->sector;
01493         int off=0; // used only for pres1,pres2, and postshower
01494         StMcCalorimeterHit * stMcHit=0;
01495         StMcEmcHitCollection::EAddHit addRet=StMcEmcHitCollection::kErr;
01496         if (Debug()>1) printf("StMcAdd:detID=%d iHit=%d de=%g\n",detId,i,h->de);
01497         
01498         switch(detId) { // covers all 6 layers of EEMCS
01499             
01500         case EEmcMCData::kEEmcMCTowerId:
01501             Beta=h->tower.eta;
01502             Bsub=h->tower.ssec;
01503             stMcHit = new StMcCalorimeterHit(Bmodule,Beta,Bsub,h->de,tr);
01504             addRet = eemcColl->addHit(stMcHit);
01505             if(addRet==StMcEmcHitCollection::kNew) tr->addEemcHit(stMcHit);
01506             break;
01507             
01508       // pre/post hits are stored in the same container
01509         case EEmcMCData::kEEmcMCPostShowerId: off++;
01510         case EEmcMCData::kEEmcMCPreShower2Id: off++;
01511         case EEmcMCData::kEEmcMCPreShower1Id:
01512             Beta=h->tower.eta;
01513             Bsub=h->tower.ssec +5*off;
01514             stMcHit = new StMcCalorimeterHit(Bmodule,Beta,Bsub,h->de,tr);
01515             addRet = eprsColl->addHit(stMcHit);
01516             if(addRet==StMcEmcHitCollection::kNew) tr->addEprsHit(stMcHit);
01517             break;
01518 
01519         case EEmcMCData::kEEmcMCSmdUStripId:
01520             Beta=h->strip;
01521             Bsub=1;
01522             stMcHit = new StMcCalorimeterHit(Bmodule,Beta,Bsub,h->de,tr);
01523             addRet  = esmduColl->addHit(stMcHit);
01524             if(addRet==StMcEmcHitCollection::kNew) tr->addEsmduHit(stMcHit);
01525             break;
01526 
01527         case EEmcMCData::kEEmcMCSmdVStripId:
01528             Beta=h->strip;
01529             Bsub=1; 
01530             stMcHit = new StMcCalorimeterHit(Bmodule,Beta,Bsub,h->de,tr);      
01531             addRet  = esmdvColl->addHit(stMcHit);
01532             if(addRet==StMcEmcHitCollection::kNew) tr->addEsmdvHit(stMcHit); 
01533             break;
01534         default:
01535             gMessMgr->Warning()<<GetName() <<"::fillEemc(), wrong detectorId=" << detId << " from EEmcMCHit. iHit="<<i<<"\n It is fatal - bug in the code, fix it, JB"<<endm;
01536             assert(1==2);
01537             break;
01538         }
01539         
01540         // garbage  collector
01541         switch (addRet ){
01542         case StMcEmcHitCollection::kNew: break;
01543         case StMcEmcHitCollection::kAdd: delete stMcHit; break; // delete used but not needed anymore hit
01544         case StMcEmcHitCollection::kErr:
01545             gMessMgr->Warning()<<"<E> Bad hit in Eemc collection" << endm;
01546             delete stMcHit;
01547             break;
01548         default:
01549             gMessMgr->Warning()<<"<E> Funny return value! EAddHit = " << static_cast<int>(addRet) << endm;
01550             delete stMcHit;
01551             break;
01552         }
01553     }// end of loop over GEANT hits
01554     if (Debug()) gMessMgr->Info()<<GetName() <<"::fillEemc()   done, nHit="<<nHit<<endm;
01555 }
01556 
01557 //_____________________________________________________________________________
01558 void StMcEventMaker::fillFpd(St_g2t_emc_hit* g2t_fpd_hitTablePointer)
01559 {
01560   if (!g2t_fpd_hitTablePointer) {
01561     if (Debug()) cout << "No FPD Hits in this event" << endl;
01562     return;
01563   }
01564 
01565   g2t_emc_hit_st* fpdHitTable = g2t_fpd_hitTablePointer->GetTable();
01566   StMcEmcHitCollection* fpdColl = mCurrentMcEvent->fpdHitCollection();
01567 
01568   for (int iHit = 0; iHit < g2t_fpd_hitTablePointer->GetNRows(); ++iHit) {
01569     g2t_emc_hit_st& hit = fpdHitTable[iHit];
01570     // volume_id = ew*10000+nstb*1000+ch
01571     int volume_id = hit.volume_id;
01572     int ch = volume_id % 1000;
01573     volume_id /= 1000;
01574     int nstb = volume_id % 10;
01575     int ew = volume_id / 10;
01576     StMcTrack* track = ttemp[hit.track_p - 1];
01577     // Store ew in module, nstb in sub, and ch in eta
01578     StMcCalorimeterHit* fpdHit = new StMcCalorimeterHit(ew,ch,nstb,hit.de,track);
01579     StMcEmcHitCollection::EAddHit fpdNew = fpdColl->addHit(fpdHit);
01580     if (fpdNew == StMcEmcHitCollection::kNew) {
01581       track->addFpdHit(fpdHit);
01582     } else {
01583       delete fpdHit; fpdHit = 0;
01584     }
01585   }
01586 
01587   if (Debug()) cout << "Filled " << mCurrentMcEvent->fpdHitCollection()->numberOfHits() << " FPD Hits" << endl;
01588 }
01589 
01590 //_____________________________________________________________________________
01591 void StMcEventMaker::fillFsc(St_g2t_emc_hit* g2t_fsc_hitTablePointer)
01592 {
01593   if (!g2t_fsc_hitTablePointer) {
01594     if (Debug()) cout << "No FSC Hits in this event" << endl;
01595     return;
01596   }
01597 
01598   g2t_emc_hit_st* fscHitTable = g2t_fsc_hitTablePointer->GetTable();
01599   StMcEmcHitCollection* fscColl = mCurrentMcEvent->fscHitCollection();
01600 
01601   for (int iHit = 0; iHit < g2t_fsc_hitTablePointer->GetNRows(); ++iHit) {
01602     g2t_emc_hit_st& hit = fscHitTable[iHit];
01603 //                      volume_id = ew*10000+nstb*1000+ch
01604     int volume_id = hit.volume_id;
01605 //      std::cout << "FSC vol_id: " << volume_id << "\n";
01606     int module = 1;             // only one FSC exists at the moment :)
01607     int eta = floor(float(volume_id) / float(80)); // X coordinate
01608     int sub = volume_id % 80;   // Y coordinate
01609     StMcTrack* track = ttemp[hit.track_p - 1];
01610     // Store ew in module, nstb in sub, and ch in eta
01611     StMcCalorimeterHit* fscHit = new StMcCalorimeterHit(module,eta,sub,hit.de,track);
01612     StMcEmcHitCollection::EAddHit fscNew = fscColl->addHit(fscHit);
01613     if (fscNew == StMcEmcHitCollection::kNew) {
01614       track->addFscHit(fscHit);
01615     } else {
01616       delete fscHit; fscHit = 0;
01617     }
01618   }
01619 
01620   if (Debug()) cout << "Filled " << mCurrentMcEvent->fscHitCollection()->numberOfHits() << " FSC Hits" << endl;
01621 }
01622 
01623 
01624 //_____________________________________________________________________________
01625 void StMcEventMaker::printEventInfo()
01626 {
01627     cout << "---------------------------------------------------------" << endl;
01628     cout << "StMcEvent at " << (void*) mCurrentMcEvent                  << endl;
01629     cout << "---------------------------------------------------------" << endl;
01630     if (! mCurrentMcEvent) return;
01631     cout << *mCurrentMcEvent << endl;
01632     mCurrentMcEvent->Print("");
01633 }