SFHMarana: TMarEvent.C Source File

00001 
00002 // Class TMarEvent
00003 //
00004 // Author     : Emmanuel Sauvan/Matti Peez/F. Cassol Brunner
00005 // Created    : 02/2002
00006 // Last update: $Date: 2005/06/30 17:08:00 $
00007 //          by: $Author: sauvan $
00008 // Comment: Contains definition of event class
00010 
00011 #ifndef __TMAREVENT
00012 #include "Marana/TMarEvent.h"
00013 #endif
00014 
00015 #define TMAREVENTDEBUG 0
00016 
00017 
00018 
00019 
00020 // Import class
00021 //ClassImp(TMarEvent)
00022 
00023 TMarEvent::TMarEvent(const TRun *run)
00024 {
00025 
00026 //... check the if the DataType and DataYear are corrected
00027 
00028   if(run->Type.find("Data")!=string::npos) RunType=Data;
00029   else if(run->Type.find("EE_Grape")!=string::npos)RunType=EE_Grape;
00030   else if(run->Type.find("MM_Grape")!=string::npos)RunType=MM_Grape;
00031   else if(run->Type.find("TT_Grape")!=string::npos)RunType=TT_Grape;
00032   else if(run->Type.find("NC_Rap")!=string::npos)RunType=NC_Rap;
00033   else if(run->Type.find("NC_Djo")!=string::npos)RunType=NC_Djo; 
00034   else if(run->Type.find("CC_Djo")!=string::npos)RunType=CC_Djo; 
00035   else if(run->Type.find("EG_Wabgen")!=string::npos)RunType=EG_Wabgen; 
00036   else if(run->Type.find("W_Epvec")!=string::npos)RunType=W_Epvec;
00037   else if(run->Type.find("WH_EpvecH")!=string::npos)RunType=W_Epvec;
00038   else if(run->Type.find("WC_Epvec")!=string::npos)RunType=W_Epvec;
00039   else if(run->Type.find("GP_Pythia")!=string::npos)RunType=GP_Pythia;
00040   else if(run->Type.find("QCDINS")!=string::npos)RunType=QCD_Ins;
00041   else if(run->Type.find("Htm")!=string::npos)RunType=Htm; //--- H++ MC
00042   else if(run->Type.find("Anotop")!=string::npos)RunType=Anotop; //--- single top, anotop
00043   else if(run->Type.find("PsCC")!=string::npos)RunType=PsCC; //--- pseudo CC
00044   else if(run->Type.find("POUDS_Rap")!=string::npos)RunType=POUDS; //--- pomeron proton dissocie uds
00045   else if(run->Type.find("RAPD_UDS")!=string::npos)RunType=RAPD_UDS; //--- pomeron proton dissocie uds
00046   else {
00047     RunType=UNKNOWN;
00048     cout << "TMarEvent::TMarEvent : Your DataType ("<<run->Type
00049          << ") is unknown. " << endl;
00050   }
00051   if(run->Year == "9497") RunYear= Y9497;
00052   else if(run->Year == "9899") RunYear= Y9899;
00053   else if(run->Year == "9900") RunYear= Y9900;  
00054   else if(run->Year == "0203") RunYear= Y0203;  
00055   else if(run->Year == "0304") RunYear= Y0304;  
00056   else if(run->Year == "0304R") RunYear= Y0304; //--- for separation of left and right polar periods
00057   else if(run->Year == "0304L") RunYear= Y0304;  
00058   else if(run->Year == "0304L1") RunYear= Y0304;  
00059   else if(run->Year == "0304L2") RunYear= Y0304;  
00060   else if(run->Year == "0304R2") RunYear= Y0304;  
00061   else if(run->Year == "05") RunYear= Y05;   //--- e- 2005 (if some one day ...)
00062   else if(run->Year == "05L") RunYear= Y05;   //--- e- 2005 (if some one day ...)
00063   else {
00064     RunYear= Y9900;  
00065     cout << "TMarEvent::TMarEvent : Your DataYear ("<<run->Year
00066          << ") is unknown. We put 9900 as default." << endl;
00067    
00068   }   
00069   numMax=run->nMaxEvt;
00070   
00071 // pointers to braches 
00072   static H1PartMuonArrayPtr ModsPartMuon("Muons"); //pointers on muODS particles
00073   static H1PartEmArrayPtr ModsPartEm("EmParticles"); //pointers on muODS particles 
00074   static H1PartJetArrayPtr ModsPartJet("KtJets");  
00075   static H1PartCandArrayPtr ModsPartCand("Particles");  // all candidats particl
00076   static H1PartSelTrackArrayPtr ModsPartSelTrack("SelectedTracks");        // 
00077 //  static H1ExclHfsIterator ModsPartHFS;
00078 
00079   static H1PartMCArrayPtr ModsPartMC("MCParticles");
00080 
00081 #if USERTREE
00082   static H1PartJetArrayPtr ModsGenKtJets("GenKtJets");
00083   static H1PartSelTrackArrayPtr ModsDTRATracks("DTRATracks");
00084   static H1PartSelTrackArrayPtr ModsDTNVTracks("DTNVTracks");
00085 #else
00086   static H1PartJetArrayPtr ModsGenKtJets("KtJets");
00087   static H1PartSelTrackArrayPtr ModsDTRATracks("SelectedTracks");
00088   static H1PartSelTrackArrayPtr ModsDTNVTracks("SelectedTracks");
00089 #endif
00090 
00091   
00092 //... if necessary  inizialise the lumi and MC cut on generated variables
00093 
00094   if(run->lumi_file!="not_given" )
00095            
00096            lumi= new TMarLumi(run->lumi_file.c_str(),run->lumi_file_bad_runs.c_str());
00097   else lumi=NULL;
00098   
00099 //... init the  genCut list   
00100   genCut =run->genCut;
00101 
00102 //... init the  lumi weight
00103   if(run->Lumi!=0)LumiWeight=1/run->Lumi;
00104   else throw string("TMarEvent::TMarEvent = you inizialize LumiWeight with Lumi=0"); 
00105 
00106 //... init the hfs calibration class    
00107   if(run->HfsCalib==1) {
00108     JetCalib = new H1JetCalibration;
00109     JetCalib->InitHadroo2KtJetCalibration((int)RunYear,(int)RunType);
00110     HfsCalib=NULL;
00111   }else if(run->HfsCalib==2) {
00112     HfsCalib = new H1HfsCalibEmHad;
00113     HfsCalib->InitHfsCalibration((int)RunYear,(int)RunType);
00114     JetCalib=NULL;
00115   } else {
00116     JetCalib=NULL;
00117     HfsCalib=NULL;
00118   }
00119   
00120 //... Init Electron calibration, only possible for new data at the moment
00121   if(run->ElCalib==1  && (RunYear==Y0203 || RunYear==Y0304 || RunYear==Y05) ) {
00122     printf(">>>>>  Init TMarH1ElecCalibration for Hera2  <<<<<\n");
00123     Bool_t bDoStackWiseCalibration = kTRUE;
00124     Bool_t bDoZimpactWiseCalibration = kTRUE;
00125     Bool_t bDoResolutionSmearing = kTRUE;
00126     H1ElecCalib = NULL;
00127     ElecCalib =  new TMarH1ElecCalibration(
00128                         bDoStackWiseCalibration, bDoZimpactWiseCalibration, bDoResolutionSmearing);
00129   }else if(run->ElCalib==2 && (RunYear==Y0203 || RunYear==Y0304 || RunYear==Y05) ) {//--- use Local H1ElecCalibration
00130     printf(">>>>>  Init H1ElecCalibration for Hera2  <<<<<\n");
00131     Bool_t bDoStackWiseCalibration = kTRUE;
00132     Bool_t bDoZimpactWiseCalibration = kTRUE;
00133     Bool_t bDoResolutionSmearing = kTRUE;
00134     ElecCalib = NULL;
00135     H1ElecCalib =  new H1ElecCalibration(
00136                         bDoStackWiseCalibration, bDoZimpactWiseCalibration, bDoResolutionSmearing);
00137   } else {
00138     ElecCalib=NULL;
00139     H1ElecCalib=NULL;
00140   }
00141 
00142 //... init trigger reweight class
00143   fTrigReweight = new TTriggerReweight();
00144   fMCReweight = new TMCXSectionReweight();
00145 
00146 //... Initialisation des listes de particules
00147   
00148   electron = new TClonesArray("TMarBody",1);
00149   photon = new TClonesArray("TMarBody",1);
00150   muon = new TClonesArray("TMarBody",1);
00151   jet = new TClonesArray("TMarBody",1);
00152   np = new TClonesArray("TMarBody",1);
00153 
00154 
00155 
00156   GenElectron = new TClonesArray("TMarBody",1);
00157   GenPhoton = new TClonesArray("TMarBody",1);
00158   GenMuon = new TClonesArray("TMarBody",1);
00159   GenTau = new TClonesArray("TMarBody",1);
00160   GenNp = new TClonesArray("TMarBody",1);
00161   GenHfs = new TClonesArray("TMarBody",1);
00162 
00163 
00164 //--- enable thight finders by default
00165   find_ElePho=true;
00166   find_Muon=true;
00167   find_Jet=true;
00168   find_NP=true;
00169 
00170 //--- default cuts for thight particles finders
00171   PtMinEle=5;
00172   ThetaMinEle=20;
00173   ThetaMaxEle=170;
00174   ConeEle=0.5;
00175 
00176   PtMinMu=0.0;
00177   ThetaMinMu=0;
00178   ThetaMaxMu=180;
00179   DTrTrMinMu=0.5;
00180 
00181   PtMinNp=5;
00182   ThetaMinNp=20;
00183   ThetaMaxNp=150;
00184 
00185   PtMinJet=2;
00186   ThetaMinJet=5;
00187   ThetaMaxJet=170;
00188   RMinJet=0.1;
00189   EmFracJet=0.9;
00190  fQCDWeight=1;
00191 
00192 // Clear event variables
00193   Clear();
00194 
00195 }
00196 
00197 // Standard destructor
00198 TMarEvent::~TMarEvent()
00199 {
00200 
00201  // Delete objects
00202  if(fTrigReweight) delete fTrigReweight;
00203  if(fMCReweight) delete fMCReweight;
00204 
00205  if (electron){
00206    electron->Delete();
00207    delete electron;
00208  }
00209  if (photon){
00210    photon->Delete();
00211    delete photon;
00212  } 
00213  if (muon){
00214    muon->Delete();
00215    delete muon;
00216  }
00217  if (jet){
00218    jet->Delete();
00219    delete jet;
00220  }
00221  if (np){
00222    np->Delete();
00223    delete np;
00224  }
00225 
00226 
00227  if (GenElectron){
00228    GenElectron->Delete();
00229    delete GenElectron;
00230  }
00231  if (GenPhoton){
00232    GenPhoton->Delete();
00233    delete GenPhoton;
00234  } 
00235  if (GenMuon){
00236    GenMuon->Delete();
00237    delete GenMuon;
00238  }
00239  if (GenTau){
00240    GenTau->Delete();
00241    delete GenTau;
00242  }
00243  if (GenNp){
00244    GenNp->Delete();
00245    delete GenNp;
00246  }
00247  if (GenHfs){
00248    GenHfs->Delete();
00249    delete GenHfs;
00250  }
00251  
00252 
00253  PartCalEm.clear();
00254  PartCalMuon.clear();
00255  PartCalJet.clear();
00256  PartScaled[0].clear();
00257  PartScaled[1].clear();
00258  PartScaled[2].clear();
00259  
00260 
00261  delete lumi;
00262 
00263  
00264 };
00265 
00266 void TMarEvent::SelectParticles()
00267 {
00268 
00269 
00270 //...define electrons and photons
00271   if(find_ElePho)FindElePho();
00272 
00273 //...define muons
00274   if(find_Muon)FindMuon();
00275 
00276 //...define Jets
00277   if(find_Jet)FindJet();
00278 
00279 //...define neutral particles
00280   if(find_NP)FindNP();
00281 
00282 }
00283 
00284 void TMarEvent::CalculateDerivedVariables()
00285 {
00286  if(TMAREVENTDEBUG) printf("Enter CalculateDerivedVariables ...\n");
00287 
00288 //... Hadronic system variables
00289   const Double_t GenS = 4*GenPl*GenPp;
00290 
00291   TLorentzVector HadTotVec=HadNoJetVec+HadJetVec;
00292   Gammah = 0;
00293 
00294 //  Gammah = TMath::ACos((HadTotVec.Pt()*HadTotVec.Pt() - pow(HadTotVec.E()-HadTotVec.Pz(),2))/(HadTotVec.Pt()*HadTotVec.Pt() + pow(HadTotVec.E()-HadTotVec.Pz(),2))); 
00295   Eh =  HadTotVec.E();  
00296   Pzh = HadTotVec.Pz(); 
00297   Pth = HadTotVec.Pt(); 
00298   Phh = HadTotVec.Phi();        
00299 
00300   if (Pth>0) Thh = 2*TMath::ATan((Eh-Pzh)/Pth); 
00301   Gammah = Thh;
00302   yh  = (Eh-Pzh)/(2*GenPl);     
00303   Q2h = (Pth*Pth)/(1-yh);
00304   if(yh!=0) xh = Q2h/(yh*GenS);
00305 
00306 //... Loop over e.m. particles 
00307   IndexElScat = -1; 
00308   
00309   for(vector<TMarBody>::iterator iele=PartEm->begin();iele != PartEm->end();iele++) {
00310     // Selection for highest isolated pt electron or highest pt if only one electron
00311     IndexElScat++;
00312 
00313 //... If scattered electron found
00314     if (ModsPartEm[(*iele).GetIndexID()]->IsIsolatedLepton()){ 
00315       
00316 
00317       // Get 4-mom of electron
00318       Pte = (*iele).GetPt();    
00319       Ee =  (*iele).GetE();    
00320       Pze = (*iele).GetPz();    
00321       The = (*iele).GetTheta();    
00322       Phe = (*iele).GetPhi();    
00323      
00324       // Calc Q2e and ye
00325       ye = 1-(Ee/GenPl)*sin(The/2.)*sin(The/2.);
00326       Q2e = 4*GenPl*Ee*cos(The/2.)*cos(The/2.);
00327       if(ye!=0) xe = Q2e/(ye*GenS);
00328    
00329       // Double-angle method => Ptda
00330       Q2da = 4*pow(GenPl,2)*sin(Gammah)*(1+TMath::Cos(The))/(TMath::Sin(Gammah)+TMath::Sin(The)-TMath::Sin(Gammah+The)); 
00331       xda = (GenPl/GenPp)*(TMath::Sin(Gammah)+TMath::Sin(The)+TMath::Sin(Gammah+The))/(TMath::Sin(Gammah)+TMath::Sin(The)-TMath::Sin(Gammah+The)); 
00332       yda = Q2da / (xda*GenS);
00333       Ptda = 2*GenPl/(TMath::Tan(The/2)+((Eh-Pzh)/Pth));
00334       Eda=2*GenPl*sin(Gammah)/(sin(Gammah)+sin(The)-sin(Gammah+The));
00335  
00336       break;
00337     }
00338    
00339     
00340   } 
00341 
00342 
00343 //...The sigma method
00344   Double_t Sigma = HadTotVec.E()-HadTotVec.Pz();
00345   if(Ee!=0) ys = Sigma/(Sigma+Ee*(1-cos(The)));
00346   if(ys!=1) Q2s = pow(Ee*sin(The),2)/(1-ys);
00347   if(ys!=0) xs = Q2s/(ys*GenS);
00348 
00349 //...electron-sigma method
00350   Q2es=Q2e;
00351   xes=xs;
00352   if(xes!=0) yes=Q2es/(xes*GenS);
00353 
00354 
00355 //...calculate the total four momentum
00356  // add the hadron system
00357  TotalVec=HadTotVec;
00358 
00359  // add the isolated muons ...
00360  for (vector<TMarBody>::iterator imu=PartMuon->begin();imu != PartMuon->end();imu++){
00361      
00362       
00363      if(ModsPartMuon[(*imu).GetIndexID()]->IsIsolatedLepton()){
00364        TotalVec+=(*imu).GetFourVector();
00365      }
00366  }
00367    
00368  // add the isolated electrons ...
00369  for (vector<TMarBody>::iterator iele=PartEm->begin();iele != PartEm->end();iele++){
00370      
00371      if(ModsPartEm[(*iele).GetIndexID()]->IsIsolatedLepton()){
00372        TotalVec+=(*iele).GetFourVector(); 
00373      }
00374  }
00375 
00376  Epz=TotalVec.E()-TotalVec.Pz();
00377  Ptmiss=TotalVec.Pt();
00378 
00379 
00380 
00381 }
00382 
00383 
00384 Bool_t TMarEvent::FillFromHat()
00385 
00386 {
00388 //     Fill variables of TMarEvent from HAT information
00390 //----  Make pointers to object stored on HAT 
00391  
00392  const Double_t GenS = 4*GenPl*GenPp;
00393 
00394 
00395 //...fill variables
00396  static H1FloatPtr PtrXVert("VtxX");    // x vertex
00397  Vertex[0]=*PtrXVert;
00398 
00399  static H1FloatPtr PtrYVert("VtxY");    // y vertex
00400  Vertex[1]=*PtrYVert;
00401 
00402  static H1FloatPtr PtrZVert("VtxZ");    // Z vertex
00403  Vertex[2]=*PtrZVert;
00404 
00405  static H1FloatPtr PtrGenVtxZ("GenVtxZ");    // Z vertex
00406  GenVtxZ=*PtrGenVtxZ;
00407 
00408  static H1FloatPtr PtrEpz("Epz");        // E-Pz 
00409  Epz=*PtrEpz;   
00410 
00411  static H1FloatPtr PtrPtMiss("PtMiss"); // Ptmiss
00412  Ptmiss=*PtrPtMiss;
00413 
00414  static H1FloatPtr PtrEt("ScalEt");  //E transverse
00415  Et=*PtrEt;
00416 
00417 //---- Energy sums in calo 
00418 
00419 
00420  static H1FloatPtr PtrLArtotE("HfsClusLarE");
00421  static H1FloatPtr PtrSpatotE("HfsClusSpaE");
00422  EtotCalo= *PtrLArtotE + *PtrSpatotE;
00423  
00424  static H1FloatPtr PtrQ2e("Q2e");        // electron method 
00425  Q2e=*PtrQ2e;
00426 
00427  static H1FloatPtr Ptrye("Ye"); 
00428  ye=*Ptrye;
00429 
00430  static H1FloatPtr PtrQ2h("Q2h");        //jb method
00431  Q2h=*PtrQ2h;
00432 
00433  static H1FloatPtr Ptryh("Yh"); 
00434  yh=*Ptryh;
00435 
00436  static H1FloatPtr PtrQ2da("Q2da");      //da method
00437  Q2da=*PtrQ2da;
00438  static H1FloatPtr Ptryda("Yda"); 
00439  yda=*Ptryda;
00440 
00441  static H1FloatPtr PtrQ2s("Q2s");      //sigma method
00442  Q2s=*PtrQ2s;
00443  static H1FloatPtr Ptrys("Ys"); 
00444  ys=*Ptrys;
00445  if(ys) xs=Q2s/(ys*GenS);
00446  Pts=sqrt((1-ys)*Q2s);
00447  
00448 //...electron-sigma method
00449   Q2es=Q2e;
00450   xes=xs;
00451   if(xes!=0) yes=Q2es/(xes*GenS);
00452 
00453  
00454  // Timing
00455  static H1FloatPtr timingCJC("T0_CJC");
00456  tCJC = *timingCJC;
00457 
00458  static H1FloatPtr timingLar("T0_LAr");
00459  tLAr = *timingLar;
00460 
00461 
00462 //---- added temporary for calib ... to remove from gal declaration ?
00463  static H1FloatPtr YhGen("YhGen");
00464  yhgen=*YhGen;
00465 
00466  static H1ShortPtr genrad("GenRad");
00467  GenRad=*genrad;
00468 
00469  // Weighting factor for MC
00470  static H1FloatPtr MCPoids1("Weight1");
00471  static H1FloatPtr MCPoids2("Weight2");
00472 
00473  if(MCFlag==0) { //if Data
00474      MCWeight = *MCPoids1 ;
00475  }
00476  else { //if MC
00477 
00478   if (!*MCPoids1 || !*MCPoids2){
00479      printf("ATTENTION: Weight = 0 => Will be set to 1\n");
00480      MCWeight = LumiWeight;
00481   }
00482   else
00483 
00484      MCWeight = *MCPoids1 * (*MCPoids2)*LumiWeight;
00485      ApplyMCXSectionReweight();
00486  }
00487  
00488  
00489 // cout << "MCWeight "<< MCWeight << " "<< *MCPoids1 * (*MCPoids2)<< endl;
00490 
00491 //--- trigger
00492  static H1BytePtr PtrIl1ac("Il1ac"); 
00493  static H1BytePtr PtrIl1rw("Il1rw"); 
00494  static H1BytePtr PtrIl1te("Il1te"); 
00495  for(int i=0;i<128;i++) {
00496      Il1ac[i]=PtrIl1ac[i];
00497      Il1rw[i]=PtrIl1rw[i];
00498  }
00499  for(int i=0;i<192;i++) {
00500      Il1te[i]=PtrIl1te[i];
00501 
00502  }
00503 
00504 //... scattered electron quantities
00505  static H1FloatPtr ElecEPtr("ElecE");   
00506  Ee =  *ElecEPtr; 
00507 
00508  static H1FloatPtr ElecThetaPtr("ElecTheta");          
00509  The = *ElecThetaPtr;  
00510  
00511  static H1FloatPtr ElecPhiPtr("ElecPhi");          
00512  The = *ElecPhiPtr;  
00513 
00514   
00515  static H1FloatPtr VantiPtr("Vanti");
00516  Vanti=*VantiPtr;
00517  
00518  static H1FloatPtr VparlPtr("Vparl");
00519  Vparl=*VparlPtr;
00520 
00521 
00522 //... calorimeter quantities
00523   static H1FloatPtr PtrSpatotPt("HfsClusSpaPt");
00524   static H1FloatPtr PtrSpatotTheta("HfsClusSpaTheta");
00525   static H1FloatPtr PtrSpatotPhi("HfsClusSpaPhi");
00526   
00527   static H1FloatPtr PtrLartotPt("HfsClusLarPt");
00528   static H1FloatPtr PtrLartotTheta("HfsClusLarTheta");
00529   static H1FloatPtr PtrLartotPhi("HfsClusLarPhi");
00530   
00531 
00532   static H1FloatPtr PtrIrontotPt("HfsClusIronPt");
00533   static H1FloatPtr PtrIrontotTheta("HfsClusIronTheta");
00534   static H1FloatPtr PtrIrontotPhi("HfsClusIronPhi");
00535   static H1FloatPtr PtrIrontotE("HfsClusIronE");
00536 
00537   TLorentzVector Lar(0.,0.,0.,0.);
00538   TLorentzVector Spacal(0.,0.,0.,0.);
00539   TLorentzVector Iron(0.,0.,0.,0.);
00540   Spacal.SetPtEtaPhiE(*PtrSpatotPt,-log(tan(*PtrSpatotTheta/2)),*PtrSpatotPhi,*PtrSpatotE);
00541   Lar.SetPtEtaPhiE(*PtrLartotPt,-log(tan(*PtrLartotTheta/2)),*PtrLartotPhi,*PtrLArtotE);
00542   Iron.SetPtEtaPhiE(*PtrIrontotPt,-log(tan(*PtrIrontotTheta/2)),*PtrIrontotPhi,*PtrIrontotE);
00543 
00544   TLorentzVector All=Spacal+Lar+Iron;
00545 
00546   ptIronfrac=Iron.Pt()/All.Pt();
00547   EIronfrac= Iron.E()/All.E();
00548   ptSpacalfrac=Spacal.Pt()/All.Pt();
00549   ESpacalfrac=Spacal.E()/All.E();
00550  
00551  return(kTRUE);
00552 }
00553 
00554 
00555 Bool_t TMarEvent::MCGenSelection(TMarCutMCList& List)
00556 {
00557     
00558 //... loop on one cut list (AND on the cuts). Return false
00559 //... if at least one of the cut is not fulfilled. Change the
00560 //... event weight otherwise   
00561 
00562      TLorentzVector BeamElectron(0.0,0.0,-1.*GenPl,GenPl);
00563      TLorentzVector BeamProton(0.0,0.0,GenPp,sqrt(GenPp*GenPp+0.880354323));
00564      Double_t lumi_weight=0;
00565      Double_t event_value=-1;   
00566          
00567      for(list<TMarCutMC>::iterator it= List.CutList.begin();it != List.CutList.end(); it++)
00568      {
00569    
00570         TMarCutMC& cut=(*it); 
00571        
00572 //...cut on electron-photon system (compton)
00573         if(cut.GetVariable()=="MELEPHO"){ 
00574           if(GenPhoton->GetEntries()==1&&GenElectron->GetEntries()==1){
00575             Float_t GenMass2ElePho=0;
00576             TMarBody *Ele=(TMarBody *)GenElectron->At(0);
00577             TMarBody *Pho=(TMarBody *)GenPhoton->At(0);
00578             GenMass2ElePho = pow(Pho->GetE()+Ele->GetE(),2)-pow(Pho->GetPx()+Ele->GetPx(),2)
00579                            -pow(Pho->GetPy()+Ele->GetPy(),2)-pow(Pho->GetPz()+Ele->GetPz(),2);
00580             if (GenMass2ElePho>0) {
00581                event_value= sqrt(GenMass2ElePho);
00582                lumi_weight=cut.Evaluate(event_value);
00583                
00584             } else throw string("TMarEvent::FillFromMods = Attention generated mass of system  photon-electron < 0");
00585           }  else throw string("TMarEvent::FillFromMods = Attention no system photon-electron in generated particles"); 
00586         }
00587 
00588 //... cut on mu-mu system
00589        else if(cut.GetVariable()=="MMUMU"){
00590           
00591           if(GenMuon->GetEntries()>1){
00592             Float_t GenMass2MuMu=0;
00593             TMarBody *Mu1=(TMarBody *)GenMuon->At(0);
00594             TMarBody *Mu2=(TMarBody *)GenMuon->At(1);            
00595             GenMass2MuMu = pow(Mu1->GetE()+Mu2->GetE(),2)-pow(Mu1->GetPx()+Mu2->GetPx(),2)
00596                      -pow(Mu1->GetPy()+Mu2->GetPy(),2)-pow(Mu1->GetPz()+Mu2->GetPz(),2);
00597             if(GenMass2MuMu>0){ 
00598                event_value = sqrt(GenMass2MuMu); 
00599                lumi_weight=cut.Evaluate(event_value);
00600             } else throw string("TMarEvent::FillFromMods = Attention generated mass of system mu-mu < 0");
00601             
00602           } else {
00603                 cout << "TMarEvent::FillFromMods = Attention no system mu-mu in generated particles"<< endl;
00604                 cout << "event num "<< num <<endl;
00605                 for (Int_t iPartMC=0;iPartMC<ModsPartMC.GetEntries();iPartMC++){
00606                 
00607                 ModsPartMC[iPartMC]->Print();
00608                 
00609         //           cout << "particle " << ModsPartMC[iPartMC]->Print()<< endl;
00610                 }
00611          }   
00612        } 
00613        
00614 //... cut on Q2
00615        else if(cut.GetVariable()=="Q2"){
00616           static H1FloatPtr PtrQ2Gen("Q2Gki"); // Generated Q2 from Gki bank
00617           event_value = (*PtrQ2Gen);
00618           lumi_weight=cut.Evaluate(event_value);
00619 //        printf("cut Q2=%f w=%f \n",*PtrQ2Gen,lumi_weight);
00620        }
00621 //... cut on new Q2 (for RAPGAP files)
00622        else if(cut.GetVariable()=="RAPQ2"){ 
00623           event_value = -(BeamElectron-GenScatEle).M2();
00624           lumi_weight=cut.Evaluate(event_value);
00625           
00626        }
00627 //... cut on Y
00628        else if(cut.GetVariable()=="Y"){
00629           static H1FloatPtr PtryGen("YGki"); 
00630           event_value = (*PtryGen);
00631           lumi_weight=cut.Evaluate(event_value);
00632        }
00633 //... cut on new Y
00634        else if(cut.GetVariable()=="RAPY"){
00635           // NewY
00636           event_value = (BeamProton*(BeamElectron-GenScatEle))/(BeamProton*BeamElectron);
00637           lumi_weight=cut.Evaluate(event_value);
00638        }
00639 //... cut on Phat for g-P files
00640        else if(cut.GetVariable()=="PHAT"){
00641           static H1FloatPtr PtrPt2HatGki("Pt2HatGki"); 
00642           event_value = sqrt(*PtrPt2HatGki);
00643           lumi_weight=cut.Evaluate(event_value);  
00644 //        printf("cut PHAT=%f w=%f \n",event_value,lumi_weight);
00645        }
00646        else {
00647           throw string("TMarEvent::FillFromMods = Error in the generator variable cut name! - stop");
00648        }
00649 
00650 //... verify the cut
00651 
00652        if(lumi_weight==0) return kFALSE;
00653        else LumiWeight=lumi_weight*LumiWeight;
00654   
00655   }
00656 
00657    return kTRUE;
00658 }
00659 
00660 Bool_t TMarEvent::FillGeneratedParticlesAndCut()
00661 {
00663 //     Fill generated quantities and evaluate the MC cut
00664 //     return false if one of the cut is not respected
00666   if(TMAREVENTDEBUG) printf("Enter FillGeneratedParticlesAndCut ...\n");
00667 
00668   Bool_t MCCutOk=kTRUE;
00669 
00670 
00671   if(MCFlag>0) {
00672 
00673     FindGeneratedParticles();   //--- retrieve generated parts from mods
00674 
00675 //... if cuts on generated variables are defined:
00676     if(genCut.size()>0){
00677 
00678 
00679 //... cut on generated variable for MC. Loop on all the possible cut  lists
00680       LumiWeight=1;
00681       for(list<TMarCutMCList>::iterator il= genCut.begin();il != genCut.end(); il++)
00682       {
00683           MCCutOk=MCGenSelection(*il);
00684 
00685 //... stop as soon as a cut configuration is fulfilled
00686           if(MCCutOk==kTRUE) break;
00687 
00688       }
00689 
00690     }
00691  }
00692 
00693  return MCCutOk;
00694 }
00695 
00696 
00697 Int_t TMarEvent::FindElePho()
00698 {
00699 // Find ele/pho in cand list and return:
00700 // Clean version (eliminate noise runs)
00701 // * 0 if ok
00702 
00703   Int_t nEm=PartEm->size();
00704   
00705   if(nEm<=0)return(kFALSE);
00706 
00707   // ID of found particle
00708   Int_t ID = -1,IDNow=-1;
00709 
00710   // Vars of candidates
00711   Float_t PtPartNow,EtaPartNow,ThetaPartNow,PhiPartNow,EPartNow;
00712 
00713   // For e/LW track matching
00714   Float_t EtaTr,PhiTr;
00715   Float_t EtaMatchTr=0,PhiMatchTr=0;
00716   Int_t NMatchLWTr,NMatchAllTr,NMatchDTNVTr;
00717  
00718 
00719 
00720 
00721   // Distance 
00722   Float_t R2Particles,RParticles,DiffEta,DiffPhi;
00723 
00724   
00725   for (vector<TMarBody>::iterator iele=PartEm->begin();iele!=PartEm->end();iele++) {
00726 
00727     // Init
00728     
00729     NMatchLWTr = 0;
00730     NMatchAllTr = 0;
00731     NMatchDTNVTr = 0;
00732     ID = -1;  
00733     
00734   
00735     
00736     
00737     IDNow=(*iele).GetIndexID();
00738     EPartNow = (*iele).GetE();
00739     PtPartNow  = (*iele).GetPt();
00740     ThetaPartNow = (*iele).GetTheta();
00741     EtaPartNow  = -log(TMath::Tan(ThetaPartNow/2));
00742     PhiPartNow  = (*iele).GetPhi();
00743           
00744     if (PtPartNow > PtMinEle && ModsPartEm[IDNow]->IsIsolatedLepton() == kTRUE){
00745 
00746 //... If e-cluster in central region ask for LW or DTRA track linked to it and
00747 //... no other LW track in same region
00748       if (ThetaPartNow*MYR2D > ThetaMinEle && ThetaPartNow*MYR2D < ThetaMaxEle){ 
00749         
00750          if (ModsPartEm[IDNow]->GetTrType() == 3)
00751             NMatchDTNVTr = ModsPartEm[IDNow]->GetNbMatchingTr(); 
00752          
00753 //... Check if DTRA or LW track pointing to cluster
00754          if (ModsPartEm[IDNow]->GetTrType() == 1 || ModsPartEm[IDNow]->GetTrType() == 2) {
00755             NMatchAllTr = ModsPartEm[IDNow]->GetNbMatchingTr();
00756             PhiMatchTr = ModsPartEm[IDNow]->GetTrPhi();
00757             EtaMatchTr = ModsPartEm[IDNow]->GetEta();
00758           }         
00759 //...look for LW track in a cone around the electron
00760 //... Loop over all LW tracks
00761         
00762           for(Int_t iTr=0;iTr<ModsPartSelTrack.GetEntries();iTr++){
00763 
00764 //... Get Eta/phi of track iTr
00765              EtaTr = ModsPartSelTrack[iTr]->GetEta();
00766              PhiTr =   ModsPartSelTrack[iTr]->GetPhi();
00767             
00768 //... Calculate distance to: - e-cluster for photon
00769 //...                        - e-track for ele
00770 
00771 //... Distance to e-track
00772              if (NMatchAllTr > 0){
00773                DiffEta = EtaMatchTr-EtaTr;
00774                DiffPhi = GetRealPhi(PhiMatchTr-PhiTr);
00775              }
00776 //...Distance to e-cluster
00777              else{
00778                DiffEta = ModsPartEm[IDNow]->GetEta()-EtaTr;
00779                DiffPhi = GetRealPhi(ModsPartEm[IDNow]->GetTrPhi()-PhiTr);
00780              }
00781 //... Final distance in eta-phi plane
00782              R2Particles = DiffEta*DiffEta + DiffPhi*DiffPhi;
00783              if (R2Particles > 0.0001)
00784                RParticles = TMath::Sqrt(R2Particles);           
00785              else
00786                RParticles = 0;
00787            
00788            
00789 // Count tracks that match in this region R < 0.5
00790 // => Isolation
00791              
00792              if (RParticles < ConeEle) NMatchLWTr++;
00793          
00794          }
00795 //... central cluster must have an isolated  LW  track
00796          h1cerr(2,"   - NMatchAllTr " << NMatchAllTr <<
00797                      ", NMatchLWTr "  << NMatchLWTr <<
00798                      ", NMatchDTNVTr "  << NMatchDTNVTr); 
00799             if (NMatchAllTr == 1 && NMatchLWTr == 1)
00800                 ID = IDELE;    
00801 //... central gamma must have no track
00802             else if (NMatchDTNVTr == 0 && NMatchLWTr == 0 && NMatchAllTr == 0)
00803                 ID = IDPHO;
00804        }
00805     
00806   
00807        if (ID == IDPHO ){
00808       
00809         AddParticleToList(photon,IDNow,ID,EPartNow,PtPartNow,EtaPartNow,PhiPartNow);
00810       }
00811       if (ID == IDELE ){
00812       
00813         AddParticleToList(electron,IDNow,ID,EPartNow,PtPartNow,EtaPartNow,PhiPartNow);
00814       }
00815 
00816     } // end if (PtPartNow
00817   } // end of loop PartEm
00818   
00819 
00820   return(kTRUE);
00821 
00822 
00823 }
00824  
00825 void TMarEvent::SortPt(TClonesArray *list)
00826 {
00827   Int_t NbBodies=list->GetEntries();
00828   
00829   if(NbBodies<=1)return;
00830 
00831   TClonesArray* oldlist=new TClonesArray("TMarBody",1);
00832   
00833   Int_t* tkIndex= new Int_t [NbBodies];
00834   Double_t* tkPt= new Double_t [NbBodies];
00835 
00836   for(Int_t itk=0;itk<NbBodies;itk++) {
00837              tkPt[itk]=((TMarBody*)list->At(itk))->GetPt();
00838              tkIndex[itk]=itk;
00839              TMarBody *body = (TMarBody *)list->AddrAt(itk);
00840              int n=oldlist->GetLast()+1;
00841              new (oldlist->AddrAt(n)) TMarBody(body);   
00842   }
00843 //   list->Clear();
00844   list->Delete();
00845 
00846 
00847   TMath::Sort(NbBodies,tkPt,tkIndex,kTRUE);
00848 //   printf("\n Nbodies=%d :",NbBodies);
00849   for(Int_t itk=0;itk<oldlist->GetEntries();itk++) {
00850      Int_t id=tkIndex[itk];
00851      TMarBody *Body = (TMarBody*) oldlist->UncheckedAt(id);
00852      int n=list->GetLast()+1;
00853 //      printf("itk=%d tkIndex=%d Pt=%f %f n=%d/ ",itk,tkIndex[itk],((TMarBody*)oldlist->At(id))->GetPt(),tkPt[tkIndex[itk]],n);
00854      new (list->AddrAt(n)) TMarBody(Body);   
00855             
00856   }
00857 
00858   oldlist->Delete();
00859   delete oldlist;
00860   delete tkIndex;
00861   delete tkPt;
00862 
00863 }
00864   
00865 Int_t TMarEvent::FindMuon()
00866 {
00867 // Find muon
00868   if(PartMuon->size()<=0)return(kFALSE);
00869   // Vars of candidates
00870   Float_t PtPartNow,EtaPartNow,ThetaPartNow,PhiPartNow,EPartNow;
00871 
00872 
00873   // ID of particle
00874   Int_t ID = -1;
00875 
00876   for (vector<TMarBody>::iterator iMu =PartMuon->begin() ;iMu!=PartMuon->end();iMu++) {
00877 
00878     // Init 
00879     ID = -1;
00880 
00881 
00882     EPartNow = (*iMu).GetE();
00883     PtPartNow  = (*iMu).GetPt();
00884     ThetaPartNow = (*iMu).GetTheta();
00885     EtaPartNow  = -log(TMath::Tan(ThetaPartNow/2));
00886     PhiPartNow  = (*iMu).GetPhi();
00887 
00888 
00889     if (PtPartNow > PtMinMu ){
00890 
00891 // If muonic particle check:
00892 // - Theta-angle 
00893 // - Tracks (vertex fitted) around muon
00894     
00895         if (ThetaPartNow*MYR2D > ThetaMinMu && ThetaPartNow*MYR2D < ThetaMaxMu){
00896           // Get Muon info
00897           if (ModsPartMuon[(*iMu).GetIndexID()]->GetDTrTr() > DTrTrMinMu)
00898             ID = IDMU;
00899         }
00900 
00901       // Check ID 
00902       if (ID == IDMU){
00903         AddParticleToList(muon,(*iMu).GetIndexID(),ID,EPartNow,PtPartNow,EtaPartNow,PhiPartNow);
00904       }
00905 
00906     } // End of if (PtPartNow
00907   }
00908 
00909 
00910   return(0);
00911 
00912 }
00913  
00914   
00915 Int_t TMarEvent::FindNP()
00916 {
00917 // Find neutrino
00918 
00919   // Vars of candidates
00920   Float_t PtPartNow,EtaPartNow,ThetaPartNow,PhiPartNow,EPartNow;
00921 
00922   // Count NP's
00923   if (Ptmiss >PtMinNp && (55.2-Epz)>0.1){
00924 
00925     // Init
00926     ThetaPartNow = 2*TMath::ATan((55.2-Epz)/Ptmiss);
00927     if (ThetaPartNow*MYR2D >  ThetaMinNp&& ThetaPartNow*MYR2D < ThetaMaxNp){
00928  
00929        PtPartNow = Ptmiss;
00930        EPartNow = sqrt(PtPartNow*PtPartNow+(55.2-Epz)*(55.2-Epz));
00931        EtaPartNow = -log(TMath::Tan(ThetaPartNow/2));
00932        PhiPartNow = TMath::ATan2(-TotalVec.Py(),-TotalVec.Px());;
00933 
00934     // Add to particle list
00935 
00936        AddParticleToList(np,INDEXNP,IDNP,EPartNow,PtPartNow,EtaPartNow,PhiPartNow);
00937     }
00938   }
00939 
00940   return(kTRUE);
00941 
00942 }
00943 
00944   
00945 Int_t TMarEvent::FindJet()
00946 {
00947 // Find jet
00948 // Clean version (eliminate noise runs)
00949   if(PartJet->size()<=0)return(kFALSE);
00950   // Vars of candidates
00951   Float_t PtPartNow,EtaPartNow,ThetaPartNow,PhiPartNow,EPartNow;
00952 
00953 
00954   // For jet anal: jet size and duaghter particle vars
00955   Float_t RJet;
00956   Int_t NCountLWTrJets;
00957   Float_t ECandDaughter,EtaCandDaughter,PhiCandDaughter,TheCandDaughter;
00958   Float_t Phi2DiffDaughter,Eta2DiffDaughter;
00959 
00960 
00961   // Loop over Jets
00962   for(vector<TMarBody>::iterator iJets=PartJet->begin();iJets!=PartJet->end();iJets++){
00963 
00964     // Get pointer to jet 
00965     
00966     TMarBody* Part=&(*iJets);
00967     
00968     // Get kinematics of jet 
00969     PtPartNow = Part->GetPt();
00970     EPartNow = Part->GetE();
00971     ThetaPartNow = Part->GetTheta();
00972     EtaPartNow = - log(TMath::Tan(ThetaPartNow/2));
00973     PhiPartNow = Part->GetPhi();
00974 
00975 
00976 
00977     // Calculus of jet size (calib only for EJet)
00978     RJet = 0.0;
00979     NCountLWTrJets = 0;
00980    
00981     const TObjArray *ListDaughters = ModsPartJet[Part->GetIndexID()]->GetParticles();
00982     Int_t NumDaughters = ListDaughters->GetEntries();
00983     for (Int_t iDaughter=0;iDaughter<NumDaughters;iDaughter++){
00984 
00985       // Get det object of particle candidate
00986 
00987       const H1PartCand *candDaughter = (H1PartCand *) ListDaughters->At(iDaughter);
00988 
00989       // For each daughter=hadronic object get eta,phi
00990       ECandDaughter = candDaughter->GetE();
00991       PhiCandDaughter = candDaughter->GetPhi();
00992       TheCandDaughter = candDaughter->GetTheta();
00993       EtaCandDaughter = -log(TMath::Tan(TheCandDaughter/2));
00994       Phi2DiffDaughter = GetRealPhi(PhiCandDaughter-PhiPartNow)*GetRealPhi(PhiCandDaughter-PhiPartNow);
00995       Eta2DiffDaughter = (EtaCandDaughter-EtaPartNow)*(EtaCandDaughter-EtaPartNow);
00996       RJet += ECandDaughter*(TMath::Sqrt(Phi2DiffDaughter+Eta2DiffDaughter));
00997 
00998       // Count tracks in jet
00999       if (candDaughter->IsSelTrack() == kTRUE)
01000         NCountLWTrJets++;
01001 
01002  
01003     }
01004     
01005 
01006     // Take R jet from detector variable
01007     RJet = ModsPartJet[Part->GetIndexID()]->GetJetSize();
01008 
01009 
01010     // Preselect jet if above min pt, in LAr, and size over 0.1 in cracks
01011     // In add EM fraction must be less 0.9
01012     if (PtPartNow > PtMinJet && (ThetaPartNow*MYR2D > ThetaMinJet && ThetaPartNow*MYR2D < ThetaMaxJet)){
01013       if (RJet > RMinJet || ModsPartJet[Part->GetIndexID()]->GetEmFraction() < EmFracJet){
01014 
01015             AddParticleToList(jet,Part->GetIndexID(),IDJET,EPartNow,PtPartNow,EtaPartNow,PhiPartNow); 
01016 
01017       }
01018     } // End of if (PtPartCalib
01019   }   
01020 
01021 
01022 
01023   return(kTRUE);
01024 
01025 }
01026   
01027   
01028 Float_t TMarEvent::MinDisttoJets(Float_t Theta,Float_t Phi)
01029 {
01031 //     calculate min distance to jets in eta-phi plan
01033 
01034   Float_t ThetaJ,PhiJ,EtaJ,DiffEta,DiffPhi,R2Particles,RParticles;
01035   Float_t MinR=99999.;
01036   
01037   for(Int_t iJ=0;iJ<ModsPartJet->GetEntries();iJ++){
01038 //... Get Eta/phi of track iTr
01039         ThetaJ = ModsPartJet[iJ]->GetTheta();
01040         PhiJ =   ModsPartJet[iJ]->GetPhi();
01041         EtaJ = -log(TMath::Tan(ThetaJ/2));
01042             
01043 //... Calculate distance to: - e-cluster for photon
01044 //...                        - e-track for ele
01045 //...Distance to e-cluster
01046         DiffEta = -log(TMath::Tan(Theta/2))-EtaJ;
01047         DiffPhi = TMath::ACos(TMath::Cos(Phi-PhiJ));
01048 
01049 //... Final distance in eta-phi plane
01050         R2Particles = DiffEta*DiffEta + DiffPhi*DiffPhi;
01051         if (R2Particles > 0.0001) RParticles = TMath::Sqrt(R2Particles);           
01052         else RParticles = 0;
01053 // Count tracks that match in this region R < 0.5
01054 
01055 // => Isolation
01056         if (RParticles < MinR && RParticles!=0) MinR=RParticles;
01057   }
01058   
01059   return MinR;
01060 }
01061 
01062 
01063 
01064 Float_t TMarEvent::MinDisttoTracks(Float_t Theta,Float_t Phi)
01065 {
01067 //     calculate min distance to LW tracks in eta-phi plan
01069 
01070   Float_t ThetaJ,PhiJ,EtaJ,DiffEta,DiffPhi,R2Particles,RParticles;
01071   Float_t MinR=99999.;
01072   
01073   for(Int_t iJ=0;iJ<ModsPartSelTrack->GetEntries();iJ++){
01074 //... Get Eta/phi of track iTr
01075         ThetaJ = ModsPartSelTrack[iJ]->GetTheta();
01076         PhiJ =   ModsPartSelTrack[iJ]->GetPhi();
01077         EtaJ = -log(TMath::Tan(ThetaJ/2));
01078             
01079 //... Calculate distance to: - e-cluster for photon
01080 //...                        - e-track for ele
01081 //...Distance to e-cluster
01082         DiffEta = -log(TMath::Tan(Theta/2))-EtaJ;
01083         DiffPhi = TMath::ACos(TMath::Cos(Phi-PhiJ));
01084 
01085 //... Final distance in eta-phi plane
01086         R2Particles = DiffEta*DiffEta + DiffPhi*DiffPhi;
01087         if (R2Particles > 0.0001) RParticles = TMath::Sqrt(R2Particles);           
01088         else RParticles = 0;
01089 // Count tracks that match in this region R < 0.5
01090 // => Isolation
01091         if (RParticles < MinR && RParticles!=0) MinR=RParticles;
01092   }
01093 
01094   
01095   return MinR;
01096 }
01097 
01098 
01099 
01100 Float_t TMarEvent::CalcDistDTNV(const H1PartEm* EMParticle)
01101 {
01102 // Calc. distance to nearest DTNV track from particle Particle (cm)
01103 
01104   // Init
01105   Float_t ThetaDTNV = 0;
01106   Float_t PhiDTNV = 0;
01107   Float_t PhiExtraDTNV = 0;
01108   Float_t XExtraOut = 0;
01109   Float_t YExtraOut = 0;
01110   Float_t ZExtraOut = 0;    
01111   Float_t DiffDCA = 0.0;
01112   Float_t PI = TMath::Pi();
01113   Float_t Alpha=0,Dlon=0;
01114 
01115    // Get vertex of event
01116   TVector3 XStart(0.,0.,0.);
01117   if(Vertex[2] != -999999) {
01118      XStart.SetX(Vertex[0]);
01119      XStart.SetY(Vertex[1]);
01120      XStart.SetZ(Vertex[2]);
01121   }
01122 
01123 
01124 
01125   // Distance vars
01126   Float_t DiffDCAMin = 999999.9;
01127 
01128   // Get e.m.-part position
01129   TVector3 ClusPos(EMParticle->GetXClus(),EMParticle->GetYClus(),EMParticle->GetZClus());
01130 
01131 //  printf("Clus(x,y,z) = %f %f %f \n",ClusPos.X(),ClusPos.Y(),ClusPos.Z());
01132 
01133   // Loop over all DTNV tracks and calc distance
01134   Int_t NumDTNV = ModsDTNVTracks->GetEntries();
01135   for (Int_t iTr = 0;iTr < NumDTNV;iTr++){
01136 
01137     // Get DTNV track
01138     H1PartSelTrack *DTNVTr = (H1PartSelTrack *) ModsDTNVTracks->At(iTr);
01139 
01140     // Get theta/phi of DTNV track
01141     if (isnan(DTNVTr->GetTheta()))
01142       continue;
01143     ThetaDTNV = DTNVTr->GetTheta();
01144     PhiDTNV = DTNVTr->GetPhi();
01145 
01146     // Extrapolate track into calo
01147     CalcExtraPolTrack(XStart,DTNVTr,XExtraOut,YExtraOut,ZExtraOut,ThetaDTNV,PhiExtraDTNV);
01148 
01149     // Get end position
01150     TVector3 EndPos(XExtraOut,YExtraOut,ZExtraOut);
01151 
01152     // Calc DCA
01153     //   printf("End(x,y,z) = %f %f %f\n",EndPos.X(),EndPos.Y(),EndPos.Z());
01154     CalcDCA(ClusPos,EndPos,ThetaDTNV,PhiExtraDTNV,DiffDCA,Alpha,Dlon);
01155 //    printf("Found DTNV(%d): Pt = %f Th = %f Ph = %f PhExtra = %f RJet = %f\n",iTr,DTNVTr->GetPt(),DTNVTr->GetTheta()*MYR2D,DTNVTr->GetPhi()*MYR2D,PhiExtraDTNV*MYR2D,DiffDCA);
01156 
01157     // Compare DCA to already found track 
01158     // => track must be in direction of cluster, not opposit
01159     if(DiffDCA < DiffDCAMin && Alpha<3*PI/8){
01160       DiffDCAMin = DiffDCA;
01161     }
01162   }
01163 
01164   // Return distance to nearest DTNV tracks
01165   return(DiffDCAMin);
01166 }
01167 
01168 
01169 Float_t TMarEvent::CalcDistTracks(H1PartSelTrackArrayPtr Tracks,const H1PartEm* EMParticle)
01170 {
01171 // Calc. distance to nearest DTNV track from particle Particle (cm)
01172 
01173   // Init
01174   Float_t ThetaDTNV = 0;
01175   Float_t PhiDTNV = 0;
01176   Float_t PhiExtraDTNV = 0;
01177   Float_t XExtraOut = 0;
01178   Float_t YExtraOut = 0;
01179   Float_t ZExtraOut = 0;    
01180   Float_t DiffDCA = 0.0;
01181   Float_t PI = TMath::Pi();
01182   Float_t Alpha=0,Dlon=0;
01183 
01184    // Get vertex of event
01185    TVector3 XStart(0.,0.,0.);
01186    if(Vertex[2] != -999999) {
01187      XStart.SetX(Vertex[0]);
01188      XStart.SetY(Vertex[1]);
01189      XStart.SetZ(Vertex[2]);
01190    }
01191 
01192 
01193 
01194   // Distance vars
01195   Float_t DiffDCAMin = 999999.9;
01196 
01197   // Get e.m.-part position
01198   TVector3 ClusPos(EMParticle->GetXClus(),EMParticle->GetYClus(),EMParticle->GetZClus());
01199 
01200 //  printf("Clus(x,y,z) = %f %f %f \n",ClusPos.X(),ClusPos.Y(),ClusPos.Z());
01201 
01202   // Loop over all DTNV tracks and calc distance
01203   Int_t NumDTNV = Tracks->GetEntries();
01204   for (Int_t iTr = 0;iTr < NumDTNV;iTr++){
01205 
01206     // Get DTNV track
01207     H1PartSelTrack *DTNVTr = (H1PartSelTrack *) Tracks->At(iTr);
01208 
01209     // Get theta/phi of DTNV track
01210     if (isnan(DTNVTr->GetTheta()))
01211       continue;
01212     ThetaDTNV = DTNVTr->GetTheta();
01213     PhiDTNV = DTNVTr->GetPhi();
01214 
01215     // Extrapolate track into calo
01216     CalcExtraPolTrack(XStart,DTNVTr,XExtraOut,YExtraOut,ZExtraOut,ThetaDTNV,PhiExtraDTNV);
01217 
01218     // Get end position
01219     TVector3 EndPos(XExtraOut,YExtraOut,ZExtraOut);
01220 
01221     // Calc DCA
01222     //   printf("End(x,y,z) = %f %f %f\n",EndPos.X(),EndPos.Y(),EndPos.Z());
01223     CalcDCA(ClusPos,EndPos,ThetaDTNV,PhiExtraDTNV,DiffDCA,Alpha,Dlon);
01224 //    printf("Found DTNV(%d): Pt = %f Th = %f Ph = %f PhExtra = %f RJet = %f\n",iTr,DTNVTr->GetPt(),DTNVTr->GetTheta()*MYR2D,DTNVTr->GetPhi()*MYR2D,PhiExtraDTNV*MYR2D,DiffDCA);
01225 
01226     // Compare DCA to already found track 
01227     // => track must be in direction of cluster, not opposit
01228     if(DiffDCA < DiffDCAMin && Alpha<3*PI/8){
01229       DiffDCAMin = DiffDCA;
01230     }
01231   }
01232 
01233   // Return distance to nearest DTNV tracks
01234   return(DiffDCAMin);
01235 }
01236 
01237 
01238 void TMarEvent::CalcDCA(TVector3 ClusStart, TVector3 TrackEnd,Float_t Theta, Float_t Phi,Float_t &Dca,Float_t &Alpha,Float_t &Dlon) const
01239 {
01240 // Find DCA between track (given by TrackEnd) and cluster (given by ClusStart)
01241 //  Double_t Dlon,Alpha;
01242 
01243   // Vector speci. direction along track
01244   Float_t Rhat[3];
01245 
01246   // Vector between ClusStart and TrackEnd
01247   Float_t StaEnd[3];
01248 
01249   // Norm and Norm squared of StaEnd
01250   Float_t NormStaEnd, NormSqStaEnd;
01251 
01252   // Looping variable
01253   Int_t Index;
01254 
01255   // Init
01256   for(Index=0;Index<3;Index++){
01257       Rhat[Index] = 0;
01258       StaEnd[Index] = 0;
01259   }
01260   NormStaEnd = 0;
01261   NormSqStaEnd = 0;
01262 
01263   // Set up vector speci. direction along track
01264   Rhat[0] = TMath::Sin(Theta)*TMath::Cos(Phi);
01265   Rhat[1] = TMath::Sin(Theta)*TMath::Sin(Phi);
01266   Rhat[2] = TMath::Cos(Theta);
01267   // Set up vector between ClusStart and TrackEnd
01268   StaEnd[0] = ClusStart.X() - TrackEnd.X();
01269   StaEnd[1] = ClusStart.Y() - TrackEnd.Y();
01270   StaEnd[2] = ClusStart.Z() - TrackEnd.Z();
01271   NormSqStaEnd = (StaEnd[0]*StaEnd[0])+(StaEnd[1]*StaEnd[1])+(StaEnd[2]*StaEnd[2]);
01272 
01273   // Cluster and track are well separated
01274   if (NormSqStaEnd>0.0001){
01275      NormStaEnd = TMath::Sqrt(NormSqStaEnd);
01276      // Find sclar product between Rhat and StaEnd
01277      Dlon = (StaEnd[0]*Rhat[0])+(StaEnd[1]*Rhat[1])+(StaEnd[2]*Rhat[2]);
01278      Alpha = TMath::ACos(Dlon/NormStaEnd);
01279      Dca = TMath::Sin(Alpha)*NormStaEnd;
01280     } else{  // Cluster is at exact position of track entry point
01281       Alpha = 0;
01282       Dca = 0;
01283       Dlon = 0;
01284     }
01285 }
01286 
01287 
01288 
01289  
01290 void TMarEvent::FindGeneratedParticles()
01291 {
01293 //     To rtrieve generated particles from MODS  
01295   if(TMAREVENTDEBUG) printf("Enter FindGeneratedParticles ...\n");
01296 
01297   // Generated quantities
01298   Float_t GenE=0,GenPt=0,GenEta=0,GenPhi=0,GenTheta=0;
01299 
01300   // Loop over generated particles
01301 
01302   for (Int_t iPartMC=0;iPartMC<ModsPartMC.GetEntries();iPartMC++){     
01303 
01304      if(TMath::Abs(ModsPartMC[iPartMC]->GetPDG())>=81
01305           && TMath::Abs(ModsPartMC[iPartMC]->GetPDG())<=100) continue;//--- to remove MC's special parts.
01306      
01307      if(ModsPartMC[iPartMC]->GetStatus() == 0){ //--- fill stable particles
01308  //     cout << "stable particle " << ModsPartMC[iPartMC]->GetPDG()<< endl;
01309       
01310        GenE = ModsPartMC[iPartMC]->GetE();
01311        GenTheta = ModsPartMC[iPartMC]->GetTheta();
01312        GenPt = ModsPartMC[iPartMC]->GetPt();
01313        if(GenE==0. || GenTheta==0. || GenPt==0.) continue; //--- remove part with Pt=0
01314 //        GenEta = ModsPartMC[iPartMC]->GetEta();
01315        GenEta = -log(TMath::Tan((GenTheta)/2));
01316        GenPhi = ModsPartMC[iPartMC]->GetPhi();
01317 //         printf("%d %f %f \n",ModsPartMC[iPartMC]->GetPDG(),GenPt,GenEta);
01318 
01319       // get stable elecron quantities
01320       if (TMath::Abs(ModsPartMC[iPartMC]->GetPDG()) == 11 ){
01321  
01322         AddParticleToList(GenElectron,iPartMC,IDGELE,GenE,GenPt,GenEta,GenPhi); 
01323 
01324         if(ModsPartMC[iPartMC]->GetMother1()==0){ // get scattered electron (Mo1=Pmo%10000-1)
01325             GenScatEle.SetPx(ModsPartMC[iPartMC]->GetPx());  
01326             GenScatEle.SetPy(ModsPartMC[iPartMC]->GetPy());  
01327             GenScatEle.SetPz(ModsPartMC[iPartMC]->GetPz());  
01328             GenScatEle.SetE(ModsPartMC[iPartMC]->GetE());  
01329 
01330         }
01331       }
01332 
01333       // Get stable photon quantitites
01334       else if (ModsPartMC[iPartMC]->GetPDG() == 22 ){
01335   
01336      //--- Test isolation to other generated hadrons ....
01337 
01338        Float_t NumGenHFSCone=0;
01339        for (Int_t ih=0;ih<ModsPartMC.GetEntries();ih++){
01340 
01341               if(ih==iPartMC) continue;
01342               if (ModsPartMC[ih]->GetStatus()==0&&ModsPartMC[ih]->GetE()!=0.&&ModsPartMC[ih]->GetTheta()!=0.
01343                                                 &&ModsPartMC[ih]->GetPt()!=0.){
01344                 if (abs(ModsPartMC[ih]->GetPDG())<11 || abs(ModsPartMC[ih]->GetPDG())>16){ //---photons are also in ..
01345                    Float_t ThisEta = -log(TMath::Tan((ModsPartMC[ih]->GetTheta())/2));
01346                    Float_t ThisPhi = ModsPartMC[ih]->GetPhi();
01347                    Float_t DiffEta = GenEta-ThisEta;
01348                    Float_t DiffPhi = TMath::ACos(TMath::Cos(GenPhi-ThisPhi));
01349                    Float_t RParticles;
01350                    Float_t R2Particles = DiffEta*DiffEta + DiffPhi*DiffPhi;
01351                    if (R2Particles > 0.0001) RParticles = TMath::Sqrt(R2Particles);
01352                    else RParticles = 0;
01353                    if(RParticles<0.5) {//---- cone of 0.1 around generated
01354                      NumGenHFSCone+=1;
01355                    }
01356                 }
01357               }
01358         }
01359 
01360 //---- Add GenPhoton only if isolated to other gen hadrons by 0.5     
01361         if(NumGenHFSCone==0) { 
01362            AddParticleToList(GenPhoton,iPartMC,IDGELE,GenE,GenPt,GenEta,GenPhi);
01363         }
01364         else {
01365            AddParticleToList(GenHfs,iPartMC,IDGELE,GenE,GenPt,GenEta,GenPhi); 
01366         } 
01367      }
01368      // Get stable muon quantitites
01369      else if (TMath::Abs(ModsPartMC[iPartMC]->GetPDG()) == 13){
01370 //        cout << "particle " << ModsPartMC[iPartMC]->GetPDG()<< endl;
01371         AddParticleToList(GenMuon,iPartMC,IDGMU,GenE,GenPt,GenEta,GenPhi); 
01372      }
01373      else {
01374         AddParticleToList(GenHfs,iPartMC,0,GenE,GenPt,GenEta,GenPhi); 
01375 //        cout << PartMC[iPartMC]->GetPDG() << endl;
01376 //        PartMC[iPartMC]->GetTParticlePDG()->Print();
01377      }
01378    }else{//--- for unstable generated ones
01379 //      cout << "unstable particle " << ModsPartMC[iPartMC]->GetPDG()<< endl;
01380       
01381        GenE = ModsPartMC[iPartMC]->GetE();
01382        GenTheta = ModsPartMC[iPartMC]->GetTheta();
01383        GenPt = ModsPartMC[iPartMC]->GetPt();
01384        if(GenE==0. || GenTheta==0. || GenPt==0.) continue; //--- remove part with Pt=0
01385        GenEta = -log(TMath::Tan((GenTheta)/2));
01386        GenPhi = ModsPartMC[iPartMC]->GetPhi();
01387 
01388      // Get tau quantitites
01389      if (TMath::Abs(ModsPartMC[iPartMC]->GetPDG()) == 15&&ModsPartMC[iPartMC]->GetStatus()!=3){
01390 //      printf("Tau found status=%d !\n",ModsPartMC[iPartMC]->GetStatus());
01391         AddParticleToList(GenTau,iPartMC,IDGTAU,GenE,GenPt,GenEta,GenPhi); 
01392 //              printf("Nb Tau %d!\n",GenTau->GetEntries());
01393      }
01394    }
01395   }//--- end loop over generated particles
01396  
01397 } 
01398 
01399 
01400   
01401 Int_t TMarEvent::AddParticleToList(TClonesArray *List,Int_t IndexID,Int_t ID,Float_t E,Float_t Pt,Float_t Eta,Float_t Phi)
01402 {
01403 
01404 // Add particle to list of preselected bodies ListOfUnSortedBodies
01405 // Return position/index in list
01406 
01407   // NumBodies = 0 if no entry;
01408   Int_t NumBodies = List->GetLast()+1;
01409   
01410   // Store index to particle candidate in array and create body
01411   
01412   // Store index to particle candidate in array and create body
01413   new (List->AddrAt(NumBodies)) TMarBody(ID,IndexID);
01414   
01415   TMarBody *body = (TMarBody *)List->AddrAt(NumBodies);
01416   body->SetEPtEtaPhi(E,Pt,Eta,Phi);
01417 
01418   
01419   return(NumBodies);
01420 
01421 }
01422 Int_t TMarEvent::AddParticleToList(TClonesArray *List,TMarBody *body)
01423 {
01424 
01425 // Add particle to list of preselected bodies ListOfUnSortedBodies
01426 // Return position/index in list
01427 
01428   // NumBodies = 0 if no entry;
01429   Int_t NumBodies = List->GetLast()+1;
01430   
01431   
01432   // Store index to particle candidate in array and create body
01433   new (List->AddrAt(NumBodies)) TMarBody(body);
01434   
01435   
01436   return(NumBodies);
01437 
01438 }
01439 Bool_t TMarEvent::Next(TRunHisto& run)
01440 {
01441 
01442 
01443 //... set  the next requested syst. scale 
01444   Syst_Type=run.Syst.Next();
01445   Syst_Sign=run.Syst.Sign;
01446 
01447 //... if NOSCALE, initialize the new event
01448   if(run.Syst.Type==0) {
01449      if(!Init()){
01450 
01451 //... write on file before going out of the event loop
01452        run.output_file_root->Write();
01453        return kFALSE;
01454      }
01455 
01456   }
01457 //... otherwise clear lists of selected particles only
01458   else{
01459      ClearSelectedParticles();
01460   }
01461 
01462 
01463 //... then calculate systematics (if requested)
01464   ApplySystematicShifts(run.Syst);
01465 
01466 //   printf("Event %d Sys_type=%d Sys_Num=%d\n",num,Syst_Type,Syst_Num);
01467 
01468 //... calculate particle dependent variables
01469   CalculateDerivedVariables();
01470 
01471 //... define more strict particles
01472   SelectParticles();
01473 
01474 //... Trigger efficiencies
01475   ApplyTriggerEfficiencies();
01476 
01477 //...  set the new output directory
01478   if(!run.output_file_root->cd(run.Syst.Name.c_str())){
01479     cout << "TMarEvent::Next= error in cd to "<< run.Syst.Name<< endl;
01480     return kFALSE;
01481 
01482   }
01483 
01484   return kTRUE;
01485 
01486 
01487 
01488 }
01489 
01490 Bool_t TMarEvent::Init()
01491 {
01492 
01493 //...Init event
01494 //...there is a local Event loop with cuts for runlist, HV, filling
01495 
01496   Bool_t initOk = kFALSE;
01497   Int_t  badLumi = 0;
01498   Int_t  badHV = 0;
01499   Int_t  badFillGen = 0;
01500 
01501   static H1ShortPtr dstType("RunType");
01502 
01503   while ( !initOk ) {
01504 
01505     if ( num == numMax ) {
01506       cout << "--> TMarEvent::Init : Reached the maximum number of events" << endl;
01507       return false;
01508     }
01509 
01510     if ( !gH1Tree->Next() ) {
01511       cout << "--> TMarEvent::Init : gH1Tree->Next()= false" << endl;
01512       return false;
01513     }
01514 
01515 //...Init type of run
01516     MCFlag = *dstType;
01517 
01518 //...re-set run and event number
01519 
01520     if ( RunNumber != gH1Tree->GetRunNumber() ) { NewRun = kTRUE; } else { NewRun = kFALSE; }
01521 
01522 //... increment event counter
01523     num++;
01524 
01525     RunNumber = gH1Tree->GetRunNumber();
01526     EventNumber = gH1Tree->GetEventNumber();
01527 
01528     //cout << RunNumber << " ev " << EventNumber << endl;
01529 
01530     if ( RunNumber == 0 ) {
01531       cout << "--> TMarEvent::Init = RunNumber=0!" << endl;
01532 
01533       Char_t ctemp[300];
01534       gH1Tree->GetCurFilename(3,ctemp);
01535       cout << "--> current file HAT " << ctemp << endl;
01536       gH1Tree->GetCurFilename(2,ctemp);
01537       cout << "--> current file MODS " << ctemp << endl;
01538 
01539       return(kFALSE);
01540     }
01541 
01542 //... if lumi is given
01543     if ( RunType == Data && lumi ) {
01544 
01545 //... in case of new run verify whether the run is selected
01546 
01547       //if ( NewRun == kTRUE ) {  cout << RunNumber << endl; }
01548 
01549         if ( !lumi->Update(RunNumber) ) {
01550           badLumi++;
01551           continue;
01552         }
01553 
01554 //... verify whether the HV condition are ok, otherwise go to the next event
01555         if ( !lumi->detStatus->IsOn() ) {
01556           badHV++;
01557           continue;
01558         }
01559     }
01560 
01561 //...clear event
01562    Clear();
01563 
01564 //... fills event variables from generated quantities and test generated cuts
01565 //... go to the next event if the present is out of cuts
01566 
01567     if ( !FillGeneratedParticlesAndCut() ) {
01568       badFillGen++;
01569       continue;
01570     }
01571 
01572     initOk = kTRUE;
01573 //...end local event loop
01574   }
01575 
01576   if ( badLumi ) { Info( "TMarEvent::Init 01", "Skipped %d events (not in runlist)", badLumi); }
01577   if ( badHV ) {   Info( "TMarEvent::Init 02", "Skipped %d events (bad HV settings)", badHV); }
01578   if ( badFillGen ) { Info( "TMarEvent::Init 03", "Gen. Filling for last %d events not passed ", badFillGen); }
01579 
01580 //...Fill first from hat
01581    if ( !FillFromHat() ) {
01582      cout << "--> TMarEvent::Init = something wrong in FillFromHat" << endl;
01583      return kFALSE;
01584    }
01585 
01586 
01587 //... define input (eventually calibrated) particles (PartCalxxx)
01588   DefineInputParticles();
01589 
01590   return(kTRUE);
01591 }
01592 
01593 
01594 void TMarEvent::ApplySystematicShifts(TMarSyst& Syst)
01595 {
01596   if(TMAREVENTDEBUG) printf("Enter ApplySystematicShifts ...\n");
01597 
01598 
01599 
01600 //... copy the particles and return if no scale is requested
01601  
01602   PartEm=&PartCalEm;
01603   PartMuon=&PartCalMuon;
01604   PartJet=&PartCalJet;
01605   HadJetVec=HadCalJetVec;
01606   HadNoJetVec=HadCalNoJetVec;
01607 
01608 //... return if no scale to be done
01609   if(Syst_Type==0) {
01610         Syst_Num=0; //--- init again number of systematic to be proceed to 0
01611         return;
01612   }
01613   
01614   Syst_Num+=1; //--- increment systematic number
01615   
01616 //... otherwise perform the scale 
01617   PartScaled[0].clear();  
01618   PartScaled[1].clear();  
01619   PartScaled[2].clear();  
01620   
01621 
01622   switch (Syst_Type) 
01623   {
01624      case 1:  // electron Energy
01625      {
01626         PartScaled[0]=PartCalEm;
01627         Syst.ScaleEnergyPartEm(PartScaled[0],Vertex[2],RunYear,0);
01628         PartEm=&PartScaled[0]; 
01629         break;
01630      }
01631      case 16:  // electron Energy correlated
01632      {
01633         PartScaled[0]=PartCalEm;
01634         Syst.ScaleEnergyPartEm(PartScaled[0],Vertex[2],RunYear,1);
01635         PartEm=&PartScaled[0]; 
01636         break;
01637      }
01638      case 17:  // electron Energy uncorrelated
01639      {
01640         PartScaled[0]=PartCalEm;
01641         Syst.ScaleEnergyPartEm(PartScaled[0],Vertex[2],RunYear,2);
01642         PartEm=&PartScaled[0]; 
01643         break;
01644      }
01645      case 2:  // electron Theta
01646      {
01647         PartScaled[0]=PartCalEm;
01648         Syst.ScaleThetaPartEm(PartScaled[0],RunYear);
01649         PartEm=&PartScaled[0];    
01650         break;
01651      }
01652      case 3:  // Jet Energy
01653      {
01654         PartScaled[0]=PartCalJet;
01655         Syst.ScaleEnergyPartJet(PartScaled[0],HadJetVec,HadNoJetVec,RunYear,0);
01656         PartJet=&PartScaled[0];    
01657         break;
01658      }
01659      case 18:  // Jet Energy correlated
01660      {
01661         PartScaled[0]=PartCalJet;
01662         Syst.ScaleEnergyPartJet(PartScaled[0],HadJetVec,HadNoJetVec,RunYear,1);
01663         PartJet=&PartScaled[0];    
01664         break;
01665      }
01666       case 19:  // Jet Energy uncorrelated
01667      {
01668         PartScaled[0]=PartCalJet;
01669         Syst.ScaleEnergyPartJet(PartScaled[0],HadJetVec,HadNoJetVec,RunYear,2);
01670         PartJet=&PartScaled[0];    
01671         break;
01672      }
01673     case 4:  // Jet Theta
01674      {
01675         PartScaled[0]=PartCalJet;
01676         Syst.ScaleThetaPartJet(PartScaled[0],HadJetVec);
01677         PartJet=&PartScaled[0];    
01678         break;
01679      }          
01680      case 5:  // Muon Pt
01681      {
01682         PartScaled[0]=PartCalMuon;
01683         Syst.ScalePtPartMuon(PartScaled[0]);
01684         PartMuon=&PartScaled[0];    
01685         break;
01686      }
01687      case 6:  // Muon Theta
01688      {
01689         PartScaled[0]=PartCalMuon;
01690         Syst.ScaleThetaPartMuon(PartScaled[0],RunYear);
01691         PartMuon=&PartScaled[0];    
01692         break;
01693      } 
01694      case 7: // electron, Jet, Muon Theta
01695      {
01696        
01697         PartScaled[0]=PartCalEm;
01698         Syst.ScaleThetaPartEm(PartScaled[0],RunYear);
01699         PartEm=&PartScaled[0]; 
01700         
01701         PartScaled[1]=PartCalJet;
01702         Syst.ScaleThetaPartJet(PartScaled[1],HadJetVec);
01703         PartJet=&PartScaled[1];
01704 
01705         PartScaled[2]=PartCalMuon;
01706         Syst.ScaleThetaPartMuon(PartScaled[2],RunYear);
01707         PartMuon=&PartScaled[2];     
01708  
01709         break;
01710      }    
01711      case 8: // electron et Jet Energy, Muon Pt
01712      {
01713         PartScaled[0]=PartCalEm;
01714         Syst.ScaleEnergyPartEm(PartScaled[0],Vertex[2],RunYear);
01715         PartEm=&PartScaled[0];  
01716 
01717         PartScaled[1]=PartCalJet;
01718         Syst.ScaleEnergyPartJet(PartScaled[1],HadJetVec,HadNoJetVec,RunYear);
01719         PartJet=&PartScaled[1];
01720 
01721         PartScaled[2]=PartCalMuon;
01722         Syst.ScalePtPartMuon(PartScaled[2]);
01723         PartMuon=&PartScaled[2];  
01724     
01725         break;
01726      }
01727      case 15:  // Noise effect
01728      {
01729         Syst.ScaleNoiseHFS(HadNoJetVec);
01730         break;
01731      }
01732     
01733   }  
01734   return;
01735 
01736 }
01737 
01738 void TMarEvent::DefineInputParticles()
01739 {
01740 
01741   Float_t PtPartNow,EtaPartNow,PhiPartNow,EPartNow;
01742 
01743 //...electrons (no calibration for the moment):
01744 
01745   for(Int_t iEm=0;iEm<ModsPartEm.GetEntries();iEm++){
01746     EPartNow = ModsPartEm[iEm]->GetE();
01747     PtPartNow  = ModsPartEm[iEm]->GetPt();
01748     if(PtPartNow==0){
01749       EtaPartNow  = -log(TMath::Tan((ModsPartEm[iEm]->GetTheta())/2));
01750       cout << "--> TMarEvent::DefineInputParticles : Warning, electron with Pt=0!"<< endl;
01751     }
01752     else EtaPartNow  = ModsPartEm[iEm]->GetEta();
01753     
01754     PhiPartNow  = ModsPartEm[iEm]->GetPhi();
01755 
01756 //... if electron calibrations: apply them    
01757       if( ElecCalib!=NULL && (ModsPartEm[iEm]->GetType()==1||ModsPartEm[iEm]->GetType()==11)) {
01758             TLorentzVector calibratedElectronCluster(0,0,0,0);
01759             TLorentzVector electronTrack(0,0,0,0);
01760             TLorentzVector uncalibratedElectronCluster(0,0,0,0);
01761             Int_t linkedTrackType=ModsPartEm[iEm]->GetTrType();  // H1PartEm::GetTrType() - LW, DTRA etc.
01762                 
01763             uncalibratedElectronCluster.SetPtEtaPhiE(ModsPartEm[iEm]->GetEnUncalib()*sin(ModsPartEm[iEm]->GetTheta()),
01764                         EtaPartNow,PhiPartNow,ModsPartEm[iEm]->GetEnUncalib());
01765             electronTrack.SetPtEtaPhiE(ModsPartEm[iEm]->GetTrMom()*sin(ModsPartEm[iEm]->GetTrTheta()),
01766                         ModsPartEm[iEm]->GetTrTheta(),ModsPartEm[iEm]->GetTrPhi(),ModsPartEm[iEm]->GetTrMom());
01767             //---- specific case for DTNV, for 2.5.x releases
01768             if(ModsPartEm[iEm]->GetNbNonVertFitTr()>0&&ModsPartEm[iEm]->GetNbSelTr()==0
01769                 &&ModsPartEm[iEm]->GetNbVertFitTr()==0) {
01770                         linkedTrackType=3;
01771                         electronTrack.SetPtEtaPhiE(ModsPartEm[iEm]->GetTrNvMom()*sin(ModsPartEm[iEm]->GetTrNvTh()),
01772                         ModsPartEm[iEm]->GetTrNvTh(),ModsPartEm[iEm]->GetTrNvPh(),ModsPartEm[iEm]->GetTrNvMom());
01773             }           
01774             Int_t charge=(Int_t)ModsPartEm[iEm]->GetTrCharge(); // of the electron track, H1PartEm::GetTrCharge()
01775 
01776             calibratedElectronCluster = ElecCalib->GetElecCalibration(uncalibratedElectronCluster,
01777                                     electronTrack,
01778                                     charge,Vertex[2],MCFlag,(int)RunYear,
01779                                     linkedTrackType);
01780 //          printf("PtUncalib=%f Ptcalib9900=%f Ptcalib=%f \n",ModsPartEm[iEm]->GetEnUncalib()*sin(ModsPartEm[iEm]->GetTheta())
01781 //                              ,PtPartNow,calibratedElectronCluster.Pt());
01782             EPartNow=calibratedElectronCluster.E();
01783             PtPartNow=calibratedElectronCluster.Pt();
01784             EtaPartNow=-log(tan(calibratedElectronCluster.Theta()/2));
01785             PhiPartNow=calibratedElectronCluster.Phi();
01786       }else if(H1ElecCalib!=NULL && (ModsPartEm[iEm]->GetType()==1||ModsPartEm[iEm]->GetType()==11)) {
01787       //--- use of H1ElecCalibrations
01788             TLorentzVector calibratedElectronCluster(0,0,0,0);
01789             TLorentzVector electronTrack(0,0,0,0);
01790             TLorentzVector uncalibratedElectronCluster(0,0,0,0);
01791             Int_t linkedTrackType=ModsPartEm[iEm]->GetTrType();  // H1PartEm::GetTrType() - LW, DTRA etc.
01792                 
01793             uncalibratedElectronCluster.SetPtEtaPhiE(ModsPartEm[iEm]->GetEnUncalib()*sin(ModsPartEm[iEm]->GetTheta()),
01794                         EtaPartNow,PhiPartNow,ModsPartEm[iEm]->GetEnUncalib());
01795             electronTrack.SetPtEtaPhiE(ModsPartEm[iEm]->GetTrMom()*sin(ModsPartEm[iEm]->GetTrTheta()),
01796                         ModsPartEm[iEm]->GetTrTheta(),ModsPartEm[iEm]->GetTrPhi(),ModsPartEm[iEm]->GetTrMom());
01797             //---- specific case for DTNV, for 2.5.x releases
01798             if(ModsPartEm[iEm]->GetNbNonVertFitTr()>0&&ModsPartEm[iEm]->GetNbSelTr()==0
01799                 &&ModsPartEm[iEm]->GetNbVertFitTr()==0) {
01800                         linkedTrackType=3;
01801                         electronTrack.SetPtEtaPhiE(ModsPartEm[iEm]->GetTrNvMom()*sin(ModsPartEm[iEm]->GetTrNvTh()),
01802                         ModsPartEm[iEm]->GetTrNvTh(),ModsPartEm[iEm]->GetTrNvPh(),ModsPartEm[iEm]->GetTrNvMom());
01803             }           
01804             Int_t charge=(Int_t)ModsPartEm[iEm]->GetTrCharge(); // of the electron track, H1PartEm::GetTrCharge()
01805 
01806             TVector3 Vtx(Vertex[0],Vertex[1],Vertex[2]);
01807             Float_t CalFactor = H1ElecCalib->GetElecCalibration(uncalibratedElectronCluster,
01808                                     electronTrack,
01809                                     charge,Vtx,RunYear,linkedTrackType);
01810             EPartNow=uncalibratedElectronCluster.E()*CalFactor;
01811             PtPartNow=uncalibratedElectronCluster.Pt()*CalFactor;
01812             EtaPartNow=-log(tan(uncalibratedElectronCluster.Theta()/2));
01813             PhiPartNow=uncalibratedElectronCluster.Phi();
01814       }
01815 
01816     PartCalEm.push_back(TMarBody(IDELE,iEm,EPartNow,PtPartNow,EtaPartNow,PhiPartNow));
01817   }
01818   
01819 //...muons: 
01820   for (Int_t iMu = 0;iMu<ModsPartMuon->GetEntries();iMu++) {
01821     EPartNow = ModsPartMuon[iMu]->GetE();
01822     PtPartNow  = ModsPartMuon[iMu]->GetPt();
01823     if(PtPartNow==0||EPartNow>1e8){
01824       EtaPartNow  = -log(TMath::Tan((ModsPartMuon[iMu]->GetTheta())/2));
01825 //      cout << "--> TMarEvent::DefineInputParticles : Warning, Muon with Pt=0!"<< endl;
01826     }
01827     else EtaPartNow  = ModsPartMuon[iMu]->GetEta();
01828     PhiPartNow  = ModsPartMuon[iMu]->GetPhi();
01829     
01830     if(EPartNow<1e8||EPartNow<-1e8) 
01831         PartCalMuon.push_back(TMarBody(IDMU,iMu,EPartNow,PtPartNow,EtaPartNow,PhiPartNow)); 
01832 //    else cout << "--> TMarEvent::DefineInputParticles : Warning, Muon with E=inf !"<< endl;
01833   }
01834   
01835  
01836 //... jets:
01837   
01838   H1ExclHfsIterator ModsPartHFS;
01839   if( JetCalib!=NULL ) {
01840 
01841 //. Apply Calibration if requested: calculate the calibrated HadCalJetVec and HadCalTotVec and fill PartCalJet
01842     ApplyHfsCalibration(); 
01843   
01844   } else if( HfsCalib!=NULL ) {
01845  //. Apply low Q2 Calibration if requested: calculate the calibrated HadCalJetVec and HadCalTotVec and fill PartCalJet
01846     ApplyLowQ2HfsCalibration(); 
01847   
01848   }else { 
01849 
01850     TLorentzVector HadCalTotVec(0,0,0,0);
01851 
01852 //. otherwise copy the ModsPartJet to PartCalJet and calculate the corresponding HadCalJetVec and HadNoJetVec 
01853     for(Int_t iJets=0;iJets<ModsPartJet->GetEntries();iJets++){
01854 
01855       // Get pointer to jet 
01856       H1PartJet *Jet = ModsPartJet[iJets];
01857 
01858       // Get kinematics of jet 
01859 
01860       PtPartNow = Jet->GetPt();
01861       if(PtPartNow==0){
01862         EtaPartNow  = -log(TMath::Tan((ModsPartJet[iJets]->GetTheta())/2));
01863         cout << "--> TMarEvent::DefineInputParticles : Warning, Jet with Pt=0!"<< endl;
01864       }
01865       else EtaPartNow  = ModsPartJet[iJets]->GetEta();
01866       EtaPartNow = Jet->GetEta();
01867       EPartNow = Jet->GetE();
01868       
01869       PhiPartNow = Jet->GetPhi();
01870       HadCalJetVec+=Jet->GetFourVector();
01871  
01872       PartCalJet.push_back(TMarBody(IDJET,iJets,EPartNow,PtPartNow,EtaPartNow,PhiPartNow)); 
01873     }
01874     for(Int_t ihad=0;ihad<ModsPartHFS.GetEntries();ihad++){
01875       
01876        HadCalTotVec += ModsPartHFS[ihad]->GetFourVector();
01877     }
01878 
01879     HadCalNoJetVec=HadCalTotVec-HadCalJetVec;
01880 
01881   }
01882     
01883 
01884 };
01885 
01886 void TMarEvent::ApplyHfsCalibration()
01887 {
01888 
01889  Int_t nJets=ModsPartJet.GetEntries();
01890 // if(nJets<=0) return;
01891 
01892 
01893  Float_t fPthCalib=0;
01894  Float_t fGammahCalib=0;
01895  Float_t fEmpzhCalib=0;
01896 
01897  TLorentzVector* jetscal= new TLorentzVector[nJets];
01898 
01899  TLorentzVector HadCalTotVec=JetCalib->GetHadroo2KtJetCalibration(fPthCalib,fGammahCalib,fEmpzhCalib,jetscal);
01900  for(int iJets=0;iJets<nJets;iJets++){
01901     TLorentzVector* jet=(jetscal+iJets);
01902     PartCalJet.push_back(TMarBody(IDJET,iJets,jet->E(),jet->Pt(),jet->Eta(),jet->Phi()));
01903     HadCalJetVec+=(*jet);
01904  }
01905  
01906  HadCalNoJetVec=HadCalTotVec-HadCalJetVec;
01907 
01908  
01909  delete [] jetscal;
01910 
01911 }
01912 
01913 void TMarEvent::ApplyLowQ2HfsCalibration()
01914 {
01915 //--- Apply the low Q2 hfs calibration 
01916 
01917  HfsCalib->ApplyHfsCalibration();    //--- apply calibrations once 
01918 
01919  //--- 4-vector of the calibrated Hadronic Final State (exclude isolated leptons)       
01920  TLorentzVector CalHFS=HfsCalib->GetCalibratedHFSFourVector();                          
01921  //--- array of 4-vector of all calibrated jets                                                 
01922  TObjArray* CalJetList=(TObjArray*) HfsCalib->GetCalibratedJets();                              
01923 
01924  for(Int_t iJets=0; iJets< CalJetList->GetEntries(); iJets++){                                    
01925     TLorentzVector* jet =  (TLorentzVector*) CalJetList->UncheckedAt(iJets);              
01926 //      printf("PtCal jet =%f \n",jetvect->Pt());                                                 
01927     PartCalJet.push_back(TMarBody(IDJET,iJets,jet->E(),jet->Pt(),jet->Eta(),jet->Phi()));
01928     HadCalJetVec+=(*jet);
01929  }                      
01930                                                                   
01931  TLorentzVector HadCalTotVec=CalHFS;
01932  HadCalNoJetVec=HadCalTotVec-HadCalJetVec;
01933 
01934 }
01935 
01936 void TMarEvent::GetPthEhTracksClusters(double& PthTracks,double& EhTracks,double& PthClusters,double& EhClusters){
01937 
01938    TLorentzVector HadronVecTracks(0,0,0,0);
01939    TLorentzVector HadronVecClusters(0,0,0,0);
01940    for(Int_t ihad=0;ihad<ModsPartHFS.GetEntries();ihad++){
01941       if(ModsPartHFS[ihad]->IsHFSChargedGoodTrk())HadronVecTracks += ModsPartHFS[ihad]->GetFourVector();
01942       else HadronVecClusters += ModsPartHFS[ihad]->GetFourVector();    
01943    } 
01944 
01945 
01946    PthTracks=HadronVecTracks.Pt();
01947    PthClusters=HadronVecClusters.Pt();
01948    EhTracks=HadronVecTracks.E();
01949    EhClusters=HadronVecClusters.E();   
01950 
01951 }
01952 
01953 
01954 void TMarEvent::ClearSelectedParticles()
01955 {
01956 // Clear up all list of selected particles
01957 
01958 //   electron->Clear();
01959 //   photon->Clear();
01960 //   muon->Clear();
01961 //   jet->Clear();
01962 //   np->Clear();
01963   electron->Delete();
01964   photon->Delete();
01965   muon->Delete();
01966   jet->Delete();
01967   np->Delete();
01968 
01969 }
01970 
01971 void TMarEvent::Clear()
01972 {
01973 // Clear up all allocated lists (ListOfBodies...)
01974 // Clear all event flags up
01975 
01976   electron->Delete();
01977   photon->Delete();
01978   muon->Delete();
01979   jet->Delete();
01980   np->Delete();
01981 
01982   GenElectron->Delete();
01983   GenPhoton->Delete();
01984   GenMuon->Delete();
01985   GenTau->Delete();
01986   GenNp->Delete();
01987   GenHfs->Delete();
01988 
01989   PartCalEm.clear();
01990   PartCalMuon.clear();
01991   PartCalJet.clear();
01992 
01993   PartEm=NULL;  
01994   PartMuon=NULL;  
01995   PartJet=NULL;  
01996 
01997   Vertex[0]=0;
01998   Vertex[1]=0;
01999   Vertex[2]=0;
02000  
02001   Epz=0;     
02002   Ptmiss=0;
02003   Et=0; 
02004 
02005   IndexElScat=-1;  
02006   Pte=0;   
02007   Ee=0;
02008   Pze=0;
02009   The=0;
02010   Phe=0;
02011   Q2e=0;
02012   ye=0;
02013   xe=0;
02014 
02015   Pth=0;    
02016   Thh=0;
02017   Phh=0; 
02018   Eh=0;
02019   Pzh=0;
02020   Q2h=0;
02021   yh=0;
02022   xh=0;
02023 
02024   Ptda=0;  
02025   Q2da=0;
02026   yda=0;
02027   xda=0;
02028 
02029   xes=0;
02030   yes=0;
02031   Q2es=0;
02032 
02033   tCJC=0;
02034   tLAr=0;
02035   Vparl=0;
02036   Vanti=0;
02037 
02038   polar=0;
02039   MCWeight=0;
02040 //  LumiWeight=1;
02041 
02042   HadJetVec.SetPxPyPzE(0,0,0,0);
02043   HadNoJetVec.SetPxPyPzE(0,0,0,0);
02044   HadCalNoJetVec.SetPxPyPzE(0,0,0,0);
02045   HadCalJetVec.SetPxPyPzE(0,0,0,0);
02046   TotalVec.SetPxPyPzE(0,0,0,0);
02047 
02048   GenScatEle.SetPxPyPzE(0,0,0,0);
02049   Syst_Type=0;
02050   Syst_Sign=0;
02051 
02052   RejectReason.clear();
02053 
02054 
02055 }
02056 
02057 
02058 
02059 Float_t TMarEvent::ApplyCCBackgroundCuts()
02060 {
02061 //---- All anti-background cuts applied in CC analysis
02062 //---- Also correct MC for additional inefficiency introduced in data
02063 
02064 //---- Timing cut
02065  if(!IsTimingOk()) return 0.;
02066 
02067  static H1IntPtr Ibg ("Ibg");     //standard array of background finders
02068  Int_t ibg=*Ibg;
02069  if (ibg!=0) {
02070         return 0.;
02071  }
02072 
02073 
02074 //---- Now Ibgfm, do not use bits 5,6,8
02075   static  H1IntPtr Ibgfm ("Ibgfm"); // Add array of background finders
02076   for(Int_t iBit=0;iBit<16;iBit++){
02077     if (((*Ibgfm >> iBit) & 1) == 1){
02078       //  if (iBit!=5 && iBit!=6 && iBit !=8){
02079         return 0.;
02080     }
02081   }
02082 
02083 //--- Now Ibgam, do not use bits 3,4,8
02084   static  H1IntPtr Ibgam ("Ibgam");
02085   for(Int_t iBit=0;iBit<9;iBit++){
02086     if (((*Ibgam >> iBit) & 1) == 1){
02087        if (iBit!=3 && iBit!=4 && iBit !=8){
02088              return 0.;
02089        }
02090     }
02091   }
02092 
02093 //--- coherent noise (inefficiency 1/75000)
02094   static  H1IntPtr Iqn ("Iqn");     //standard array of background finders
02095   Int_t iqn=*Iqn;
02096   if (iqn==10){
02097         return 0.;
02098   }
02099 
02100 
02101 
02102 //--- new bg finders cuts (Benjamin Portheault)
02103   static H1FloatPtr PtrLarClusTheta("HfsClusLarTheta");
02104   Float_t thetaLAr= (*PtrLarClusTheta)*MYR2D;
02105   Short_t nLW=ModsPartSelTrack.GetEntries();
02106   static H1ShortPtr ndtnv("NumNonvertexTracks");
02107   Short_t nDTNV=(*ndtnv);
02108   static H1ShortPtr ndtra("NumLooseTracks");
02109   Short_t nDTRA=(*ndtra);
02110 
02111   if (nLW==0 && MYR2D*Gammah > 40.) {
02112         return 0.;
02113   }
02114   if (nLW==0 && thetaLAr > 20.) {
02115         return 0.;
02116   }
02117   if (nDTRA!=0){
02118     if (nDTNV/nDTRA > 20) {
02119         return 0.;
02120     }
02121     if (nDTNV/nDTRA >=2 && nDTRA <3 && thetaLAr>40 ) {
02122         return 0.;
02123     }
02124   }
02125 
02126 // //--- now distance between any track and any cluster (from hat variables,2.7.7)
02127 //   if ( (dtrkcls>0.5 && thetacls*MYR2D >20) || dtrkcls > 1.) {
02128 //      return 0.; //---- !!! was disabled before ????
02129 //   }
02130 
02131 //--- Now for MC apply additional inefficiency
02132  if(MCFlag==0) return 1.;
02133  Float_t MCIneffWeight=1.;
02134 
02135 
02136  if(xh<=0 || Q2h<=0) return 1.0;
02137 
02138 
02139  Float_t logxh=log10(xh);
02140  Float_t logQ2h=log10(Q2h);
02141  Bool_t foundx=kFALSE;
02142  Bool_t foundq2=kFALSE;
02143  Int_t ibinx=7;
02144  Int_t ibinq2=7;
02145 //---- [x][y] y== colonne, x==ligne
02146  const Float_t bgtimingconst[8][8]={
02147  {0.971463,0.995788,1.0,1,1,1 ,1 ,1,},
02148  {0.986663,0.987707,0.989387,0.955399,1, 1, 1, 1},
02149  {0.996249,0.987958, 0.989286, 0.989491, 0.997307, 1, 1, 1},
02150  {0.991483, 0.985337, 0.985433, 0.977568, 0.997002, 0.970574, 1.0, 1},
02151  {0.936412, 0.987353, 0.981143, 0.985565, 0.972558, 0.980058, 0.992863, 1},
02152  {1, 1.0, 0.939658, 0.97497, 0.98069, 0.933589, 1.0, 0.960095},
02153  {1, 1, 1, 1, 1, 1, 1, 1},
02154  {1, 1, 1, 1, 1, 1, 1, 1}};
02155  const Float_t xl[8]={-2.33,-2.0,-1.67,-1.33,-1.0,-0.75,-0.5,-0.25};
02156  const Float_t xh[8]={-2.0,-1.67,-1.33,-1.0,-0.75,-0.5,-0.25,-0.0457};
02157  const Float_t q2l[8]={2.35,2.6,2.8,3.1,3.35,3.6,3.85,4.1};
02158  const Float_t q2h[8]={2.6,2.8,3.1,3.35,3.6,3.85,4.1,4.4};
02159 
02160  for(Int_t i=0;i<8;i++){
02161    if(xl[i] < logxh &&  logxh <= xh[i]) {
02162      ibinx=i;
02163      foundx=kTRUE;
02164    }
02165  if(q2l[i] < logQ2h &&  logQ2h <= q2h[i]) {
02166      ibinq2=i;
02167      foundq2=kTRUE;
02168    }
02169  }
02170 
02171  if(foundx==kFALSE) {
02172    if(logxh <= xl[0]) ibinx=0;
02173    if(logxh > xh[7]) ibinx=7;
02174  }
02175 
02176  if(foundq2==kFALSE) {
02177    if(logQ2h <= q2l[0]) ibinq2=0;
02178    if(logQ2h > q2h[7]) ibinq2=7;
02179  }
02180 
02181   MCIneffWeight=bgtimingconst[ibinx][ibinq2];
02182 
02183 
02184   return MCIneffWeight;
02185 }
02186 
02187 
02188 
02189 
02190 Float_t TMarEvent::CCTriggerEfficiency()
02191 {
02192 //--- To reweight for trigger CC efficiencies
02193 //    (hera-1 and hera-2 available)
02194 //    for hera-1, from isolated lepton analysis, but only 98-00 available ??
02195 Float_t TriggerWeight=1;
02196 
02197 
02198  if(IsDataHeraII() || (MCFlag&&RunYear==TMarEvent::Y0304) || (MCFlag&&RunYear==TMarEvent::Y05)){//--- hera2
02199    if(MCFlag) {
02200        TriggerWeight=fTrigReweight->ApplyCCTRIGGEFFreweightHERA2(xh,Q2h);
02201    }else{
02202      if( !(   (IsL1Trigg(66) && IsL4Trigg(66))
02203            || (IsL1Trigg(67)&& IsL4Trigg(67))
02204            || (IsL1Trigg(77) && IsL4Trigg(77)) )
02205         ) {//--- do not pass triggers
02206            Char_t towrite[200];
02207            sprintf(towrite,"Do not pass hera-2 CC triggers \n");
02208            RejectReason.push_back(towrite);
02209            TriggerWeight=0;
02210      }
02211    }
02212  }else{//---- hera-1 
02213    if(MCFlag) {
02214        TriggerWeight=fTrigReweight->CCTriggerEff(Pth,Thh,RunYear);
02215    }else{
02216      if(IsL4Trigg(66) == kFALSE 
02217        && IsL4Trigg(67) == kFALSE 
02218        && IsL4Trigg(71) == kFALSE 
02219        && IsL4Trigg(75) == kFALSE 
02220        && IsL4Trigg(77) == kFALSE
02221        &&( IsData2000() 
02222           || IsData1999p()
02223           || IsData1999e()
02224           || IsData1998() )
02225         ) {//--- do not pass triggers
02226            Char_t towrite[200];
02227            sprintf(towrite,"Do not pass hera-1 CC triggers \n");
02228            RejectReason.push_back(towrite);
02229            TriggerWeight=0;
02230      }
02231    }
02232  }
02233 
02234 
02235 return TriggerWeight;
02236 }
02237 
02238 
02239 Float_t TMarEvent::CCAndLeptonTriggerEfficiency()
02240 {
02241 //--- To reweight for trigger CC + leptons efficiencies (for isolated lepton analysis)
02242 Float_t TriggerWeight=1;
02243 Float_t CCTriggerWeight=1;
02244 Float_t MuTriggerWeight=1;
02245   
02246 //---- 100% efficiency if one electron E>1  
02247   if(ModsPartEm->GetEntries()>0&&ModsPartEm[0]->GetE()>10) return TriggerWeight;
02248 
02249 //--- Ptmiss triggers
02250   if(Ptmiss>8) CCTriggerWeight=CCTriggerEfficiency();
02251 //--- Muon triggers
02252   if(ModsPartMuon->GetEntries()>0&&ModsPartMuon[0]->GetGrade()<4) MuTriggerWeight=MuTriggerEfficiency();
02253 
02254 //--- combine both
02255   TriggerWeight=max(CCTriggerWeight,MuTriggerWeight);
02256   
02257   return TriggerWeight;
02258 }
02259 
02260 
02261 
02262 
02263 Float_t TMarEvent::JJTriggerEfficiency()
02264 {
02265 //--- To reweight for trigger JJ efficiencies
02266 //    (only determined for hera-1, all years. Set to 1 for hera-2)
02267 Float_t TriggerWeight=1;
02268 
02269    if(MCFlag) {
02270      TMarBody *Jet1 = (TMarBody*) jet->At(0);
02271      H1PartJet *ModsJet1 = ModsPartJet[Jet1->GetIndexID()];
02272      TriggerWeight=fTrigReweight->DiJetsTriggerEff(ModsJet1->GetPt(),ModsJet1->GetTheta(),RunYear);
02273    }else{
02274     if(IsL4Trigg(64) == kFALSE 
02275        && IsL4Trigg(67) == kFALSE 
02276        && IsL4Trigg(75) == kFALSE 
02277        && IsL4Trigg(77) == kFALSE
02278        && !IsDataHeraII()) {//--- do not pass triggers
02279            Char_t towrite[200];
02280            sprintf(towrite,"Do not pass hera-1 di-jet triggers \n");
02281            RejectReason.push_back(towrite);
02282            TriggerWeight=0;
02283        }
02284    }
02285 
02286 return TriggerWeight;
02287 }
02288 
02289 
02290 Float_t TMarEvent::MuTriggerEfficiency()
02291 {
02292 //--- To reweight for Muon trigger efficiencies
02293 Float_t TriggerWeight=1;
02294 
02295  if(IsDataHeraII() || (MCFlag&&RunYear==TMarEvent::Y0304) || (MCFlag&&RunYear==TMarEvent::Y05) ){//--- hera2
02296     if(IsL1Trigg(18) == kFALSE) TriggerWeight=0;
02297  }else{//--- hera1
02298     if(IsL1Trigg(19) == kFALSE && IsL1Trigg(22) == kFALSE
02299        && IsL1Trigg(34) == kFALSE  && IsL1Trigg(56) == kFALSE) TriggerWeight=0;
02300  }
02301 
02302 
02303 return TriggerWeight;
02304 }
02305 
02306 
02307 
02308 void TMarEvent::ApplyTriggerEfficiencies()
02309 {
02310 //--- To reweight for trigger efficiencies
02311 //---  Here is an exemple function ! to apply trigger efficiencies, 
02312 //     copy this function to your own event and un-comment the needed trigger efficiency
02313 if(Syst_Type!=0) return; //--- don't do it again for systematics
02314 
02315 Float_t TriggerWeight=1;
02316 
02317 //  TriggerWeight=JJTriggerEfficiency();
02318 //  TriggerWeight=CCTriggerEfficiency();
02319 //  TriggerWeight=MuTriggerEfficiency();
02320 //  TriggerWeight=CCAndLeptonTriggerEfficiency();
02321 
02322 //--- modify MCWeight
02323    fTriggerWeight=TriggerWeight;
02324    MCWeight*=TriggerWeight;
02325 
02326 return;
02327 
02328 }
02329 
02330 void TMarEvent::ApplyMCXSectionReweight()
02331 {
02332 //---- to reweight MC xsection to polarisation for hera2 NC/CC or NLO for epvec
02333 
02334  if(Syst_Type!=0)return; //--- don't apply for systematics
02335  if(RunType!=NC_Djo&&RunType!=CC_Djo) return;
02336  if(RunYear!=Y0304&&RunYear!=Y05)return;
02337 
02338  static H1FloatPtr prtQ2Gen("Q2hGen");
02339  static H1FloatPtr prtyGen("YhGen");
02340  const Double_t GenS = 4*GenPl*GenPp;
02341 
02342  Double_t Q2gen=(*prtQ2Gen);
02343  Double_t ygen=(*prtyGen);
02344  Double_t xgen=0;
02345 
02346  Short_t LeptonCharge=1;
02347  if(RunYear==Y05) LeptonCharge=-1;
02348 
02349  if(ygen!=0) xgen = Q2gen/(ygen*GenS);
02350 
02351 //--- NC/CC reweight only for hera2
02352 
02353 
02354 //---define the polarisation value
02355  if(lumi==NULL){
02356 
02357 // set value on the base of average polarition
02358       if(RunYear==Y0304) polar=-0.02;
02359       else if(RunYear==Y05) {polar=-0.20;}
02360 
02361  } else {
02362 
02363 // get the polarization from the data distribution
02364       polar=(Double_t)lumi->GetPolar();
02365 
02366  }
02367 
02368  if(RunType==CC_Djo) {
02369    MCWeight *= (1+polar*LeptonCharge);
02370  }
02371 
02372  //fQCDWeight=fMCReweight->GetQCDWeight(xgen,Q2gen,polar,LeptonCharge,(Short_t)RunType);
02373  //MCWeight*=fQCDWeight;
02374 
02375 }
02376 
02377 
02378 void TMarEvent::PrintEventInfo()
02379 {
02380 // Print informations on event
02381 //
02382 
02383 
02384  printf("*****************************************************\n");
02385  printf("Run/Event %d %d \n",RunNumber,EventNumber);
02386  
02387  Bool_t rejected=kFALSE;
02388  
02389  for(list<string>::iterator it=RejectReason.begin();it!=RejectReason.end(); it++) {
02390         string& reason=(*it);
02391         rejected=kTRUE;
02392         printf("Rejected: %s\n",reason.c_str());
02393  }
02394  
02395 // if(rejected) return;
02396  
02397 //---- not rejected: print all infos
02398   printf(" E-Pz = %f Ptmiss = %f \n",Epz,Ptmiss); 
02399   printf(" Pth = %f Gammah = %f \n",Pth,Gammah); 
02400   
02401 //---- print electron infos
02402   printf("Nb electron = %d Nb muons =%d Nb jets =%d Nb neutrino =%d photon = %d \n",
02403                 electron->GetEntries(),muon->GetEntries(),jet->GetEntries(),
02404                 np->GetEntries(),photon->GetEntries());
02405 
02406  for(Int_t itk=0;itk<electron->GetEntries();itk++) {
02407      TMarBody *Body = (TMarBody*) electron->At(itk);
02408      printf("El %d: Pt=%f E=%f Theta=%f Phi=%f \n",itk,Body->GetPt(),Body->GetE(),
02409                 Body->GetTheta()*MYR2D,Body->GetPhi()*MYR2D);
02410  }
02411  for(Int_t itk=0;itk<muon->GetEntries();itk++) {
02412      TMarBody *Body = (TMarBody*) muon->At(itk);
02413      printf("Mu %d: Pt=%f E=%f Theta=%f Phi=%f \n",itk,Body->GetPt(),Body->GetE(),
02414                 Body->GetTheta()*MYR2D,Body->GetPhi()*MYR2D);
02415  }
02416  for(Int_t itk=0;itk<jet->GetEntries();itk++) {
02417      TMarBody *Body = (TMarBody*) jet->At(itk);
02418      printf("Jet %d: Pt=%f E=%f Theta=%f Phi=%f \n",itk,Body->GetPt(),Body->GetE(),
02419                 Body->GetTheta()*MYR2D,Body->GetPhi()*MYR2D);
02420  }
02421  for(Int_t itk=0;itk<photon->GetEntries();itk++) {
02422      TMarBody *Body = (TMarBody*) photon->At(itk);
02423      printf("Pho %d: Pt=%f E=%f Theta=%f Phi=%f \n",itk,Body->GetPt(),Body->GetE(),
02424                 Body->GetTheta()*MYR2D,Body->GetPhi()*MYR2D);
02425  }
02426  for(Int_t itk=0;itk<np->GetEntries();itk++) {
02427      TMarBody *Body = (TMarBody*) np->At(itk);
02428      printf("Np %d: Pt=%f E=%f Theta=%f Phi=%f \n",itk,Body->GetPt(),Body->GetE(),
02429                 Body->GetTheta()*MYR2D,Body->GetPhi()*MYR2D);
02430  }
02431   
02432  
02433  
02434 
02435 }
02436 
02437   
02438 
02439