StRoot: StEEmcPool/StEEmcIUPi0/StMcOutputMaker.cxx Source File

00001 //*-- Author :Weihong He
00002 // 
00003 // $Id: StMcOutputMaker.cxx,v 1.3 2011/04/11 19:35:41 fisyak Exp $
00004 #include <TH2.h>
00005 
00006 #include "StMcOutputMaker.h"
00007 #include "StChain.h"
00008 #include "St_DataSetIter.h"
00009 #include "StMcEventMaker/StMcEventMaker.h"
00010 #include "StMcEvent/StMcEvent.hh"
00011 #include "StMcEvent/StMcVertex.hh"
00012 #include "StMcEvent/StMcTrack.hh"
00013 #include "StEEmcUtil/EEmcGeom/EEmcGeomSimple.h"
00014 //MuDst
00015 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
00016 #include "StMuDSTMaker/COMMON/StMuDst.h"
00017 #include "StMuDSTMaker/COMMON/StMuTrack.h"
00018 
00019 #include "StarClassLibrary/StParticleDefinition.hh"
00020 
00021 
00022 
00023 ClassImp(StMcOutputMaker)
00024 
00025 
00026 StMcOutputMaker::StMcOutputMaker(const char *name):StMaker(name){
00027   memset(hA,0,sizeof(hA));
00028   HList=0;
00029   geomE= new EEmcGeomSimple();
00030 
00031 }
00032 
00033 
00034 StMcOutputMaker::~StMcOutputMaker(){
00035   //
00036 }
00037 
00038 
00039 //_____________________________________________________________________________
00040 //
00041 Int_t 
00042 StMcOutputMaker::Init(){
00043 
00044 
00045   memset(hA,0,sizeof(hA));
00046   hA[0]=new TH2F("gEEmcEtaPhi", "Generated pi0 track; phi / deg;EEmc Eta",360,-180.,180,20,1.,2.);
00047   //hA[1]=new TH2F("rEtaPhi", "Recon pi0; track eta; phi / deg",20,0.5,2.5,20,0,90);
00048   hA[1]=new TH1F("gPt", "Generated pi0 Pt ",50,0.,25.);
00049   hA[2]=new TH1F("gEta", "Generated pi0 eta ",50,0.,2.5);
00050   hA[3]=new TH1F("gEnergy","generated Mc pi0 energy",80,0.,80.);
00051   hA[4]=new TH1F("gZgg","generated Mc pi0 energy sharing",100,0.,1.);
00052   hA[5]=new TH1F("gZgg01","generated Mc pi0 energy sharing[5,10]",100,0.,1.);
00053   hA[6]=new TH1F("gZgg02","generated Mc pi0 energy sharing[10,15]",100,0.,1.);
00054   hA[7]=new TH1F("gZgg03","generated Mc pi0 energy sharing[15,20]",100,0.,1.);
00055   hA[8]=new TH1F("gZgg04","generated Mc pi0 energy sharing[20,30]",100,0.,1.);
00056   hA[9]=new TH1F("gZgg05","generated Mc pi0 energy sharing[30,60]",100,0.,1.);
00057   hA[10]=new TH1F("gPhi", "Generated pi0 Phi ",360,-180.,180.);
00058   hA[11]=new TH1F("gEEmcEta", "Generated pi0 eta in EEMC ",50,0.,2.5);
00059   hA[12]=new TH1F("DeltaEta", "genEta - EEmcEta ",40,-2.,2.);
00060   hA[13]=new TH1F("gZvert", "generated Z vertex ",300,-150.,150.);
00061   //hA[3]=new TH1F("rEta", "Recon pi0 eta ",40,0.5,2.5);
00062   //hA[4]=new TH1F("rNP", "nFitPoints,  Recon pi0; nFitP ",50,-0.5,49.5);
00063   //hA[5]=new TH1F("rrpt", "Pt resolution, ptGen=6 GeV ; gen-reco/gen",100,-0.5,0.5);
00064  
00065   assert(HList);
00066   int i;
00067   for(i=0;i<mxHa;i++) if(hA[i]) HList->Add(hA[i]);
00068 
00069    return StMaker::Init();
00070 }
00071 
00072 //________________________________________________
00073 //________________________________________________
00074 void
00075 StMcOutputMaker::Clear(const Option_t*){
00076   gTr.clear();
00077   geemcEta.clear();
00078   genZZ.clear();
00079   genXX.clear();
00080   genYY.clear();
00081 }
00082 
00083 //________________________________________________
00084 //________________________________________________
00085 Int_t 
00086 StMcOutputMaker::Make(){
00087   //
00088   //PrintInfo();
00089   //
00090   int fflag=0;
00091   float hlow=270.0*tan(15.0*3.14159265/180.0);
00092   float hhigh=270.0*tan(40.0*3.14159265/180.0);
00093   //std::cout<<"hlow="<<hlow<<"hhigh="<<hhigh<<std::endl;
00094 
00095 
00096   probTr= StLorentzVectorF();
00097   
00098   // **************************************
00099   // get pointer to mcEventMaker, Event *
00100   // **************************************
00101   StMcEvent* mMcEvent = 0;  
00102   mMcEvent = (StMcEvent*) StMaker::GetChain()->GetDataSet("StMcEvent"); 
00103   assert(mMcEvent); 
00104 
00105   //#if 0
00106   //StMcVertex*  PV=  mMcEvent->       primaryVertex();
00107   //printf("%s:: nPrimV=%p\n", GetName(), (void*)V);
00108   //assert(PV);
00109   //int nPD=PV-> numberOfDaughters();
00110   //#endif 
00111   //printf("n prim Doughters=%d\n",nPD);
00112   //all vertices
00113   const StSPtrVecMcVertex &VL=mMcEvent->vertices();
00114   //printf("%s::nVert=%d\n",GetName(),VL.size());
00115   //  int ipr=1;
00116   UInt_t i;
00117   zgg=0;
00118   for( i=0;i<VL.size();i++) { 
00119     StMcVertex* V=VL[i];
00120     
00121     // Rxy vx z vertex to remove fake electrons
00122     float x=V->position().x();
00123     float y=V->position().y();
00124     float zz=V->position().z();
00125     float Rxy=TMath::Sqrt(x*x+y*y);
00126     
00127     //printf("vx=%f,vy=%f,vz=%f Rxy=%f\n",x,y,z,Rxy);
00128     //cout<<*V<<endl<<endl;
00129     int nD=V-> numberOfDaughters();
00130     //cout <<Form("#daughters=%d\n",nD)<<endl;
00131     for(int it=0; it<nD;it++ ) {
00132       StMcTrack* tr=V->daughter(it);
00133       int gid=tr->geantId();
00134       if(gid!=7)
00135         {continue;}
00136       float hHeight=0.0;
00137       StLorentzVectorF p4=tr->fourMomentum();
00138       float psRap=tr->pseudoRapidity();
00139       if( psRap<0 || psRap>3.) continue; // reject sputnik-muons
00140       //  cout<<*tr<<endl;
00141       probTr=p4;
00142       float ppt=p4.perp();
00143       if(ppt<=3.0) continue;
00144       float ppz=p4.pz();
00145       hHeight=ppt*(270.0-zz)/ppz+Rxy;
00146       
00147       if(hHeight<=hlow || hHeight >=hhigh) continue;
00148       fflag++;
00149       float etatheta=TMath::ATan(hHeight/270.0);
00150       //printf("theta=%f etatheta=%f\n",p4.theta(), etatheta);
00151       float mideta=tan(etatheta/2.0);
00152       float eemceta=-log(mideta);
00153       //printf("geta=%f eemceta=%f zv=%f\n",tr->pseudoRapidity(),eemceta,zz);
00154       //printf(" GEANT prob: gid=%d  E/GeV=%.1f theta/deg=%.1f phi/deg=%.1f eta=%.2f\n",gid,p4.e(),p4.theta()/3.1416*180,p4.phi()/3.1416*180 ,tr->pseudoRapidity());
00155       hA[0]->Fill(p4.phi()/3.1416*180,eemceta);
00156       hA[2]->Fill(tr->pseudoRapidity());
00157       hA[10]->Fill(p4.phi()/3.1416*180);
00158       hA[1]->Fill(p4.perp());
00159       hA[3]->Fill(p4.e());
00160       hA[11]->Fill(eemceta);
00161       hA[12]->Fill(tr->pseudoRapidity()-eemceta);
00162       hA[13]->Fill(zz);
00163       //printf("eta=%f eta2=%f\n",p4.pseudoRapidity(),tr->pseudoRapidity());
00164       genE=0;
00165       genEta=0;
00166       genPhi=0;
00167       if(p4.e()>=0)
00168         {
00169           genEta=tr->pseudoRapidity();
00170           genE=p4.e();
00171           genPhi=p4.phi()/3.1416*180;
00172         }
00173       gTr.push_back(tr);
00174       geemcEta.push_back(eemceta);
00175       genZZ.push_back(zz);
00176       genXX.push_back(x);
00177       genYY.push_back(y);
00178       //printf("zvertex=%f pt=%f pz=%f height=%f size=%d\n",zz,ppt,ppz,hHeight,gTr.size());
00179     }
00180     if(nD==2)
00181     {
00182       StMcTrack* tr1=V->daughter(0);
00183       int geid1=tr1->geantId();
00184       StMcTrack* tr2=V->daughter(1);
00185       int geid2=tr2->geantId();
00186       float Etot=0,Edif=0;
00187       if((geid1==1)&&(geid2==1))
00188         {
00189           for(int it=0; it<nD;it++ ) { // primary tracks
00190             StMcTrack* tr=V->daughter(it);
00191             //int geid=tr->geantId();
00192             //float phi=tr->momentum().phi(); 
00193             //float eta=tr->pseudoRapidity();              
00194             //float pt=tr->pt(); 
00195             float ener=tr->fourMomentum().e();
00196             //printf("geid=%d eta=%f phi=%f pt=%f e=%f\n",geid,eta,phi,pt,ener);
00197             Etot+=ener;
00198             Edif=ener-Edif;
00199           }
00200           zgg=fabs(Edif)/Etot;
00201           //printf("zgg=%f\n",zgg);
00202           hA[4]->Fill(zgg);
00203           if(Etot>=5.0 && Etot<=10.0)
00204             {
00205               hA[5]->Fill(zgg);
00206             }
00207           if(Etot>=10.0 && Etot<=15.0)
00208             {
00209               hA[6]->Fill(zgg);
00210             }
00211           if(Etot>=15.0 && Etot<=20.0)
00212             {
00213               hA[7]->Fill(zgg);
00214             }
00215           if(Etot>=20.0 && Etot<=30.0)
00216             {
00217               hA[8]->Fill(zgg);
00218             }
00219           if(Etot>=30.0 && Etot<=60.0)
00220             {
00221               hA[9]->Fill(zgg);
00222             }
00223         }
00224     }
00225   }
00226       //  printf("x=%f,y=%f,z=%f,Rxy=%f\n",x,y,z,Rxy);  
00227 
00228   //printf("found %d pi0\n",gTr.size());
00229   return kStOK;
00230 }
00231 
00232 
00233 
00234 
00235 
00236 
00237 
00238