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

00001 
00002 // SmearingClass.h -- ATLAS Experiment Software //2
00014 
00015 #include <cstring>
00016 #include <iostream>
00017 #include <stdio.h>
00018 #include <TROOT.h>
00019 #include <math.h>
00020 #include "TRandom3.h"
00021 
00022 #ifndef MuonMomentumCorrections_SmearingClass
00023 #define MuonMomentumCorrections_SmearingClass
00024 
00025 class SmearingClass{
00026  public:
00027   /*Constructors*/
00028   SmearingClass() {Initialize();}
00029   
00030   
00031   /************/
00032   /* Methods */    
00033   /**********/
00034   
00036   void Initialize() {    
00037     detType="";
00038     GeV=1000;
00039     gRand = new TRandom3();
00040     gRand->SetSeed(0);
00041     useTan2=true;
00042     useScale=false;
00043     ptms=0;
00044     ptid=0;
00045     ptcb=0;
00046     eta=0;
00047   }
00048   
00051   void SetSeed(int seed) {gRand->SetSeed(seed);}
00052   void SetSeed(int evtnum, int muon_index, int offset=680049) {gRand->SetSeed(offset + evtnum + muon_index*100);}
00053   void UseTan2(bool U) {useTan2=U;}/* use pT/tan(theta)^2 parameterization for CSC region */  
00054   void UseGeV() {GeV=1;}  /* if input momenta are in GeV */
00055   void UseScale(bool applyScale) {useScale=applyScale;}  /* also use momentum scale, set as default */
00056 
00057   
00058   /* MS or ID smearing only. DetType="MS","ID"*/
00059   void Event(double Pt, double Eta,std::string DetType) 
00060   {        
00061     ptms=0;
00062     ptid=0;
00063     ptcb=0;
00064     
00065     if (DetType=="MS")      ptms = Pt;
00066     else if (DetType=="ID") ptid = Pt;
00067     else std::cout<<"SmearingClass ERROR: wrong DetType in input "<<DetType<<std::endl; 
00068     eta=Eta;    
00069     detType=DetType;
00070     Event();
00071   }
00072   
00073   /* For full MS ID and CB smearing */
00074   void Event(double PtMS, double PtID, double PtCB, double Eta) 
00075   {
00076     detType="All";
00077     ptms=PtMS;
00078     ptid=PtID;
00079     ptcb=PtCB;
00080     eta=Eta;
00081     Event();
00082   }
00083 
00084   /* CB smearing only */
00085   void Event(double PtCB, double Eta)
00086   {
00087     detType="All";
00088     ptms=PtCB;
00089     ptid=PtCB;
00090     ptcb=PtCB;
00091     eta=Eta;
00092     Event();
00093   }
00094   
00095   void  Event()  {
00096     smearMS=0;
00097     smearID=0;
00098     smearCB=0;
00099     
00100     /* Detector Region */
00101     if (fabs(eta)<1.05)     detRegion=0;
00102     else if (fabs(eta)<1.7) detRegion=1;
00103     else if (fabs(eta)<2.0) detRegion=2;
00104     else if (fabs(eta)<2.5) detRegion=3;
00105     else detRegion=-1;
00106     
00107     g1 = gRand->Gaus(0,1);  
00108     g2 = gRand->Gaus(0,1); 
00109     g3 = gRand->Gaus(0,1); 
00110     g4 = gRand->Gaus(0,1); 
00111     
00113     smearMS=Smearing("MS");
00114     smearID=Smearing("ID");
00115     if (detType=="All") smearCB=Combine(smearMS,smearID);
00116     ErrorMatrix();    
00117   }
00118   
00119 
00120   double Smearing(std::string DetType) {
00121     
00122     if (detRegion<0||detRegion>3) return 0;
00123     
00124     if (DetType=="MS") {return (MS_MS[detRegion]*g1+MS_AL[detRegion]*g2*ptms/GeV);}
00125     else if(DetType=="ID") 
00126       {
00127         if (useTan2&&detRegion==3) {return (ID_MS[detRegion]*g3+ID_ALTAN[detRegion]*g4*ptid/GeV*sinh(eta)*sinh(eta));}
00128         else {return (ID_MS[detRegion]*g3+ID_AL[detRegion]*g4*ptid/GeV);}
00129       } 
00130     else std::cout<<"Error:: DetType not defined "<<DetType<<std::endl;
00131     return 0;
00132   }
00133   
00134   double Combine(double SmearMS,double SmearID) {
00137     if (detRegion<0||detRegion>3) return 0;
00138     if (ptcb==0) {std::cout<<"Warning: ptcb==0"<<std::endl;return 0;}
00139     
00140     double SigmaMS=pow(
00141                        pow(MC_MS_CALO[detRegion]/ptcb*GeV,2)+
00142                        pow(MC_MS_MS[detRegion],2)+
00143                        pow(MC_MS_AL[detRegion]*ptcb/GeV,2),0.5);
00144     
00145     double SigmaID=pow(
00146                        pow(MC_ID_MS[detRegion],2)+
00147                        pow(MC_ID_AL[detRegion]*ptcb/GeV,2),0.5);
00148     
00149     if (detRegion==3&&useTan2) {
00150       SigmaID=pow(
00151                   pow(MC_ID_MS[detRegion],2)+
00152                   pow(MC_ID_ALTAN[detRegion]*ptcb/GeV*sinh(eta)*sinh(eta),2),0.5);
00153     }
00154     
00155     
00156     wMS=1./SigmaMS/SigmaMS;
00157     wID=1./SigmaID/SigmaID;
00158     return (SmearMS*wMS+SmearID*wID)/(wMS+wID);
00159   }
00160   
00161   void ErrorMatrix() {
00162     
00163     vms=0;
00164     vid=0;
00165     corr=0;
00166     
00167     if (detRegion<0||detRegion>3) return;
00168     
00169     
00170     if (!useTan2||detRegion!=3){
00171       double s1=pow(E_MS_MS[detRegion]*E_MS_MS[detRegion]+S_MS_MS[detRegion]*S_MS_MS[detRegion],0.5);
00172       double s2=pow(E_MS_AL[detRegion]*E_MS_AL[detRegion]+S_MS_AL[detRegion]*S_MS_AL[detRegion],0.5);
00173       double s3=pow(E_ID_MS[detRegion]*E_ID_MS[detRegion]+S_ID_MS[detRegion]*S_ID_MS[detRegion],0.5);
00174       double s4=pow(E_ID_AL[detRegion]*E_ID_AL[detRegion]+S_ID_AL[detRegion]*S_ID_AL[detRegion],0.5);
00175       
00176       vms=g1*g1*s1*s1+
00177         g2*g2*s2*s2*ptms/GeV*ptms/GeV+
00178         2.0*g1*g2*ptms/GeV*s1*s2*CorrMat[detRegion][0];
00179       vms=pow(vms,0.5);
00180       
00181       vid=g3*g3*s3*s3+g4*g4*s4*s4*ptid/GeV*ptid/GeV+2.0*g3*g4*s3*s4*ptid/GeV*CorrMat[detRegion][5];
00182       vid=pow(vid,0.5);
00183       
00184       if (vms*vid!=0)
00185         corr=(g1*s1*(g3*s3*CorrMat[detRegion][1]+g4*s4*ptid/GeV*CorrMat[detRegion][2])+
00186               g2*s2*ptms/GeV*(g3*s3*CorrMat[detRegion][3]+g4*s4*ptid/GeV*CorrMat[detRegion][4]))/vms/vid;
00187       
00188     } else {
00189       double s1=pow(E_MS_MS[detRegion]*E_MS_MS[detRegion]+S_MS_MS[detRegion]*S_MS_MS[detRegion],0.5);
00190       double s2=pow(E_MS_AL[detRegion]*E_MS_AL[detRegion]+S_MS_AL[detRegion]*S_MS_AL[detRegion],0.5);
00191       double s3=pow(E_ID_MS[detRegion]*E_ID_MS[detRegion]+S_ID_MS[detRegion]*S_ID_MS[detRegion],0.5);
00192       double s4=pow(E_ID_ALTAN[detRegion]*E_ID_ALTAN[detRegion]+S_ID_ALTAN[detRegion]*S_ID_ALTAN[detRegion],0.5);
00193       
00194       vms=g1*g1*s1*s1+
00195         g2*g2*s2*s2*ptms/GeV*ptms/GeV+
00196         2.0*g1*g2*ptms/GeV*s1*s2*CorrMatTan[detRegion][0];
00197       vms=pow(vms,0.5);
00198       
00199       vid=g3*g3*s3*s3+g4*g4*s4*s4*ptid/GeV*ptid/GeV*sinh(eta)*sinh(eta)*sinh(eta)*sinh(eta)
00200         +2.0*g3*g4*s3*s4*ptid/GeV*sinh(eta)*sinh(eta)*CorrMatTan[detRegion][5];
00201       vid=pow(vid,0.5);
00202       
00203       if (vms*vid!=0)
00204         corr=(g1*s1*(g3*s3*CorrMatTan[detRegion][1]+g4*s4*ptid/GeV*sinh(eta)*sinh(eta)*CorrMatTan[detRegion][2])+
00205               g2*s2*ptms/GeV*(g3*s3*CorrMatTan[detRegion][3]+g4*s4*ptid/GeV*sinh(eta)*sinh(eta)*CorrMatTan[detRegion][4]))/vms/vid;
00206       
00207     }
00208     
00209   }
00210 
00211   
00212   double pTMS() {return ptms*(1+smearMS);} /* smeared muon spectrometer pT */
00213   double pTMS(double SmearMS) {return ptms*(1+SmearMS);} /* smeared muon spectrometer pT */
00214   
00215   double pTID() {return ptid*(1+smearID);}  /* smeared inner detector pT */
00216   double pTID(double SmearID) {return ptid*(1+SmearID);}  /* smeared inner detector pT */
00217   
00218   double pTCB() {
00219     if (useScale && detRegion>-1 && detRegion<4) return scale_CB[detRegion]*ptcb*(1+smearCB);
00220     else return ptcb*(1+smearCB);
00221   }  /* smeared combined pT*/
00222 
00223   double pTCB(double SmearCB) {
00224     if (useScale && detRegion>-1 && detRegion<4) return scale_CB[detRegion]*ptcb*(1+SmearCB);
00225     else return ptcb*(1+SmearCB);
00226   }  /* smeared combined pT*/
00227   
00228   
00229   double pTCBScaleUp() {
00230     return (scale_CB[detRegion]+scaleSyst_CB[detRegion])*ptcb*(1+smearCB);
00231   }
00232   
00233   double pTCBScaleDown() {
00234     return (scale_CB[detRegion]-scaleSyst_CB[detRegion])*ptcb*(1+smearCB);
00235   }
00236 
00237   double SMS() {return smearMS;}
00238   double SID() {return smearID;}
00239   double SCB() {return smearCB;}
00240   
00241   double VMS()  {return vms;}  /* error smearMS */
00242   double VID()  {return vid;}  /* error smearID */
00243   double Corr() {return corr;} /* correlation between smearMS and smearID */
00244   
00245   
00246   void MSUP(double &PtMS) {
00247     double SmearMS=smearMS+Sign(smearMS)*vms;
00248     PtMS=pTMS(SmearMS);
00249   }
00250   
00251   void MSUP(double &PtMS,double &PtID,double &PtCB) {
00252     double SmearMS=smearMS+Sign(smearMS)*vms;
00253     double SmearID=smearID+Sign(smearMS)*vid*corr;
00254     
00255     PtMS=pTMS(SmearMS);
00256     PtID=pTID(SmearID);
00257     if (detType=="All") {
00258       double SmearCB=Combine(SmearMS,SmearID);
00259       PtCB=pTCB(SmearCB);
00260     }
00261   }
00262   
00263   void MSLOW(double &PtMS) {
00264     double SmearMS=smearMS-Sign(smearMS)*vms;
00265     PtMS=pTMS(SmearMS);
00266   }
00267   
00268   void MSLOW(double &PtMS,double &PtID,double &PtCB) {
00269     double SmearMS=smearMS-Sign(smearMS)*vms;
00270     double SmearID=smearID-Sign(smearMS)*vid*corr;
00271     
00272     PtMS=pTMS(SmearMS);
00273     PtID=pTID(SmearID);
00274     if (detType=="All") {
00275       double SmearCB=Combine(SmearMS,SmearID);
00276       PtCB=pTCB(SmearCB);
00277     }
00278   }
00279 
00280 
00281   void IDUP(double &PtID) {
00282     double SmearID=smearID+Sign(smearID)*vid;
00283     PtID=pTID(SmearID);
00284   }
00285   
00286   void IDUP(double &PtMS,double &PtID,double &PtCB) {
00287     double SmearMS=smearMS+Sign(smearID)*vms*corr;
00288     double SmearID=smearID+Sign(smearID)*vid;
00289     
00290     PtMS=pTMS(SmearMS);
00291     PtID=pTID(SmearID);
00292     if (detType=="All") {
00293       double SmearCB=Combine(SmearMS,SmearID);
00294       PtCB=pTCB(SmearCB);
00295     }
00296   }
00297 
00298   void IDLOW(double &PtID) {
00299     double SmearID=smearID-Sign(smearID)*vid;
00300     PtID=pTID(SmearID);
00301   }
00302   
00303   void IDLOW(double &PtMS,double &PtID,double &PtCB) {
00304     double SmearMS=smearMS-Sign(smearID)*vms*corr;
00305     double SmearID=smearID-Sign(smearID)*vid;
00306     
00307     PtMS=pTMS(SmearMS);
00308     PtID=pTID(SmearID);
00309     if (detType=="All") {
00310       double SmearCB=Combine(SmearMS,SmearID);
00311       PtCB=pTCB(SmearCB);
00312     }
00313   }
00314   
00315 
00316 
00317   void PTVar(double &Pt,std::string fun) {
00318     if (fun=="IDUP") IDUP(Pt);
00319     else if (fun=="IDLOW") IDLOW(Pt);
00320     else if (fun=="MSUP") MSUP(Pt);
00321     else if (fun=="MSLOW") MSLOW(Pt);
00322   }
00323   
00324   void PTVar(double &PtMS,double &PtID,double &PtCB,std::string fun) {
00325     if (fun=="IDUP") IDUP(PtMS,PtID,PtCB);
00326     else if (fun=="IDLOW") IDLOW(PtMS,PtID,PtCB);
00327     else if (fun=="MSUP") MSUP(PtMS,PtID,PtCB);
00328     else if (fun=="MSLOW") MSLOW(PtMS,PtID,PtCB);
00329   }
00330   
00331 
00332   double Sign(double x){return (x<0? -1 : 1);}
00333 
00334 
00337   double ptms_orig() {return ptms;}
00338   double ptid_orig() {return ptid;}
00339   double ptcb_orig() {return ptcb;}
00340   double Eta() {return eta;}  
00341   int DetRegion() {return detRegion;}
00342 
00343 
00344   /*members*/
00345   TRandom3* gRand;
00346   double ptms,ptid,ptcb,eta;
00347   double vms,vid,corr;
00348   double smearMS,smearID,smearCB;
00349   bool   useTan2;
00350   std::string detType;
00351   int detRegion;
00352   double GeV;
00353   double g1,g2,g3,g4;
00354   double wMS,wID;
00355   bool useScale;
00356 
00357   static const double scale_CB[4];   
00358   static const double scaleSyst_CB[4];  
00359 
00360   static const double ID_MS[4];   
00361   static const double ID_AL[4];   
00362   static const double ID_ALTAN[4];
00363   static const double MS_MS[4];   
00364   static const double MS_AL[4];   
00365 
00366   static const double E_ID_MS[4];   
00367   static const double E_ID_AL[4];   
00368   static const double E_ID_ALTAN[4];
00369   static const double E_MS_MS[4];   
00370   static const double E_MS_AL[4];   
00371   
00372   static const double S_ID_MS[4];   
00373   static const double S_ID_AL[4];   
00374   static const double S_ID_ALTAN[4];
00375   static const double S_MS_MS[4];   
00376   static const double S_MS_AL[4];   
00377 
00378   static const double MC_ID_MS[4];   
00379   static const double MC_ID_AL[4];   
00380   static const double MC_ID_ALTAN[4];
00381   static const double MC_MS_CALO[4]; 
00382   static const double MC_MS_MS[4];   
00383   static const double MC_MS_AL[4];   
00384   
00385   static const double CorrMat[4][6];
00386   static const double CorrMatTan[4][6];
00387   
00388 };
00389 
00390 
00392 const double SmearingClass::scale_CB[4] = {0.9997, 0.9999, 0.9990, 1.0013};
00393 const double SmearingClass::scaleSyst_CB[4] = {0.0002, 0.0006, 0.0012, 0.0007};
00394 
00396 const double SmearingClass::ID_MS[4]       = {0.00624,0.00006,0.00087,0.00087};
00397 const double SmearingClass::ID_AL[4]       = {0.000299,0.000721,0.000845,0.};
00398 const double SmearingClass::ID_ALTAN[4]    = {0,0,0,0.000048};
00399 const double SmearingClass::MS_MS[4]       = {0.02035,0.04994,0.02643,0.01705};
00400 const double SmearingClass::MS_AL[4]       = {0.000129,0.000335,0.000163,0.000443};
00402 const double SmearingClass::E_ID_MS[4]       = {0.00102,0.00438,0.00443,0.00415};
00403 const double SmearingClass::E_ID_AL[4]       = {0.000015,0.000049,0.000016,0.};
00404 const double SmearingClass::E_ID_ALTAN[4]    = {0,0,0,0.000003};
00405 const double SmearingClass::E_MS_MS[4]       = {0.00016,0.00198,0.00045,0.00383};
00406 const double SmearingClass::E_MS_AL[4]       = {0.000023,0.000012,0.000041,0.000053};
00408 const double SmearingClass::S_ID_MS[4]       = {0,0.00001,0,0};
00409 const double SmearingClass::S_ID_AL[4]       = {0.000018,0.000063,0.000017,0.};
00410 const double SmearingClass::S_ID_ALTAN[4]    = {0,0,0,0.000008};
00411 const double SmearingClass::S_MS_MS[4]       = {0.00071,0.00105,0.00526,0.00471};
00412 const double SmearingClass::S_MS_AL[4]       = {0.000011,0.000082,0.00003,0.000026};
00413 
00415 const double SmearingClass::MC_ID_MS[4]    = {0.01549,0.02524,0.03385,0.04623};
00416 const double SmearingClass::MC_ID_AL[4]    = {0.000303,0.000333,0.000431,0};
00417 const double SmearingClass::MC_ID_ALTAN[4] = {0,0,0,0.000051};
00418 const double SmearingClass::MC_MS_CALO[4]  = {0.26,0,0,0.17};
00419 const double SmearingClass::MC_MS_MS[4]    = {0.02557,0.04934,0.0383,0.0271};
00420 const double SmearingClass::MC_MS_AL[4]    = {0.0001144,0.000202,0.000084,0.000092};
00421 
00422 
00423 /***********************************/
00425 const double SmearingClass::CorrMat[4][6] = { {-0.811,-0.608, 0.932, 0.111,-0.875,-0.559}, 
00426                                               {-0.088,-0.347,-0.543, 0.078,-0.040, 0.366}, 
00427                                               { 1.217,-0.251,-1.006,-0.190,-0.884, 0.211},
00428                                               {-0.537,-0.281,-0.349, 0.072, 0.008,-0.412}};
00429 
00430 
00432 const double SmearingClass::CorrMatTan[4][6] = { {-0.811,-0.608, 0.932, 0.111,-0.875,-0.559}, 
00433                                                  {-0.088,-0.347,-0.543, 0.078,-0.040, 0.366}, 
00434                                                  { 1.217,-0.251,-1.006,-0.190,-0.884, 0.211},
00435                                                  {-0.537,-0.281,-0.349, 0.072, 0.008,-0.412}};
00436 
00437 #endif