StRoot: St_tcl_Maker/StTpcFastSimMaker.cxx Source File

00001 /* $Id: StTpcFastSimMaker.cxx,v 1.5 2012/05/07 14:54:45 fisyak Exp $
00002     $Log: StTpcFastSimMaker.cxx,v $
00003     Revision 1.5  2012/05/07 14:54:45  fisyak
00004     Add printout
00005 
00006     Revision 1.4  2011/01/04 21:40:22  fisyak
00007     Add pile-up
00008 
00009     Revision 1.3  2010/08/16 21:59:46  fisyak
00010     leave coordinates in TpcLocalCoordinate because StTpcHitMover expects that
00011 
00012     Revision 1.2  2010/05/28 16:28:44  fisyak
00013     Adjust for new TpcDb interface, remove pseudo pad rows
00014 
00015     Revision 1.1  2009/11/10 21:15:33  fisyak
00016     pams clean up
00017  
00018 */
00019 #include "Stiostream.h"
00020 #include "StTpcFastSimMaker.h"
00021 #include "StHit.h"
00022 #include "StEventTypes.h"
00023 #include "tables/St_HitError_Table.h"
00024 #include "StTpcHitCollection.h"
00025 #include "StEvent.h"
00026 #include "TMath.h"
00027 #include "TRandom.h"
00028 #include "StTpcDb/StTpcDb.h"
00029 #include "tables/St_g2t_vertex_Table.h"
00030 #include "tables/St_g2t_track_Table.h"
00031 #include "tables/St_g2t_tpc_hit_Table.h"
00032 #include "StDbUtilities/StTpcCoordinateTransform.hh"
00033 #include "StDbUtilities/StMagUtilities.h"
00034 #include "StTpcDb/StTpcDb.h"
00035 #include "TDataSet.h"
00036 #include "TDataSetIter.h"
00037 #include "StDetectorDbMaker/StiTpcInnerHitErrorCalculator.h"
00038 #include "StDetectorDbMaker/StiTpcOuterHitErrorCalculator.h"
00039 ClassImp(StTpcFastSimMaker);
00040 //____________________________________________________________
00041 Int_t StTpcFastSimMaker::Make() {
00042   mExB = gStTpcDb->ExB();
00043   if (! gRandom) gRandom = new TRandom();
00044   // Get the input data structures from StEvent
00045   StEvent *rEvent =  (StEvent*) GetInputDS("StEvent");
00046   if (! rEvent) {    LOG_WARN << "No StEvent on input, bye bye" << endm; return kStWarn; }
00047   StTpcHitCollection *rCol = rEvent->tpcHitCollection();
00048   if (!rCol) {
00049     rCol = new StTpcHitCollection; 
00050     rEvent->setTpcHitCollection(rCol); 
00051   }
00052   St_g2t_tpc_hit *g2t_tpc_hit = (St_g2t_tpc_hit *) GetDataSet("geant/g2t_tpc_hit");
00053   if (! g2t_tpc_hit) {
00054     LOG_WARN << "No g2t_tpc_hit on input, bye bye" << endm; return kStWarn;
00055    }
00056   Int_t Nhits = g2t_tpc_hit->GetNRows();
00057   if (Nhits <= 0) return kStWarn;
00058   St_g2t_track *g2t_track = (St_g2t_track *) GetDataSet("geant/g2t_track"); //  if (!g2t_track)    return kStWarn;
00059   g2t_track_st *tpc_track = 0;
00060   if (g2t_track) tpc_track = g2t_track->GetTable();
00061   St_g2t_vertex  *g2t_ver = (St_g2t_vertex *) GetDataSet("geant/g2t_vertex");// if (!g2t_ver)      return kStWarn;
00062   g2t_vertex_st     *gver = 0;
00063   if (g2t_ver) gver = g2t_ver->GetTable();
00064   g2t_tpc_hit_st *tpc_hit         = g2t_tpc_hit->GetTable();
00065   StTpcCoordinateTransform transform(gStTpcDb);
00066   for (Int_t i = 0; i < Nhits; i++)    {
00067     if (tpc_hit[i].volume_id > 100000) continue; // skip pseudo pad rows
00068     Int_t Id         = tpc_hit[i].track_p;
00069     Int_t id3 = 0;
00070     if (tpc_track) 
00071       id3        = tpc_track[Id-1].start_vertex_p;
00072     Int_t sector = (tpc_hit[i].volume_id%10000)/100;
00073     Int_t row    =  tpc_hit[i].volume_id%100;
00074     StGlobalDirection     dirG(tpc_hit[i].p[0],tpc_hit[i].p[1],tpc_hit[i].p[2]);
00075     static StTpcLocalSectorDirection dirL;
00076     transform(dirG, dirL, sector, row);
00077     StGlobalCoordinate    coorG(tpc_hit[i].x[0],tpc_hit[i].x[1],tpc_hit[i].x[2]);
00078     static StTpcLocalCoordinate  coorLT;
00079     transform(coorG,coorLT,sector,row);
00080     StTpcLocalCoordinate  coorLTD = coorLT;
00081     // ExB corrections
00082     Float_t pos[3] = {coorLTD.position().x(),coorLTD.position().y(),coorLTD.position().z()};
00083     Float_t posMoved[3];
00084     if ( mExB ) {
00085       mExB->DoDistortion(pos,posMoved);   // input pos[], returns posMoved[]
00086       StThreeVector<double> newPos(posMoved[0],posMoved[1],posMoved[2]);
00087       coorLTD.setPosition(newPos);
00088     }
00089     static StTpcLocalSectorAlignedCoordinate  coorLSA;
00090     transform(coorLTD,coorLSA); //
00091     static StTpcLocalSectorCoordinate  coorLS;
00092     transform(coorLSA,coorLS); // alignment
00093     Double_t xyzL[3] = {coorLS.position().x(),coorLS.position().y(),coorLS.position().z()};
00094     if (Debug() && TMath::Abs(xyzL[1]-transform.yFromRow(row)) > 0.1000) {
00095       cout << "Id: " << tpc_hit[i].volume_id  
00096            << "\txyzL :" << xyzL[0] << "\t" << xyzL[1] << "\t" << xyzL[2] 
00097            << "\tdR :" << xyzL[1]-transform.yFromRow(row) << endl;
00098     }
00099     Double_t Z = xyzL[2];
00100     Double_t eta = TMath::PiOver2() - TMath::Abs(dirL.position().phi());
00101     Double_t tanl = dirL.position().z()/dirL.position().perp();
00102     Double_t sigmaY2, sigmaZ2;
00103     if (row <= 13)  StiTpcInnerHitErrorCalculator::instance()->calculateError(Z,eta,tanl,sigmaY2, sigmaZ2);
00104     else            StiTpcOuterHitErrorCalculator::instance()->calculateError(Z,eta,tanl,sigmaY2, sigmaZ2);
00105     Double_t sigmaY = TMath::Sqrt(sigmaY2);
00106     Double_t sigmaZ = TMath::Sqrt(sigmaZ2);
00107     StThreeVectorF e(0, sigmaY, sigmaZ);
00108     xyzL[0] += gRandom->Gaus(0, sigmaY);
00109     xyzL[2] += gRandom->Gaus(0, sigmaZ);
00110     StThreeVectorD newPosition(xyzL);
00111     coorLS.setPosition(newPosition);
00112     static StTpcPadCoordinate Pad;
00113     transform(coorLS,Pad,kFALSE,kTRUE); // don't use T0, use Tau
00114     Double_t tof = 0;
00115     if (gver) tof = gver[id3-1].ge_tof;
00116     tof += tpc_hit[i].tof;
00117     Float_t  timebkt =  Pad.timeBucket() + 1.e6*gStTpcDb->Electronics()->samplingFrequency()*tof;
00118     if (timebkt < 0 || timebkt > 512) continue;
00119     StTpcPadCoordinate newPad(Pad.sector(),Pad.row(), Pad.pad(),timebkt );
00120     static StTpcLocalCoordinate global; // leave coordinates in TpcLocalCoordinate because StTpcHitMover expects that.
00121     transform(newPad,global,kFALSE); // alignment
00122     StThreeVectorF p(global.position().x(),global.position().y(),global.position().z());
00123     UInt_t hw = 1;   // detid_tpc
00124     hw += sector << 4;     // (row/100 << 4);   // sector
00125     hw += row    << 9;     // (row%100 << 9);   // row
00126     StTpcHit *tpcHit = new StTpcHit(p,e, 
00127                                     hw,TMath::Abs(tpc_hit[i].de), i+1,   // hw, q, c
00128                                     tpc_hit[i].track_p, 100,             // idTruth, quality
00129                                     0,                                   // id
00130                                     0,  0, 0,                            // mnpad, mxpad, mntmbk
00131                                     0, Pad.pad(), Pad.timeBucket());     // mxtmbk, cl_x , cl_t
00132     if (Debug() > 1) tpcHit->Print();
00133     rCol->addHit(tpcHit);
00134   }
00135   return kStOK;
00136 }