manTreeSFrame: /afs/hep.man.ac.uk/u/markowen/ATLAS/SFrameProof/2011EPS/topUtils/src/SmearingClass

00001 
00002 // SmearingClass.h -- ATLAS Experiment Software //
00014 #include <cstring>
00015 #include <iostream>
00016 #include <stdio.h>
00017 #include <TROOT.h>
00018 #include <math.h>
00019 #include "TRandom3.h"
00020 
00021 
00022 class SmearingClass{
00023  public:
00024   /*Constructors*/
00025   SmearingClass() {Initialize();}
00026   
00027   
00028   /************/
00029   /* Methods */    
00030   /**********/
00031   
00033   void Initialize() {    
00034     detType="";
00035     GeV=1000;
00036     gRand = new TRandom3();
00037     gRand->SetSeed(0);
00038     useTan2=true;
00039     useScale=false;
00040     ptms=0;
00041     ptid=0;
00042     ptcb=0;
00043     eta=0;
00044   }
00045   
00048   void SetSeed(int seed) {gRand->SetSeed(seed);}
00049   void SetSeed(int evtnum, int muon_index, int offset=680049) {gRand->SetSeed(offset + evtnum + muon_index*100);}
00050   void UseTan2(bool U) {useTan2=U;}/* use pT/tan(theta)^2 parameterization for CSC region */  
00051   void UseGeV() {GeV=1;}  /* if input momenta are in GeV */
00052   void UseScale(bool applyScale) {useScale=applyScale;}  /* also use momentum scale, set as default */
00053 
00054   
00055   /* MS or ID smearing only. DetType="MS","ID"*/
00056   void Event(double Pt, double Eta,std::string DetType) 
00057   {        
00058     ptms=0;
00059     ptid=0;
00060     ptcb=0;
00061     
00062     if (DetType=="MS")      ptms = Pt;
00063     else if (DetType=="ID") ptid = Pt;
00064     else std::cout<<"SmearingClass ERROR: wrong DetType in input "<<DetType<<std::endl; 
00065     eta=Eta;    
00066     detType=DetType;
00067     Event();
00068   }
00069   
00070   /* For full MS ID and CB smearing */
00071   void Event(double PtMS, double PtID, double PtCB, double Eta) 
00072   {
00073     detType="All";
00074     ptms=PtMS;
00075     ptid=PtID;
00076     ptcb=PtCB;
00077     eta=Eta;
00078     Event();
00079   }  
00080   
00081   void  Event()  {
00082     smearMS=0;
00083     smearID=0;
00084     smearCB=0;
00085     
00086     /* Detector Region */
00087     if (fabs(eta)<1.05)     detRegion=0;
00088     else if (fabs(eta)<1.7) detRegion=1;
00089     else if (fabs(eta)<2.0) detRegion=2;
00090     else if (fabs(eta)<2.5) detRegion=3;
00091     else detRegion=-1;
00092     
00093     g1 = gRand->Gaus(0,1);  
00094     g2 = gRand->Gaus(0,1); 
00095     g3 = gRand->Gaus(0,1); 
00096     g4 = gRand->Gaus(0,1); 
00097     
00099     smearMS=Smearing("MS");
00100     smearID=Smearing("ID");
00101     if (detType=="All") smearCB=Combine(smearMS,smearID);
00102     ErrorMatrix();    
00103   }
00104   
00105 
00106   double Smearing(std::string DetType) {
00107     
00108     if (detRegion<0||detRegion>3) return 0;
00109     
00110     if (DetType=="MS") {return (MS_MS[detRegion]*g1+MS_AL[detRegion]*g2*ptms/GeV);}
00111     else if(DetType=="ID") 
00112       {
00113         if (useTan2&&detRegion==3) {return (ID_MS[detRegion]*g4+ID_ALTAN[detRegion]*g3*ptid/GeV*sinh(eta)*sinh(eta));}
00114         else {return (ID_MS[detRegion]*g4+ID_AL[detRegion]*g3*ptid/GeV);}
00115       } 
00116     else std::cout<<"Error:: DetType not defined "<<DetType<<std::endl;
00117     return 0;
00118   }
00119   
00120   double Combine(double SmearMS,double SmearID) {
00123     if (detRegion<0||detRegion>3) return 0;
00124     if (ptcb==0) {std::cout<<"Warning: ptcb==0"<<std::endl;return 0;}
00125     
00126     double SigmaMS=pow(
00127                        pow(MC_MS_CALO[detRegion]/ptcb*GeV,2)+
00128                        pow(MC_MS_MS[detRegion],2)+
00129                        pow(MC_MS_AL[detRegion]*ptcb/GeV,2),0.5);
00130     
00131     double SigmaID=pow(
00132                        pow(MC_ID_MS[detRegion],2)+
00133                        pow(MC_ID_AL[detRegion]*ptcb/GeV,2),0.5);
00134     
00135     if (detRegion==3&&useTan2) {
00136       SigmaID=pow(
00137                   pow(MC_ID_MS[detRegion],2)+
00138                   pow(MC_ID_ALTAN[detRegion]*ptcb/GeV*sinh(eta)*sinh(eta),2),0.5);
00139     }
00140     
00141     
00142     wMS=1./SigmaMS/SigmaMS;
00143     wID=1./SigmaID/SigmaID;
00144     return (SmearMS*wMS+SmearID*wID)/(wMS+wID);
00145   }
00146   
00147   void ErrorMatrix() {
00148     
00149     vms=0;
00150     vid=0;
00151     corr=0;
00152     
00153     if (detRegion<0||detRegion>3) return;
00154     
00155     
00156     if (!useTan2||detRegion!=3){
00157       double s1=pow(E_MS_MS[detRegion]*E_MS_MS[detRegion]+S_MS_MS[detRegion]*S_MS_MS[detRegion],0.5);
00158       double s2=pow(E_MS_AL[detRegion]*E_MS_AL[detRegion]+S_MS_AL[detRegion]*S_MS_AL[detRegion],0.5);
00159       double s3=pow(E_ID_AL[detRegion]*E_ID_AL[detRegion]+S_ID_AL[detRegion]*S_ID_AL[detRegion],0.5);
00160       vms=g1*g1*s1*s1+
00161         g2*g2*ptms/GeV*ptms/GeV*s2*s2+
00162         2.*g1*g2*ptms/GeV*s1*s2*CorrMat[detRegion][0];
00163       vms=pow(vms,0.5);
00164       vid=g3*g3*ptid/GeV*ptid/GeV*s3*s3;
00165       vid=pow(vid,0.5);
00166       if (vms*vid!=0)
00167         corr=(g1*g3*s1*s3*ptid/GeV*CorrMat[detRegion][1]+g2*g3*ptms/GeV*ptid/GeV*s2*s3*CorrMat[detRegion][2])/vms/vid;
00168     } else {
00169       double s1=pow(E_MS_MS[detRegion]*E_MS_MS[detRegion]+S_MS_MS[detRegion]*S_MS_MS[detRegion],0.5);
00170       double s2=pow(E_MS_AL[detRegion]*E_MS_AL[detRegion]+S_MS_AL[detRegion]*S_MS_AL[detRegion],0.5);
00171       double s3=pow(E_ID_ALTAN[detRegion]*E_ID_ALTAN[detRegion]+S_ID_ALTAN[detRegion]*S_ID_ALTAN[detRegion],0.5);
00172       vms=g1*g1*s1*s1+
00173         g2*g2*ptms/GeV*ptms/GeV*s2*s2+
00174         2.*g1*g2*ptms/GeV*s1*s2*CorrMatTan[detRegion][0];
00175       vms=pow(vms,0.5);
00176       vid=g3*g3*ptid/GeV*ptid/GeV*s3*s3*sinh(eta)*sinh(eta)*sinh(eta)*sinh(eta);
00177       vid=pow(vid,0.5);
00178       if (vms*vid!=0)
00179         corr=(g1*g3*s1*s3*ptid/GeV*sinh(eta)*sinh(eta)*CorrMatTan[detRegion][1]+
00180               g2*g3*ptms/GeV*sinh(eta)*sinh(eta)*ptid/GeV*s2*s3*CorrMatTan[detRegion][2])/vms/vid;
00181     }
00182     
00183   }
00184   
00185   
00186   double pTMS() {return ptms*(1+smearMS);} /* smeared muon spectrometer pT */
00187   double pTMS(double SmearMS) {return ptms*(1+SmearMS);} /* smeared muon spectrometer pT */
00188   
00189   double pTID() {return ptid*(1+smearID);}  /* smeared inner detector pT */
00190   double pTID(double SmearID) {return ptid*(1+SmearID);}  /* smeared inner detector pT */
00191   
00192   double pTCB() {
00193     if (useScale && detRegion>-1 && detRegion<4) return scale_CB[detRegion]*ptcb*(1+smearCB);
00194     else return ptcb*(1+smearCB);
00195   }  /* smeared combined pT*/
00196 
00197   double pTCB(double SmearCB) {
00198     if (useScale && detRegion>-1 && detRegion<4) return scale_CB[detRegion]*ptcb*(1+SmearCB);
00199     else return ptcb*(1+SmearCB);
00200   }  /* smeared combined pT*/
00201   
00202   double pTCBScaleUp() {
00203     return (scale_CB[detRegion]+scaleSyst_CB[detRegion])*ptcb*(1+smearCB);
00204   }
00205 
00206   double pTCBScaleDown() {
00207     return (scale_CB[detRegion]-scaleSyst_CB[detRegion])*ptcb*(1+smearCB);
00208   }
00209 
00210   
00211   double SMS() {return smearMS;}
00212   double SID() {return smearID;}
00213   double SCB() {return smearCB;}
00214   
00215   double VMS()  {return vms;}  /* error smearMS */
00216   double VID()  {return vid;}  /* error smearID */
00217   double Corr() {return corr;} /* correlation between smearMS and smearID */
00218   
00219   
00220   void MSUP(double &PtMS) {
00221     double SmearMS=smearMS+Sign(smearMS)*vms;
00222     PtMS=pTMS(SmearMS);
00223   }
00224   
00225   void MSUP(double &PtMS,double &PtID,double &PtCB) {
00226     double SmearMS=smearMS+Sign(smearMS)*vms;
00227     double SmearID=smearID+Sign(smearMS)*vid*corr;
00228     
00229     PtMS=pTMS(SmearMS);
00230     PtID=pTID(SmearID);
00231     if (detType=="All") {
00232       double SmearCB=Combine(SmearMS,SmearID);
00233       PtCB=pTCB(SmearCB);
00234     }
00235   }
00236   
00237   void MSLOW(double &PtMS) {
00238     double SmearMS=smearMS-Sign(smearMS)*vms;
00239     PtMS=pTMS(SmearMS);
00240   }
00241   
00242   void MSLOW(double &PtMS,double &PtID,double &PtCB) {
00243     double SmearMS=smearMS-Sign(smearMS)*vms;
00244     double SmearID=smearID-Sign(smearMS)*vid*corr;
00245     
00246     PtMS=pTMS(SmearMS);
00247     PtID=pTID(SmearID);
00248     if (detType=="All") {
00249       double SmearCB=Combine(SmearMS,SmearID);
00250       PtCB=pTCB(SmearCB);
00251     }
00252   }
00253 
00254 
00255   void IDUP(double &PtID) {
00256     double SmearID=smearID+Sign(smearID)*vid;
00257     PtID=pTID(SmearID);
00258   }
00259   
00260   void IDUP(double &PtMS,double &PtID,double &PtCB) {
00261     double SmearMS=smearMS+Sign(smearID)*vms*corr;
00262     double SmearID=smearID+Sign(smearID)*vid;
00263     
00264     PtMS=pTMS(SmearMS);
00265     PtID=pTID(SmearID);
00266     if (detType=="All") {
00267       double SmearCB=Combine(SmearMS,SmearID);
00268       PtCB=pTCB(SmearCB);
00269     }
00270   }
00271 
00272   void IDLOW(double &PtID) {
00273     double SmearID=smearID-Sign(smearID)*vid;
00274     PtID=pTID(SmearID);
00275   }
00276   
00277   void IDLOW(double &PtMS,double &PtID,double &PtCB) {
00278     double SmearMS=smearMS-Sign(smearID)*vms*corr;
00279     double SmearID=smearID-Sign(smearID)*vid;
00280     
00281     PtMS=pTMS(SmearMS);
00282     PtID=pTID(SmearID);
00283     if (detType=="All") {
00284       double SmearCB=Combine(SmearMS,SmearID);
00285       PtCB=pTCB(SmearCB);
00286     }
00287   }
00288   
00289 
00290 
00291   void PTVar(double &Pt,std::string fun) {
00292     if (fun=="IDUP") IDUP(Pt);
00293     else if (fun=="IDLOW") IDLOW(Pt);
00294     else if (fun=="MSUP") MSUP(Pt);
00295     else if (fun=="MSLOW") MSLOW(Pt);
00296   }
00297   
00298   void PTVar(double &PtMS,double &PtID,double &PtCB,std::string fun) {
00299     if (fun=="IDUP") IDUP(PtMS,PtID,PtCB);
00300     else if (fun=="IDLOW") IDLOW(PtMS,PtID,PtCB);
00301     else if (fun=="MSUP") MSUP(PtMS,PtID,PtCB);
00302     else if (fun=="MSLOW") MSLOW(PtMS,PtID,PtCB);
00303   }
00304   
00305 
00306   double Sign(double x){return (x<0? -1 : 1);}
00307 
00308 
00311   double ptms_orig() {return ptms;}
00312   double ptid_orig() {return ptid;}
00313   double ptcb_orig() {return ptcb;}
00314   double Eta() {return eta;}  
00315   int DetRegion() {return detRegion;}
00316 
00317 
00318   /*members*/
00319   TRandom3* gRand;
00320   double ptms,ptid,ptcb,eta;
00321   double vms,vid,corr;
00322   double smearMS,smearID,smearCB;
00323   bool   useTan2;
00324   std::string detType;
00325   int detRegion;
00326   double GeV;
00327   double g1,g2,g3,g4;
00328   double wMS,wID;
00329   bool useScale;
00330 
00331   static const double scale_CB[4];   
00332   static const double scaleSyst_CB[4];  
00333 
00334   static const double ID_MS[4];   
00335   static const double ID_AL[4];   
00336   static const double ID_ALTAN[4];
00337   static const double MS_MS[4];   
00338   static const double MS_AL[4];   
00339 
00340   static const double E_ID_MS[4];   
00341   static const double E_ID_AL[4];   
00342   static const double E_ID_ALTAN[4];
00343   static const double E_MS_MS[4];   
00344   static const double E_MS_AL[4];   
00345   
00346   static const double S_ID_MS[4];   
00347   static const double S_ID_AL[4];   
00348   static const double S_ID_ALTAN[4];
00349   static const double S_MS_MS[4];   
00350   static const double S_MS_AL[4];   
00351 
00352   static const double MC_ID_MS[4];   
00353   static const double MC_ID_AL[4];   
00354   static const double MC_ID_ALTAN[4];
00355   static const double MC_MS_CALO[4]; 
00356   static const double MC_MS_MS[4];   
00357   static const double MC_MS_AL[4];   
00358   
00359   static const double CorrMat[4][3];
00360   static const double CorrMatTan[4][3];
00361   
00362 };
00363 
00364 
00366 const double SmearingClass::scale_CB[4] = {1.0000, 0.9933, 1.0037, 1.0033};
00367 const double SmearingClass::scaleSyst_CB[4] = {0.0006 , 0.0013, 0.0021, 0.0022};
00368 
00370 const double SmearingClass::ID_MS[4]       = {0,0,0,0};
00371 const double SmearingClass::ID_AL[4]       = {0.000255,0.00053,0.000819,0.001177};
00372 const double SmearingClass::ID_ALTAN[4]    = {0,0,0,0.0000593};
00373 const double SmearingClass::MS_MS[4]       = {0.02428,0.0678,0.0362,0.0237};
00374 const double SmearingClass::MS_AL[4]       = {0.000191,0.00019,0.00012,0.000669};
00376 const double SmearingClass::E_ID_MS[4]       = {0,0,0,0};
00377 const double SmearingClass::E_ID_AL[4]       = {0.000037,0.00013,0.000071,0.000050};
00378 const double SmearingClass::E_ID_ALTAN[4]    = {0,0,0,0.0000032};
00379 const double SmearingClass::E_MS_MS[4]       = {0.00088,0.0027,0.0019,0.0045};
00380 const double SmearingClass::E_MS_AL[4]       = {0.000028,0.00012,0.00011,0.000091};
00382 const double SmearingClass::S_ID_MS[4]       = {0,0,0,0};
00383 const double SmearingClass::S_ID_AL[4]       = {0.000015,0.00014,0.000144,0.00014};
00384 const double SmearingClass::S_ID_ALTAN[4]    = {0,0,0,0.0000072};
00385 const double SmearingClass::S_MS_MS[4]       = {0.00019,0.0016,0.0028,0.0108};
00386 const double SmearingClass::S_MS_AL[4]       = {0.000009,0.00004,0.00008,0.000167};
00387 
00389 const double SmearingClass::MC_ID_MS[4]    = {0.016,0.026,0.034,0.048};
00390 const double SmearingClass::MC_ID_AL[4]    = {0.00029,0.0003,0.00042,0.00094};
00391 const double SmearingClass::MC_ID_ALTAN[4] = {0,0,0,0.000063};
00392 const double SmearingClass::MC_MS_CALO[4]  = {0.23,0,0,0.17};
00393 const double SmearingClass::MC_MS_MS[4]    = {0.027,0.054,0.033,0.027};
00394 const double SmearingClass::MC_MS_AL[4]    = {0.00013,0.00023,0.000093,0.000074};
00395 /***********************************/
00397 const double SmearingClass::CorrMat[4][3] = { {0.248,-0.878,-0.493}, 
00398                                               {-0.698,-0.157,-0.200},
00399                                               {-0.299,-0.732,-0.024}, 
00400                                               {-0.818,-0.421,-0.051}};
00401 
00402 const double SmearingClass::CorrMatTan[4][3] = { {-0.937,-0.095,0.077}, 
00403                                                  {-0.650,-0.767,0.47},
00404                                                  {0.077,-0.553,-0.345}, 
00405                                                  {-0.326,-0.301,-0.375}};
00406 
00407
/afs/hep.man.ac.uk/u/markowen/ATLAS/SFrameProof/2011EPS/topUtils/src/SmearingClass_rel16.h
