StRoot: StEmcPool/StPhotonAnalysis/AnaCuts.cxx Source File

00001 #include <Riostream.h>
00002 #include <TString.h>
00003 #include <TMath.h>
00004 #include <assert.h>
00005 
00006 #include <StEmcPool/StPhotonCommon/MyEvent.h> 
00007 #include <StEmcPool/StPhotonCommon/MyPoint.h> 
00008 //#include <StEmcPool/StPhotonCommon/MyMcTrack.h> 
00009 
00010 #include "AnaCuts.h"
00011 
00012 ClassImp(AnaCuts)
00013 
00014 AnaCuts::AnaCuts(const char* flag)
00015 {
00016   cout<<"constructing CUTS for "<<flag<<endl;
00017 
00018   if(strcmp(flag,"dAu")==0){
00019 
00020     vertexxCUT=5.0;
00021     vertexyCUT=5.0;
00022     vertexzCUT=60.;
00023     neutralCUT=5.;//10, note this is cm's, 1 strip = 1.7cm or so.
00024     photonCUT=15.;
00025     etaHitsCUT=2;//2;
00026     phiHitsCUT=1;//2;
00027     energyRatioCUT=0.;//0
00028   
00029     etaMinCUT=.1;//.1
00030     etaMaxCUT=.9;//.9
00031     rapidityMinCUT=.1;//0.
00032     rapidityMaxCUT=.9;//1.
00033     rapPionMinCUT=.1;//.1
00034     rapPionMaxCUT=.9;//.9
00035     asymmetryCUT=0.7;//.7
00036     ratioCUT=0.8;//.8
00037 
00038     softTrigHT1=0.;
00039     softTrigHT2=0.;
00040 
00041     nPtBinsMB=8;
00042     ptBinsMB.Set(nPtBinsMB+1);
00043     for(Int_t i=0;i<=nPtBinsMB;i++) ptBinsMB[i]=(Double_t)0.5*i;
00044     ptBinsMB[7]=4.;
00045     ptBinsMB[8]=5.;
00046 
00047     nPtBinsHT1=14;
00048     ptBinsHT1.Set(nPtBinsHT1+1);
00049     for(Int_t i=0;i<=nPtBinsHT1;i++) ptBinsHT1[i]=(Double_t)1.*i;
00050     ptBinsHT1[13]=(Double_t)13.5;
00051     ptBinsHT1[14]=(Double_t)15.;
00052 
00053     nPtBinsHT2=14;
00054     ptBinsHT2.Set(nPtBinsHT2+1);
00055     for(Int_t i=0;i<=nPtBinsHT2;i++) ptBinsHT2[i]=(Double_t)1.*i;
00056     ptBinsHT2[13]=(Double_t)13.5;
00057     ptBinsHT2[14]=(Double_t)15.;
00058 
00059 
00060     nPtBinsEffMB=nPtBinsMB;//MB
00061     ptBinsEffMB=ptBinsMB;//MB
00062     nPtBinsEffHT1=nPtBinsHT1;
00063     ptBinsEffHT1=ptBinsHT1;
00064     nPtBinsEffHT2=nPtBinsHT2;
00065     ptBinsEffHT2=ptBinsHT2;
00066      
00067     nMinvBinsMB=500;
00068     mInvLowMB=0.0;
00069     mInvHighMB=5.0;
00070     nMinvBinsHT1=250;
00071     mInvLowHT1=0.0;
00072     mInvHighHT1=5.0;
00073     nMinvBinsHT2=250;
00074     mInvLowHT2=0.0;
00075     mInvHighHT2=5.0;
00076 
00077     nPtBinsEff=40;
00078     ptBinsEff.Set(nPtBinsEff+1);
00079     for(Int_t i=0;i<=nPtBinsEff;i++){
00080       ptBinsEff[i]=i*0.5;
00081     }
00082 
00083     ht1AdcCUT=287;
00084     ht2AdcCUT=447;
00085 
00086     isHot.Set(4800);
00087     isHot.Reset();
00088     Int_t badIdsDAU[]={82, 156, 240, 264, 719, 733, 818, 853, 931, 1289, 1355, 1358, 1359, 1472,
00089                        1527, 1759, 1943, 1965, 2005, 2011, 2025, 2090, 2130, 2131, 2152, 2170, 2369, -1};
00090     Int_t entryDAU=0;
00091     while(badIdsDAU[entryDAU]>-1) {isHot[badIdsDAU[entryDAU]-1]=1; entryDAU++;}
00092     
00093     timesSigma=2.0;
00094   }
00095   else if(strcmp(flag,"pp05")==0 || strcmp(flag,"pythia")==0){
00096 
00097     vertexxCUT=5.0;
00098     vertexyCUT=5.0;
00099     vertexzCUT=60.;//60.
00100     neutralCUT=5.;//10, note this is cm's, 1 strip = 1.7cm or so.
00101     photonCUT=15.;
00102     etaHitsCUT=2;//2;
00103     phiHitsCUT=1;//1;
00104     energyRatioCUT=0.;//0
00105   
00106     etaMinCUT=.2;//.1
00107     etaMaxCUT=.8;//.9
00108     rapidityMinCUT=.1;//0.
00109     rapidityMaxCUT=.9;//1.
00110     rapPionMinCUT=.1;//.1
00111     rapPionMaxCUT=.9;//.9
00112     asymmetryCUT=.7;//.7
00113     ratioCUT=1.;
00114 
00115     softTrigHT1=0.;
00116     softTrigHT2=0.;
00117 
00118     nPtBinsMB=8;
00119     ptBinsMB.Set(nPtBinsMB+1);
00120     for(Int_t i=0;i<=nPtBinsMB;i++) ptBinsMB[i]=(Double_t)0.5*i;
00121     ptBinsMB[7]=4.0;
00122     ptBinsMB[8]=5.0;
00123 
00124 
00125     nPtBinsHT1=14;
00126     ptBinsHT1.Set(nPtBinsHT1+1);
00127     for(Int_t i=0;i<=nPtBinsHT1;i++) ptBinsHT1[i]=(Double_t)1.*i;
00128     ptBinsHT1[13]=(Double_t)13.5;
00129     ptBinsHT1[14]=(Double_t)15.0;
00130 
00131     nPtBinsHT2=14;
00132     ptBinsHT2.Set(nPtBinsHT2+1);
00133     for(Int_t i=0;i<=nPtBinsHT2;i++) ptBinsHT2[i]=(Double_t)1.*i;
00134     ptBinsHT2[13]=(Double_t)13.5;
00135     ptBinsHT2[14]=(Double_t)15.0;
00136 
00137     nPtBinsEffMB=nPtBinsMB;//MB
00138     ptBinsEffMB=ptBinsMB;//MB
00139     nPtBinsEffHT1=nPtBinsHT1;
00140     ptBinsEffHT1=ptBinsHT1;
00141     nPtBinsEffHT2=nPtBinsHT2;
00142     ptBinsEffHT2=ptBinsHT2;
00143         
00144     nMinvBinsMB=500;
00145     mInvLowMB=0.0;
00146     mInvHighMB=5.0;
00147     nMinvBinsHT1=250;
00148     mInvLowHT1=0.0;
00149     mInvHighHT1=5.0;
00150     nMinvBinsHT2=250;
00151     mInvLowHT2=0.0;
00152     mInvHighHT2=5.0;
00153 
00154     nPtBinsEff=40;
00155     ptBinsEff.Set(nPtBinsEff+1);
00156     for(Int_t i=0;i<=nPtBinsEff;i++){
00157       ptBinsEff[i]=i*0.5;
00158     }
00159 
00160     ht1AdcCUT=0;//420;//prelim. rejection of obvious bad triggers
00161     ht2AdcCUT=0;//550;//but for pp05 trigger must be simulated.
00162     
00163     isHot.Set(4800);
00164     isHot.Reset();
00165     Int_t badIdsPP[]={750, 757, 762, 768, 812, -1};
00166     Int_t entryPP=0;
00167     while(badIdsPP[entryPP]>-1) {isHot[badIdsPP[entryPP]-1]=1; entryPP++;}
00168     
00169     timesSigma=2.0;
00170   }
00171   else{
00172     cout<<"error constructing cuts"<<endl;
00173     assert(0);
00174   }
00175 
00176 
00177   if(0){
00178     nPtBinsEffMB=40;
00179     ptBinsEffMB.Set(nPtBinsEffMB+1);
00180     for(Int_t i=0;i<=nPtBinsEffMB;i++) ptBinsEffMB[i]=(Double_t).5*i;
00181     nPtBinsEffHT1=40;
00182     ptBinsEffHT1.Set(nPtBinsEffHT1+1);
00183     for(Int_t i=0;i<=nPtBinsEffHT1;i++) ptBinsEffHT1[i]=(Double_t).5*i;
00184     nPtBinsEffHT2=40;
00185     ptBinsEffHT2.Set(nPtBinsEffHT2+1);
00186     for(Int_t i=0;i<=nPtBinsEffHT2;i++) ptBinsEffHT2[i]=(Double_t).5*i;
00187   }
00188 
00189 }
00190 AnaCuts::~AnaCuts()
00191 {
00192   cout<<"cuts destructed!"<<endl;
00193 }
00194 Bool_t AnaCuts::isPointOK(MyPoint *p,TVector3 vpos)
00195 {
00196   TVector3 pos=p->position();
00197   TVector3 mom=pos-vpos;
00198   mom.SetMag(p->energy());
00199   if(p->distanceToTrack()<neutralCUT) return kFALSE;
00200   if(!(p->nHitsEta()>=etaHitsCUT && p->nHitsPhi()>=phiHitsCUT))
00201     {
00202       if(!(p->nHitsEta()>=phiHitsCUT && p->nHitsPhi()>=etaHitsCUT))
00203         {
00204           return kFALSE;
00205         }
00206     }
00207   if(pos.PseudoRapidity()<etaMinCUT) return kFALSE;
00208   if(pos.PseudoRapidity()>etaMaxCUT) return kFALSE;
00209 
00210   //new
00211   if((p->energyEta()+p->energyPhi())/p->energy()<energyRatioCUT) return kFALSE;
00212  
00213   return kTRUE;
00214 }
00215 
00216 Bool_t AnaCuts::isEventOK(MyEvent *ev, const char *flag)
00217 {
00218   TVector3 vPos=ev->vertex();
00219   Float_t ratioE=TMath::Abs(ev->energyInBarrel()+ev->momentumInTpc())>0. ?
00220     ev->energyInBarrel()/(ev->energyInBarrel()+ev->momentumInTpc()) : 0.;
00221   //event cuts
00222   if(strcmp(flag,"dAu")==0){
00223     if(TMath::Abs(vPos.Z())>vertexzCUT) return kFALSE;
00224   }
00225   if(strcmp(flag,"pp05")==0){
00226     if(TMath::Abs(ev->bbcVertexZ())>vertexzCUT) return kFALSE;
00227   }
00228   if(TMath::Abs(vPos.X())>vertexxCUT) return kFALSE;
00229   if(TMath::Abs(vPos.Y())>vertexyCUT) return kFALSE;
00230 
00231   if(ev->trigger()<1) return kFALSE;
00232   //hot towers and beam bg:
00233   if(ev->trigger()&2 || ev->trigger()&4){
00234     if(ev->highTowerId()>0 && isHot[ev->highTowerId() - 1]){
00235       if(ev->trigger()&1) ev->setTrigger(1);
00236       else return false;
00237     }
00238     if(ev->highTowerEnergy()<0.001){
00239       if(ev->trigger()&1) ev->setTrigger(1);
00240       else return false;
00241     }
00242     if(ratioE>ratioCUT){
00243       if(ev->trigger()&1) ev->setTrigger(1);
00244       else return false;
00245     }
00246   }
00247   if(ev->trigger()&4 && ev->highTowerAdc()<=ht2AdcCUT){
00248     ev->setTrigger(ev->trigger()-4);
00249   }
00250   if(ev->trigger()&2 && ev->highTowerAdc()<=ht1AdcCUT){
00251     ev->setTrigger(1);
00252   }
00253 
00254   return true;
00255 }
00256 void AnaCuts::printCuts()
00257 {
00258   cout<<"----- event cuts ------"<<endl;
00259   cout<<"vertexxCUT="<<vertexxCUT<<endl;
00260   cout<<"vertexyCUT="<<vertexyCUT<<endl;
00261   cout<<"vertexzCUT="<<vertexzCUT<<endl;
00262   cout<<"ratioCUT="<<ratioCUT<<endl<<endl;
00263 
00264   cout<<"----- point cuts ------"<<endl;
00265   cout<<"neutralCUT="<<neutralCUT<<endl;
00266   cout<<"chargedCUT="<<photonCUT<<endl;
00267   cout<<"etaHitsCUT="<<etaHitsCUT<<endl;
00268   cout<<"phiHitsCUT="<<phiHitsCUT<<endl;
00269   cout<<"etaMinCUT="<<etaMinCUT<<endl;
00270   cout<<"etaMaxCUT="<<etaMaxCUT<<endl;
00271   cout<<"rapidityMinCUT="<<rapidityMinCUT<<endl;
00272   cout<<"rapidityMaxCUT="<<rapidityMaxCUT<<endl;
00273   cout<<"rapPionMinCUT="<<rapPionMinCUT<<endl;
00274   cout<<"rapPionMaxCUT="<<rapPionMaxCUT<<endl;
00275   cout<<"asymmetryCUT="<<asymmetryCUT<<endl;
00276   cout<<"ht1AdcCUT="<<ht1AdcCUT<<endl;
00277   cout<<"ht2AdcCUT="<<ht2AdcCUT<<endl<<endl;
00278 
00279 
00280   cout<<endl<<"------ binning ------"<<endl;
00281   
00282   //bins
00283   cout<<"mb bins={";
00284   for(Int_t i=0;i<=nPtBinsMB;i++) cout<<ptBinsMB[i]<<",";
00285   cout<<"};"<<endl;
00286   cout<<"ht1 bins={";
00287   for(Int_t i=0;i<=nPtBinsHT1;i++) cout<<ptBinsHT1[i]<<",";
00288   cout<<"};"<<endl;
00289   cout<<"ht2 bins={";
00290   for(Int_t i=0;i<=nPtBinsHT2;i++) cout<<ptBinsHT2[i]<<",";
00291   cout<<"};"<<endl;
00292   cout<<"Eff bins={";
00293   for(Int_t i=0;i<=nPtBinsEff;i++) cout<<ptBinsEff[i]<<",";
00294   cout<<"};"<<endl;
00295 
00296   cout<<"inv mass bins mb: "<<nMinvBinsMB<<" "<<mInvLowMB<<" "<<mInvHighMB<<endl;
00297   cout<<"inv mass bins ht1: "<<nMinvBinsHT1<<" "<<mInvLowHT1<<" "<<mInvHighHT1<<endl;
00298   cout<<"inv mass bins ht2: "<<nMinvBinsHT2<<" "<<mInvLowHT2<<" "<<mInvHighHT2<<endl;  
00299 
00300   //hot towers:
00301   cout<<endl<<endl<<"------- hot towers -------"<<endl;
00302   for(Int_t k=1;k<=4800;k++)
00303     {
00304       if(isHot[k-1]) cout<<k<<", ";
00305     }
00306   cout<<endl<<"-----------------------"<<endl<<endl<<endl;
00307 }