StRoot: StEEmcPool/FgtSandbox1/StFgtSlowSimuMaker.cxx Source File

00001 // *-- Author : J.Balewski
00002 // 
00003 // $Id: StFgtSlowSimuMaker.cxx,v 1.5 2011/04/11 19:35:38 fisyak Exp $
00004 #include <TVector3.h>
00005 #include <TH2.h>
00006 #include <TF1.h>
00007 #include <TFile.h>
00008 #include <TLine.h>
00009 #include <TPolyLine.h>
00010 #include <TCrown.h>
00011 #include <TRandom3.h>
00012 
00013 #include "StFgtSlowSimuMaker.h"
00014 #include "HexLatice.h"
00015 
00016 #include  "tables/St_g2t_fgt_hit_Table.h"
00017 
00018 ClassImp(StFgtSlowSimuMaker)
00019 
00020 //--------------------------------------------
00021 StFgtSlowSimuMaker::StFgtSlowSimuMaker(const char *name):StMaker(name){
00023   setHList(0);
00024   memset(hA,0,sizeof(hA));
00025   geom=new StFgtGeom();
00026   mRnd = new TRandom3(); // general use random generator
00027   //  mRnd->SetSeed(0); // activate, assure every set of data is different
00028 
00029   forcePerpTracks(false); // normal operation
00030   meLossTab[0]=-999; // default=uninitialized
00031 
00032   par_2DpixAmplThres=0.1; // a.u. MIP respons has maxAmpl~200-300== # of electrons
00033   par_stripAmplThres=1.0; // a.u., drop strips below it, default=1.0
00034   par_XYamplSigma=0.035 ; // cm, for 2D gauss smearing, default=0.035
00035 
00036   par_radStripGainMean=1.0;par_radStripGainSigma=0.0;// allow unequal energy sharing
00037   par_cutoffOfBichel = 9992; // meLossTab[9993-10000]=5.0E+03 to 6.0E+06 (too high) are cut off. They are replaced by meLossTab[9992]=4.78E+03 in simulation. 
00038   mRadStripRelativeGain=0;// clear pointer
00039   mPhiStripRelativeGain=0;// clear pointer
00040   hexLat=0;
00041 }
00042 
00043 //--------------------------------------------
00044 void 
00045 StFgtSlowSimuMaker::Clear(Option_t *) { 
00046   int i,j;
00047   for(i=0;i<kFgtMxDisk;i++){
00048     for(j=0;j<kFgtMxQuad;j++) {
00049       mG2tHitList[i][j].clear();
00050     }
00051     mRadAdcList[i].clear();
00052     mPhiAdcList[i].clear();
00053   }
00054   StMaker::Clear();
00055 }
00056 
00057 
00058 //--------------------------------------------
00059 StFgtSlowSimuMaker::~StFgtSlowSimuMaker(){
00060 
00061 }
00062 
00063 //_______________________________________________
00064 //________________________________________________
00065 void
00066 StFgtSlowSimuMaker::saveHisto(TString fname){
00067 
00068   CloseHisto();
00069 
00070   TString outName=fname+".hist.root";
00071   TFile f( outName,"recreate");
00072   assert(f.IsOpen());
00073   printf("%d histos are written  to '%s' ...\n",HList->GetEntries(),outName.Data());
00074 
00075   HList->Write();
00076   f.Close();
00077 
00078 }
00079 
00080 
00081 //--------------------------------------------
00082 //--------------------------------------------
00083 //--------------------------------------------
00084 Int_t 
00085 StFgtSlowSimuMaker::Finish(){
00086   LOG_INFO<<"::Finish() \n"<<  endm; 
00087 
00088 
00089   return StMaker::Finish();
00090 }
00091 
00092  
00093 
00094 //--------------------------------------------
00095 //--------------------------------------------
00096 //--------------------------------------------
00097 Int_t 
00098 StFgtSlowSimuMaker::Init(){
00099   LOG_INFO<<"::Init() "<<  endm; 
00100   assert(HList);
00101   assert(  meLossTab[0]>0); //Frank's model of # of electron pairs must be initialized 
00102   mInpEve=0;
00103 
00104   /* model transverse response in 1D
00105      Frank: However, from test beam data I calculate the RMS of each cluster, 
00106      and there I get something like 0.5 - 0.6 strips, which is around 350 um.
00107      all clusters (both in my diploma thesis and in my email) are 1D. A
00108      cluster is defined as something on ONE projection of the detector. So
00109      the FWHM and the sigma is all 1D. So you should not divide by Sqrt(2) to
00110      get the 1 D value.
00111   */
00112 
00113   amplF= new TF1("Gs","exp( -(x-[0])*(x-[0])/2./[1]/[1])",-1.,1.); // X in cm
00114   amplF->SetParNames("X0","sig");
00115   amplF ->SetParameters(0.,par_XYamplSigma);
00116   amplF ->SetLineWidth(2); amplF->SetLineColor(kMagenta);
00117 
00118   InitHisto1();
00119   InitHisto2();
00120 
00121   // initialize RAD- gains
00122   mRadStripRelativeGain=new double[geom->radStripGBLId_number()+1]; // allocate array
00123   memset(mRadStripRelativeGain,0,sizeof(mRadStripRelativeGain));
00124   int i;
00125   for(i=0;i<geom->radStripGBLId_number();i++) {
00126     double del=99999;
00127     while(fabs(del)> 2.*par_radStripGainSigma)
00128       del= mRnd->Gaus(0,par_radStripGainSigma);
00129       mRadStripRelativeGain[i]=par_radStripGainMean+del;
00130       // if(i<20) printf("rad%d %.4f\n",i,mRadStripRelativeGain[i]);
00131   }
00132   
00133   // initialize PHI- gains
00134   mPhiStripRelativeGain=new double[geom->phiStripGBLId_number()+1]; // allocate array
00135   memset(mPhiStripRelativeGain,0,sizeof(mPhiStripRelativeGain));
00136   for(i=0;i<geom->phiStripGBLId_number();i++) {
00137     double del=99999;
00138     while(fabs(del)> 2.*par_phiStripGainSigma)
00139       del= mRnd->Gaus(0,par_phiStripGainSigma);
00140       mPhiStripRelativeGain[i]=par_phiStripGainMean+del;
00141       // if(i<20) printf("phi%d %.4f\n",i,mPhiStripRelativeGain[i]);
00142   }
00143   
00144   // initialize hexagonal Gem latice
00145   if(par_hexLaticePitch >0.001) {
00146     hexLat=new  HexLatice (par_hexLaticePitch, par_hexLaticePhi1deg);
00147   } else {
00148     LOG_INFO<<Form("::Init hexGemLatice DISABLED, too small pitch")<<endm;
00149     par_hexLaticePitch = 0.0;
00150   }
00151 
00152   geom->printParam();
00153 
00154   LOG_INFO<<Form("::Init params: X,YamplSigma=%.4f cm,  2DpixAmplThres=%.2f a.u., stripAmplThres=%.2f a.u.,  forcePerpTracks=%d  useOnlyDisk=%d RadGain(m=%.2f,sig=%.2f)  PhiGain(m=%.2f,sig=%.2f, GemHexLatice: pitch/um=%.1f, phi1/deg=%.1f, transDiffusion=%.1f um/1cmOfPath, cutoffOfBichel=%d)",par_XYamplSigma,par_2DpixAmplThres,par_stripAmplThres,par_forcePerp,par_useOnlyDisk,par_radStripGainMean,par_radStripGainSigma,par_phiStripGainMean,par_phiStripGainSigma,par_hexLaticePitch*1e4, par_hexLaticePhi1deg,par_transDiffusionPerPath*1e4,par_cutoffOfBichel)<<  endm; // fix it
00155 
00156   assert(par_XYamplSigma>0);
00157 
00158   Info("Init","testing access to TGeoManager");
00159 
00160   if (gGeoManager) { // Geom already there
00161     Info("Load","TGeoManager(%s,%s) is already there"
00162             ,gGeoManager->GetName(),gGeoManager->GetTitle());
00163   } else {
00164     Warning("Init","add  TGeoManager");
00165     //   TString geo="y2009a"; // got it from Victor
00166     //   TString ts("$STAR/StarDb/VmcGeometry/");
00167     //  ts+=geo; ts+=".h"; 
00168     TString ts = "SandBox.C(0,0)";
00169     //printf("WILL execute macro=%s=\n",ts.Data()); 
00170     int ierr=0;
00171     gROOT->Macro(ts, &ierr);
00172     assert(!ierr);
00173   }
00174   assert(gGeoManager);
00175 
00176   TGeoVolume *geoFgt = gGeoManager ->FindVolumeFast("FGMO");
00177   assert(geoFgt);
00178   geoFgt -> Print();
00179   geoFgt -> InspectMaterial();
00180   geoFgt -> InspectShape();
00181   printf("XXXXXXXXXXXXXXXXXXXXXXXXX\nXXXXXXXXXXXXX\nXXXXXXXXXX\n");
00182 
00183 
00184 
00185   return StMaker::Init();
00186 }
00187 
00188 //--------------------------------------------
00189 //--------------------------------------------
00190 //--------------------------------------------
00191 Int_t 
00192 StFgtSlowSimuMaker::Make(){
00193   mInpEve++;
00194   
00195   LOG_INFO <<"::Make() inpEve="<<  mInpEve<<endm;
00196   //============ FGT ==========
00197   St_g2t_fgt_hit *fgt_hitT = (St_g2t_fgt_hit *) GetDataSet("g2t_fgt_hit");
00198   if(fgt_hitT ==0) {
00199     LOG_FATAL<<Form("g2t_Fgt table empty")<<endm;
00200     return kStOK;
00201   }
00202 
00203   /* In year 2012 only 6 disks are build and installed in STAR
00204      This implementation allows for up to 8 disk to be used in reco and
00205      additional disk #9 can exist in GSTR but sort_g2t_hits()  will dropp them.
00206      If disk #10 or more is reported the code will crash
00207   */
00208 
00209   sort_g2t_hits( fgt_hitT );
00210 
00211   int iDisk=0,iQuad=0;
00212   for(iDisk=0;iDisk<kFgtMxDisk;iDisk++) { 
00213     digRad->Reset();  digPhi->Reset();  //clear once per disk  
00214     if(par_useOnlyDisk) if(par_useOnlyDisk-1 != iDisk) continue;
00215     for(iQuad=0;iQuad<kFgtMxQuad;iQuad++) {
00216       digXY->Reset(); // needs reset for every quadrant
00217       
00218       printf(" process %d g2t hits in disk=%d quad=%d\n", mG2tHitList[iDisk][iQuad].size(),iDisk,iQuad);
00219       // accumulate response in digXY array
00220       vector<fgt_g2t_auxil> &L=mG2tHitList[iDisk][iQuad];
00221       if(L.size()<=0) continue;// drop if empty quad
00222       for(UInt_t
00223  i=0;i<L.size();i++) { // populate: digXY
00224         //responseLinearModel(h->Rloc,h->Dloc);
00225         responseFrankModel(L[i].Rloc,L[i].Dloc);
00226         hA[11+iDisk]->Fill(L[i].Rlab.x(),L[i].Rlab.y()); // monitor hit distribution
00227         printf("iQ=%d itr=%d  R/cm=%.2f  phi/deg=%.3f\n",iQuad,i,L[i].Rlab.Perp(), L[i].Rlab.Phi()/3.1416*180); // pi value corrected WMZ
00228       }
00229       // now response is stored in digXY array
00230       if(!projectQuadrant(iQuad)) continue;// drop if empty quad
00231       // strip response is generated and projected to 1D stip histos: R & Phi
00232       // break; // only 1st quad, tmp
00233     }
00234     // only once per disk
00235     exportStripPlane(digRad,mRadAdcList[iDisk]);
00236     exportStripPlane(digPhi,mPhiAdcList[iDisk]);
00237     // break; 
00238   } // end of disk
00239 
00240   printf(" Summary of fired strips\n disk  #Rad-strip  #Phi-strip \n");
00241   for(iDisk=0;iDisk<kFgtMxDisk;iDisk++) printf("%d  %d  %d \n",iDisk,mRadAdcList[iDisk].size(),mPhiAdcList[iDisk].size());
00242   
00243   return kStOK;
00244 }
00245 
00246 //--------------------------------------------
00247 //--------------------------------------------
00248 //--------------------------------------------
00249 void 
00250 StFgtSlowSimuMaker::sort_g2t_hits( St_g2t_fgt_hit *fgt_hitT){
00251   g2t_fgt_hit_st *hitPtr  = fgt_hitT->GetTable();
00252   assert(hitPtr);
00253   LOG_INFO<<Form("Sorting g2t FGT hits, size=%d",fgt_hitT->GetNRows())<<endm;
00254 
00255   int ntot=0;
00256   Int_t     nhits      = fgt_hitT->GetNRows();
00257   for(Int_t ih=0; ih<nhits; ih++,hitPtr++) {
00258     Int_t   ivid  = hitPtr->volume_id;
00259 
00260     // decode volume_id WMZ
00261     Int_t numbv1 = ivid/1000000;  
00262     Int_t numbv2 = (ivid/10000)%100;
00263     int diskID, iQuad;  // diskID and quadID, both start from 0
00264 
00265 //No Quad ID for now, it assigned to 0 for now
00266      diskID = numbv1 - 1;
00267      iQuad = 0; //temp
00268 //    
00269 /* Disk9 is not a concer for now, commented out 
00270     if(numbv2 != 0) {
00271        diskID = numbv1 - 1;
00272        iQuad = numbv2 - 1;
00273     } else {
00274        diskID = 8;
00275        iQuad = numbv1 - 1;
00276     }
00277 */
00278     cout << " Volume_id diskID QuadID: " 
00279          << ivid << " " << diskID << " " << iQuad << endl;
00280     
00281     TVector3 Rlab( hitPtr->x); // entrance point  in Lab ref 
00282     TVector3 Rloc=Rlab; 
00283 /*
00284  IQuad my not match Rlab.Phi() for hits at edges between two quadrants 
00285  (deadQuadEdge)! The following conversion would rotate Rlab with a wrong 
00286  angle if a mismatch happens. However, Rloc.x() or Rloc.y() would be still 
00287  very close to one of two edges of a local quadrant and the cuts on Rloc.x() 
00288  and Rloc.y() later would reject them.    WMZ
00289 */
00290 
00291     Rloc.RotateZ(-geom->phiQuadXaxis(iQuad));
00292 
00293     hA[3]->Fill(1);
00294     if(diskID!=8 && !geom->inDisk(Rlab)) continue;  // drop 8 ? JB
00295     hA[3]->Fill(2);
00296 
00297     if(fabs(Rloc.x()) < geom->deadQuadEdge())  continue;
00298     hA[3]->Fill(3);
00299     if(fabs(Rloc.y()) < geom->deadQuadEdge())  continue;
00300     hA[3]->Fill(4);
00301 
00302     double tof_ns=hitPtr->tof*1.e9;
00303     hA[5]->Fill(tof_ns);
00304     if(hitPtr->tof> geom-> maxTof() ) continue;
00305     hA[3]->Fill(5);
00306 
00307     TVector3 Plab(hitPtr->p);
00308     if(Plab.Mag()>0) hA[7]->Fill(log10(Plab.Mag()*1000.));
00309     if(Plab.Mag()< geom-> minPmag() ) continue;
00310     hA[3]->Fill(6);
00311 
00312 
00313     /*
00314       TVector3 RGeant=Rloc;    RGeant.RotateZ(-0.7854); 
00315       
00316       cout << " R = " <<  Rlab.Perp() << endl;
00317       cout << " Rlab_x, _y   = " << Rlab.x() << " " << Rlab.y() << endl;
00318       cout << " Rloc_x, _y   = " << Rloc.x() << " " << Rloc.y() << endl;
00319       cout << " RGeant_x, _y = " << RGeant.x() << " " << RGeant.y() << endl;
00320       cout << " phiG, phiLoc, phiLab  = " << RGeant.Phi()*180/3.14159 << " "
00321       << Rloc.Phi()*180/3.14159 << " "
00322       << Rlab.Phi()*180/3.14159 << " " << endl;
00323     */
00324     
00325     //  Compare quadID from StFgtGeom and the one decoded from Geant  WMZ
00326     int iQuadChk=geom->getQuad(Rlab.Phi());
00327     if(iQuad != iQuadChk) {
00328       cout << "  Something wrong in quadIDs!! " << endl;
00329       cout << "  iQuad from StFgtGeom and iQuad from Geant  = " 
00330            << iQuadChk << ", " << iQuad << endl;
00331     }
00332     
00333     //  hA[3]->Fill(7...9); free
00334     
00335     //....... hit is accepted for processing ...............
00336     double ds=hitPtr->ds;
00337    
00338     TVector3 verLab=Plab.Unit(); // direction of track in LAB
00339     if(par_forcePerp) verLab=TVector3(0,0,1);// for testing ONLY:  make track perp to GEM
00340     TVector3 Rloc2=Rlab+ds*verLab;   Rloc2.RotateZ(-geom->phiQuadXaxis(iQuad));
00341     TVector3 Dloc=Rloc2-Rloc; // local vector along the path
00342 
00343     fgt_g2t_auxil aux;
00344     aux.hitPtr=hitPtr;
00345     aux.Rloc=Rloc;
00346     aux.Dloc=Dloc;
00347     aux.Rlab=Rlab+ds/2.*verLab; // set it in the middle of the track segment
00348 
00349     aux.iQuad=iQuad;
00350     mG2tHitList[diskID][iQuad].push_back(aux); // ID starts from 0, WMZ
00351 
00352     ntot++;
00353 
00354     //...... only QA is below .......
00355 
00356     // QA  accpted hits
00357     hA[3]->Fill(10+5*diskID);
00358     hA[3]->Fill(10+5*diskID+iQuad+1);
00359     double de_kev=hitPtr->de*1.e6;
00360     hA[0]->Fill(de_kev);
00361     hA[1]->Fill(hitPtr->ds);
00362     hA[2]->Fill(Rlab.z());
00363     hA[4]->Fill(Rlab.x(),Rlab.y());
00364     hA[6]->Fill(Rlab.Perp());
00365     
00366   }// loop over hits
00367   LOG_INFO<<Form("Sorting g2g FGT %d hits --> accepted %d",nhits,ntot)<<endm;
00368 }
00369 
00370 
00371 //--------------------------------------------
00372 //--------------------------------------------
00373 //--------------------------------------------
00374 bool 
00375 StFgtSlowSimuMaker::projectQuadrant( int iquad) {
00376 
00377   double totAmp=digXY->Integral();
00378   LOG_INFO<<Form("::digiQuad()  totAmp=%g",totAmp)<<  endm;
00379   if(totAmp <par_stripAmplThres) return false;
00380   
00381   int nPix=0;
00382   
00383   double sumAmp=0; //for QA
00384   
00385   int nx=digXY->GetNbinsX();
00386   int ny=digXY->GetNbinsY();
00387   int bx,by;
00388   for(bx=1;bx<=nx;bx++)
00389     for(by=1;by<=ny;by++) {
00390       float amp=digXY->GetBinContent(bx,by); // this is not double because histo if float
00391       if(amp<par_2DpixAmplThres) continue; // drop too low signal
00392       nPix++;
00393       sumAmp+=amp;
00394       double x=digXY->GetXaxis()->GetBinCenter(bx); // SLOW
00395       double y=digXY->GetYaxis()->GetBinCenter(by); // SLOW
00396       int iRadID,iPhiID;// strip coordinates
00397       if( !geom->localXYtoStripID(iquad,x,y,iRadID, iPhiID)) continue;
00398 //        printf( "ampl=%f loc: x=%f y=%f iRad=%d iPhi=%d bx=%d by=%d\n",amp,x,y,iRadID, iPhiID,bx,by); 
00399      
00400       //accumulate response of individual strips
00401       digRad->Fill(iRadID,amp* mRadStripRelativeGain[iRadID]);
00402       digPhi->Fill(iPhiID,amp* mPhiStripRelativeGain[iPhiID]);      
00403 
00404       digRadAll->Fill(iRadID,amp* mRadStripRelativeGain[iRadID]);
00405       digPhiAll->Fill(iPhiID,amp* mPhiStripRelativeGain[iPhiID]);      
00406 
00407     }
00408   printf("digi: nPix=%d   sumAmp=%g\n", nPix,sumAmp); 
00409   
00410   return true;
00411 } 
00412   
00413 
00414 //--------------------------------------------
00415 //--------------------------------------------
00416 //--------------------------------------------
00417 void
00418 StFgtSlowSimuMaker::exportStripPlane(TH1F *h, vector<fgt_strip> &L) {
00419   float *y=h->GetArray();
00420   y++; // root histo counts bins from 1 - incredible silly
00421   float nx=h->GetNbinsX();
00422   for(int id=0;id<nx;id++,y++) {
00423     double ene=*y;
00424     if(ene <par_stripAmplThres) continue;
00425     fgt_strip st;
00426     st.id=id; // strips IDs are counted from zero
00427     st.adc=ene;
00428     L.push_back(st);
00429     //printf("%s add strID=%d adc=%.1f\n",h->GetName(),st.id,st.adc);
00430   }
00431 }
00432 
00433         
00434 
00435 
00436 
00439 
00440 // $Log: StFgtSlowSimuMaker.cxx,v $
00441 // Revision 1.5  2011/04/11 19:35:38  fisyak
00442 // Replace uint by UInt_t, use TMath
00443 //
00444 // Revision 1.4  2011/04/08 22:18:42  balewski
00445 // added access to TGeo
00446 //
00447 // Revision 1.3  2011/04/08 19:25:45  wzhang
00448 // Changed diskID assignment for Jan temporarily
00449 //
00450 // Revision 1.2  2011/04/08 01:14:13  balewski
00451 // removed most of FGT from ver 3
00452 //
00453 // Revision 1.1  2011/04/07 19:31:22  balewski
00454 // start
00455 //
00456 
00457 
00458  
00459 
00460