StRoot: StSpinPool/StWalgoB2009/St2009pubJSMaker.cxx Source File

00001 // $Id: St2009pubJSMaker.cxx,v 1.2 2010/03/14 22:50:31 balewski Exp $
00002 //
00003 //*-- Author :  Justin Stevens, IUCF
00004 // 
00005 
00006 #include "St2009WMaker.h"
00007 #include "St2009pubJSMaker.h"
00008 
00009 ClassImp(St2009pubJSMaker)
00010 
00011 //_______________________________________________________
00012 //_______________________________________________________
00013 St2009pubJSMaker::St2009pubJSMaker(const char *name):StMaker(name){
00014   wMK=0;HList=0;
00015   par_highEtow=5; // GeV
00016   par_awayNTrCut=2;
00017   
00018 }
00019 
00020 
00021 //_______________________________________________________
00022 //_______________________________________________________
00023 St2009pubJSMaker::~St2009pubJSMaker(){
00024   //
00025 }
00026 
00027 
00028 //_______________________________________________________
00029 //_______________________________________________________
00030 Int_t St2009pubJSMaker::Init(){
00031   assert(wMK);
00032   par_countTrPt=wMK->par_countTrPt;
00033   par_countTowEt=wMK->par_countTowEt;
00034   par_awayTotET=30; // old, clean it up, Jan
00035   assert(HList);
00036   initHistos();
00037   return StMaker::Init();
00038 }
00039 
00040 
00041 //_______________________________________________________
00042 //_______________________________________________________
00043 Int_t 
00044 St2009pubJSMaker::Make(){
00045   //printf("-----------in %s\n", GetName());
00046     
00047   doWanalysis();
00048   
00049   return kStOK;
00050 }
00051 
00052 //_______________________________________________________
00053 //_______________________________________________________
00054 void 
00055 St2009pubJSMaker::doWanalysis(){
00056   //has access to whole W-algo-maker data via pointer 'wMK'
00057  
00058   // run through W cuts to fill other histos............
00059   for(uint iv=0;iv<wMK->wEve.vertex.size();iv++) {
00060     WeveVertex &V=wMK->wEve.vertex[iv];
00061     for(uint it=0;it<V.eleTrack.size();it++) {
00062       if(it==0) etowQA(0,V.z); //if there is a reco track do etowQA
00063       WeveEleTrack &T=V.eleTrack[it];
00064       if(T.isMatch2Cl==false) continue;
00065       assert(T.cluster.nTower>0); // internal logical error
00066       assert(T.nearTotET>0); // internal logical error
00067       
00068       if(T.cluster.ET /T.nearTotET< wMK->par_nearTotEtFrac) continue; // too large nearET
00069       if(T.awayTotET> par_awayTotET) continue; // too large awayET, Jan
00070       //Transverse mass reco
00071 
00072       //hadronic recoil 
00073       TVector3 hadronicPt(T.hadronicRecoil.X(),T.hadronicRecoil.Y(),0); //transvers momentum vector
00074       hA[25]->Fill(hadronicPt.Perp()); 
00075       hA[26]->Fill(T.hadronicRecoil.Eta());
00076       hA[27]->Fill(T.cluster.ET);      
00077 
00078       TVector3 electronPt(T.cluster.position.X(),T.cluster.position.Y(),0); //transvers momentum vector
00079       electronPt.SetMag(T.cluster.ET);
00080 
00081       //neutrino reco 
00082       TVector3 neutrinoPt=-1*(hadronicPt+electronPt);
00083       hA[28]->Fill(neutrinoPt.Perp());  
00084 
00085       float recoDeltaPhi=neutrinoPt.DeltaPhi(electronPt);
00086       float Mt=sqrt(2*T.cluster.ET*neutrinoPt.Perp()*(1-cos(recoDeltaPhi)));  
00087       hA[29]->Fill(Mt); //real data
00088       hA[30]->Fill(Mt,T.cluster.ET);
00089 
00090       if(T.cluster.ET < 30) continue; //enrich W sample
00091       //Kinematics plots
00092       float eleEta=T.primP.Eta();
00093       if(eleEta < -0.8) 
00094         hA[31]->Fill(T.cluster.energy);
00095       if(eleEta > 0.8)
00096         hA[32]->Fill(T.cluster.energy);
00097       if(eleEta > -0.1 && eleEta < 0.1)
00098         hA[33]->Fill(T.cluster.energy);
00099 
00100       //Can we reco W pL?
00101       float eleTheta=T.primP.Theta();
00102       float ratioE=T.cluster.energy/40;  
00103       float pLRecoPlus=80*ratioE*((cos(eleTheta))+sqrt(cos(eleTheta)*cos(eleTheta)+sin(eleTheta)*sin(eleTheta)*(1-ratioE*ratioE)))/(ratioE*ratioE*sin(eleTheta)*sin(eleTheta));//+ solution
00104       float pLRecoMinus=80*ratioE*((cos(eleTheta))+sqrt(cos(eleTheta)*cos(eleTheta)+sin(eleTheta)*sin(eleTheta)*(1-ratioE*ratioE)))/(ratioE*ratioE*sin(eleTheta)*sin(eleTheta));//+ solution
00105       hA[34]->Fill(pLRecoPlus);
00106       hA[35]->Fill(pLRecoMinus);
00107 
00108 
00109       //Other stuff in pubJS
00110       hA[11]->Fill(T.nearNTow); 
00111       hA[12]->Fill(T.nearNTr);
00112 
00113       hA[24]->Fill(T.cluster.ET/T.nearTotET,T.cluster.position.PseudoRapidity());
00114       
00115       etowQA(1,-999); //before near cut do etowQA
00116       if(T.cluster.ET /T.nearTotET< wMK->par_nearTotEtFrac) continue; // too large nearET
00117       
00118       hA[7]->Fill(T.awayNTow); 
00119       hA[8]->Fill(T.awayNTr);
00120 
00121       etowQA(2,-999); //before away cut do etowQA
00122       if(T.awayTotET> par_awayTotET) { // too large awayET, Jan
00123         if(T.awayNTr > par_awayNTrCut)  
00124           hA[9]->Fill(T.cluster.ET);   
00125         else                           
00126           hA[10]->Fill(T.cluster.ET);  
00127         continue; 
00128       }
00129       
00130       if(T.cluster.ET < 25){ //should be all background
00131         hA[13]->Fill(V.z);
00132         hA[14]->Fill(T.awayTotET);
00133         hA[15]->Fill(T.nearNTow); 
00134         hA[16]->Fill(T.nearNTr);
00135         
00136         hA[18]->Fill(V.id);
00137         hA[19]->Fill(V.funnyRank);
00138         hA[20]->Fill(wMK->wEve.bx48);
00139         hA[21]->Fill(wMK->wEve.bx7);
00140         hA[22]->Fill(T.cluster.position.PseudoRapidity(),T.cluster.position.Phi());
00141         hA[23]->Fill(T.cluster.nTower);
00142       }
00143 
00144     }
00145   }
00146 }
00147 
00148 
00149 //________________________________________________
00150 //________________________________________________
00151 void
00152 St2009pubJSMaker::etowQA(int whichCut,float zVert){//whichCut gives where in algo you are calling etowQA
00153   
00154   //loop over all endcap towers
00155   for(int iphi=0;iphi<60;iphi++){
00156     for(int ieta=0;ieta<12;ieta++){
00157       TVector3 towerPos=wMK->positionEtow[iphi][ieta];
00158       float ene=wMK->wEve.etow.ene[iphi][ieta];
00159       int etabin=abs(ieta-11);//reverse ie. now etabin[0,11] -> phys eta[1,2]
00160       //if(iphi==0) cout<<"iphi = "<<iphi<<" phi = "<<towerPos.Phi()<<" ieta = "<<ieta<<" eta = "<<towerPos.Eta()<<" etabin = "<<etabin<<endl;
00161 
00162       //get ET for test
00163       TVector3 primP=towerPos-TVector3(0,0,zVert);
00164       primP.SetMag(ene); // it is 3D momentum in the event ref frame
00165       float ET=primP.Perp();
00166 
00167       hB[whichCut]->Fill(etabin,towerPos.Phi(),ene); 
00168       if(ET > par_highEtow) //high tower
00169         hA[whichCut]->Fill(etabin,towerPos.Phi());
00170     
00171       // Look at vertex Z dependence of energy for all reco vertices that have a track with pT > 10
00172       if(whichCut==0 && etabin==11){//only once for all events
00173         if(zVert>-100  && zVert<=-50)
00174           hA[3]->Fill(ene);
00175         else if(zVert>-50  && zVert<=0)
00176           hA[4]->Fill(ene);
00177         else if(zVert>0  && zVert<=50)
00178           hA[5]->Fill(ene);
00179         else if(zVert>50  && zVert<=100)
00180           hA[6]->Fill(ene);
00181       }
00182 
00183 
00184     }
00185   }
00186   return;
00187 }
00188 
00189 
00190 
00191 // $Log: St2009pubJSMaker.cxx,v $
00192 // Revision 1.2  2010/03/14 22:50:31  balewski
00193 // *** empty log message ***
00194 //
00195 // Revision 1.1  2009/11/23 23:00:18  balewski
00196 // code moved spin-pool
00197 //
00198 // Revision 1.1  2009/11/23 21:11:18  balewski
00199 // start
00200 //