00001
00002 #include <XANADOO/domain/interface/XANAEsdEvent.h>
00003 #include <XANADOO/XANAMcInfo/interface/XANAGeneratorEvent.h>
00004 #include <XANADOO/XANAMcInfo/interface/XANAGeantEvent.h>
00005 #include <XANADOO/XANAMcInfo/interface/XANAGeneratorParticle.h>
00006
00007 #include <XANADOO/XANAMcInfo/interface/XANAGeantTrack.h>
00008 #include <XANADOO/XANAMcInfo/interface/XANAGeantVertex.h>
00009
00010 #include <XANADOO/XANAElectronCandidate/interface/XANAElectronGSFTrack.h>
00011 #include <XANADOO/XANAClusters/interface/XANACaloRecHit.h>
00012 #include <XANADOO/XANATracks/interface/XANATrackHit.h>
00013 #include <XANADOO/XANAClusters/interface/XANACluster.h>
00014 #include <XANADOO/XANAClusters/interface/XANAEmCluster.h>
00015 #include <XANADOO/XANAClusters/interface/XANASuperCluster.h>
00016 #include <XANADOO/XANAClusters/interface/XANACaloRecHit.h>
00017 #include <XANADOO/XANAElectronCandidate/interface/XANAElectronCandidate.h>
00018 #include <XANADOO/XANAVertex/interface/XANAVertex.h>
00019
00020 #include <CLHEP/Geometry/Point3D.h>
00021
00022 #include "TFile.h"
00023 #include "TBranch.h"
00024 #include "TObjArray.h"
00025 #include "TRefArray.h"
00026 #include "TTree.h"
00027 #include "TChain.h"
00028
00029 #include <XANADOO/XANAnalysisZZ/interface/treeStruct.h>
00030
00031 #include <map>
00032 #include <vector>
00033
00034 #define PI 3.14159265358979323846
00035
00036 #define SIGMA_EB 0.04
00037 #define SIGMA_EE 0.150
00038
00039 using namespace std;
00040
00041 void globalReset(simpleEvent *e);
00042 void kineReset(simpleGeneratorEvent *k);
00043
00044 int mom(int idx, TObjArray *tracks);
00045 XANAGeantVertex *getDaughterVertex(XANAGeantTrack *gtk, TObjArray *tracks);
00046
00047
00048 void setAddr(TTree *T, XANAEsdEvent *event, XANAGeneratorEvent *generatorEvent, XANAGeantEvent *geantEvent)
00049 {
00050 T->SetBranchAddress("Event.",&event);
00051 T->SetBranchStatus("Event.",0);
00052 T->SetBranchAddress("GeneratorEvent.",&generatorEvent);
00053 T->SetBranchAddress("GeantEvent.",&geantEvent);
00054 }
00055
00056
00057
00058 int main(int argc, char** argv)
00059 {
00060 if (argc==1) {
00061 printf("Usage: elTree [XANADOO tree]\n");
00062 exit(1);
00063 }
00064
00065 XANAEsdEvent *event = new XANAEsdEvent();
00066 XANAGeneratorEvent *generatorEvent = new XANAGeneratorEvent();
00067 XANAGeantEvent *geantEvent = new XANAGeantEvent();
00068
00069 TFile *fout = new TFile("higgsTree.root","RECREATE");
00070
00071
00072 TTree *tree = new TTree("tree","tree");
00073
00074
00075
00076
00077 simpleEvent e;
00078 simpleGeneratorEvent k;
00079
00080 tree->Branch("k",&k.ievent,
00081 "ievent/I"
00082 ":ifirst/I"
00083 ":nple/I"
00084 ":pid[10]/I"
00085 ":mom[10]/I"
00086 ":x[10][4]/F"
00087 ":p[10][4]/F"
00088 ":pMag[10]/F"
00089 ":pT[10]/F"
00090 ":eta[10]/F"
00091 ":phi[10]/F"
00092 ":mH/F"
00093 ":mZ1/F"
00094 ":mZ2/F"
00095 );
00096 tree->Branch("e",&e.ievent,
00097 "ievent/I"
00098
00099 ":kn_deltaRKine/F"
00100 ":kn_pid/I"
00101 ":kn_pidMo1/I"
00102 ":kn_lineMo1/I"
00103 ":kn_pidMo2/I"
00104 ":kn_lineMo2/I"
00105 ":kn_E/F"
00106 ":kn_pT/F"
00107 ":kn_eta/F"
00108 ":kn_phi/F"
00109 ":kn_x[4]/F"
00110 ":kn_p[4]/F"
00111 ":kn_pMag/F"
00112
00113 ":sc_isc/I"
00114 ":sc_nsc/I"
00115 ":sc_Erec/F"
00116 ":sc_ErecUncorr/F"
00117 ":sc_ErecScaleFactor/F"
00118 ":sc_ErecS1/F"
00119 ":sc_ErecS4/F"
00120 ":sc_ErecS9/F"
00121 ":sc_ErecS25/F"
00122 ":sc_eta/F"
00123 ":sc_phi/F"
00124 ":sc_seedE/F"
00125 ":sc_seedEta/F"
00126 ":sc_seedPhi/F"
00127 ":sc_seedNxtals/I"
00128 ":sc_seedNxtals2sigma/I"
00129 ":sc_nBrems/I"
00130 ":sc_bremEtot/F"
00131 ":sc_bremNxtals/I"
00132 ":sc_bremNxtals2sigma/I"
00133
00134 ":tk_ntk/I"
00135 ":tk_itk/I"
00136
00137 ":tk_xetaFP/F"
00138 ":tk_xphiFP/F"
00139 ":tk_xFP[3]/F"
00140 ":tk_xetaLP/F"
00141 ":tk_xphiLP/F"
00142 ":tk_xLP[3]/F"
00143
00144 ":tk_petaFPmode/F"
00145 ":tk_pphiFPmode/F"
00146 ":tk_pFPmode[3]/F"
00147 ":tk_pFPmodeMag/F"
00148 ":tk_petaFPmean/F"
00149 ":tk_pphiFPmean/F"
00150 ":tk_pFPmean[3]/F"
00151 ":tk_pFPmeanMag/F"
00152
00153 ":tk_petaLPmode/F"
00154 ":tk_pphiLPmode/F"
00155 ":tk_pLPmode[3]/F"
00156 ":tk_pLPmodeMag/F"
00157 ":tk_petaLPmean/F"
00158 ":tk_pphiLPmean/F"
00159 ":tk_pLPmean[3]/F"
00160 ":tk_pLPmeanMag/F"
00161
00162 ":tk_petaVXmode/F"
00163 ":tk_pphiVXmode/F"
00164 ":tk_pVXmode[3]/F"
00165 ":tk_pVXmodeMag/F"
00166 ":tk_petaVXmean/F"
00167 ":tk_pphiVXmean/F"
00168 ":tk_pVXmean[3]/F"
00169 ":tk_pVXmeanMag/F"
00170
00171 ":tk_tkNhits/I"
00172 ":tk_tkNLostHits/I"
00173 ":tk_chi2OverDof/F"
00174 ":tk_chi2OverDofBack/F"
00175 ":tk_chargeRec/I"
00176 ":tk_highestWeightFP/F"
00177 ":tk_highestWeightLP/F"
00178 ":tk_highestWeightVX/F"
00179
00180 ":vx_x[3]/F"
00181 ":vx_xeDiag[3]/F"
00182 ":vx_xeOffDiag[3]/F"
00183 ":vx_chi2/F"
00184 ":vx_chi2Norm/F"
00185 ":vx_Ndof/I"
00186 ":vx_vectSum[3]/F"
00187
00188 ":el_iel/I"
00189 ":el_phi/F"
00190 ":el_eta/F"
00191 ":el_pVX[3]/F"
00192 ":el_pVXMag/F"
00193 ":el_pTVX/F"
00194 ":el_pEC[3]/F"
00195 ":el_pECMag/F"
00196 ":el_DetaScTkVx/F"
00197 ":el_DetaScTkEc/F"
00198 ":el_DetaSeedTkEc/F"
00199 ":el_DphiScTkVx/F"
00200 ":el_DphiScTkEc/F"
00201 ":el_DphiSeedTkEc/F"
00202 ":el_tkIso/F"
00203 ":el_ecIso/F"
00204 ":el_EOP/F"
00205 ":el_EseedOP/F"
00206 ":el_HOE/F"
00207
00208 ":tg_l1D/b"
00209 ":tg_hltD/b"
00210 ":tg_l1R[128]/b"
00211 ":tg_hltR[152]/b"
00212 );
00213
00214
00215 TTree *T;
00216 Int_t nevent=0;
00217 Int_t nb = 0;
00218 e.ievent=0;
00219 k.ievent=0;
00220 for (int i=1; i<argc; i++) {
00221 printf("--> %s\n",argv[i]);
00222 TFile *F = new TFile(argv[i]);
00223 if (!(F->GetNkeys())) {
00224 printf("File %s has no key --> skipping\n",argv[i]);
00225 continue;
00226 }
00227 T = (TTree*)F->Get("ESD");
00228 nevent = (Int_t)T->GetEntries();
00229 printf("Going to read %d events...\n",nevent);
00230 setAddr(T,event,generatorEvent,geantEvent);
00231
00232
00233
00234
00235
00236 TRefArray *gen=0;
00237 TRefArray *genAll=0;
00238 TObjArray *tracks=0;
00239 TObjArray *sc=0;
00240
00241 for (Int_t iev=0; iev<nevent; iev++) {
00242 nb += T->GetEntry(iev);
00243 ++e.ievent;
00244 ++k.ievent;
00245
00246 if (iev%50==1) fprintf(stdout," Event % 5d\n",iev);
00247
00248 tracks = geantEvent->getGeantTracks();
00249 gen = generatorEvent->getStableSignalParticles();
00250 genAll = generatorEvent->getSignalParticles();
00251
00252 kineReset(&k);
00253
00254
00255
00256 vector<int> plist;
00257 int ifirstel=0;
00258 int aZ=0;
00259
00260 for (int ip=0; ip<generatorEvent->getNumberOfStableSignalParticles(); ip++) {
00261 XANAGeneratorParticle *gp = (XANAGeneratorParticle*)(*gen)[ip];
00262 if (abs(gp->getPartId())!=11 && abs(gp->getPartId())!=22) continue;
00263 int mo1 = gp->getMo1();
00264 int line=0;
00265
00266 while (1) {
00267 XANAGeneratorParticle *gpMoN = (XANAGeneratorParticle*)(*genAll)[mo1-1];
00268 mo1 = gpMoN->getMo1();
00269 if ((gpMoN->getPartId()==23 && gpMoN->getStatus()==3)) {
00270 line = gpMoN->getLine();
00271 break;
00272 }
00273 if (mo1==0) break;
00274 }
00275 if (line) {
00276 if (ifirstel==0) {
00277 aZ = line;
00278 ifirstel=1;
00279 }
00280
00281 if (line==aZ) plist.push_back(ip);
00282 else plist.push_back(-ip);
00283 }
00284 }
00285 k.nple = (int)plist.size();
00286 float ptot[4];
00287 float pZ1[4];
00288 float pZ2[4];
00289 for (int j=0; j<4; j++) {
00290 ptot[j]=0;
00291 pZ1[j] =0;
00292 pZ2[j] =0;
00293 }
00294
00295 map<float,int> ptlist;
00296 for (unsigned int ip=0; ip<plist.size(); ip++) {
00297 XANAGeneratorParticle *gp = (XANAGeneratorParticle*)(*gen)[abs(plist[ip])];
00298 ptlist[gp->getMomentum().vect().perp()] = ip;
00299 }
00300
00301 map<float,int>::reverse_iterator rit;
00302 int cnt=0;
00303 for (rit=ptlist.rbegin(); rit!=ptlist.rend(); rit++) {
00304 int ip = rit->second;
00305 XANAGeneratorParticle *gp = (XANAGeneratorParticle*)(*gen)[abs(plist[ip])];
00306
00307 k.pid[cnt] = gp->getPartId();
00308 k.mom[cnt] = (int)copysign(1,plist[ip]);
00309 k.x[cnt][0] = gp->getVertex().x();
00310 k.x[cnt][1] = gp->getVertex().y();
00311 k.x[cnt][2] = gp->getVertex().z();
00312 k.x[cnt][3] = gp->getVertex().e();
00313 k.p[cnt][0] = gp->getMomentum().x();
00314 k.p[cnt][1] = gp->getMomentum().y();
00315 k.p[cnt][2] = gp->getMomentum().z();
00316 k.p[cnt][3] = gp->getMomentum().e();
00317 k.pMag[cnt] = gp->getMomentum().vect().mag();
00318 k.pT[cnt] = gp->getMomentum().vect().perp();
00319 k.eta[cnt] = gp->getMomentum().eta();
00320 k.phi[cnt] = gp->getMomentum().phi();
00321
00322 for (int j=0; j<4; j++) {
00323 ptot[j] += k.p[cnt][j];
00324 if (k.mom[cnt]>0) pZ1[j] += k.p[cnt][j];
00325 else pZ2[j] += k.p[cnt][j];
00326 }
00327 cnt++;
00328
00329
00330 }
00331 float mmH=0, mmZ1=0, mmZ2=0;
00332
00333 for (int i=0; i<3; i++) {
00334 mmH -= ptot[i]*ptot[i];
00335 mmZ1 -= pZ1[i]*pZ1[i];
00336 mmZ2 -= pZ2[i]*pZ2[i];
00337 }
00338 mmH += ptot[3]*ptot[3];
00339 mmZ1 += pZ1[3]*pZ1[3];
00340 mmZ2 += pZ2[3]*pZ2[3];
00341
00342 k.mH = sqrt(mmH);
00343 if (mmZ1>mmZ2) {
00344 k.mZ1 = sqrt(mmZ1);
00345 k.mZ2 = sqrt(mmZ2);
00346 } else {
00347 k.mZ2 = sqrt(mmZ1);
00348 k.mZ1 = sqrt(mmZ2);
00349 }
00350 plist.clear();
00351 ptlist.clear();
00352
00353
00354
00355
00356
00357
00358
00359 e.tg_l1D=0;
00360 e.tg_hltD=0;
00361 for (int i=0; i<128; i++) e.tg_l1R[i]=0;
00362 for (int i=0; i<152; i++) e.tg_hltR[i]=0;
00363
00364 XANATriggerInfo xtg = event->getTriggerInfo();
00365 e.tg_l1D = xtg.getGlobalL1Decision();
00366 e.tg_hltD = xtg.getHltDecision();
00367 cout << xtg.getGlobalL1Decision() << " " << xtg.getHltDecision() << " " << endl;
00368 TBits tb;
00369 tb = xtg.getGlobalL1Response();
00370 for (int i=0; i<16*8; i++) {
00371 e.tg_l1R[i] = tb.TestBitNumber(i);
00372 }
00373 tb = xtg.getHltResponse();
00374 for (int i=0; i<19*8; i++) {
00375 e.tg_hltR[i] = tb.TestBitNumber(i);
00376 }
00377
00378 int ifilled=0;
00379
00380 sc = event->getSuperClusters();
00381 e.sc_nsc = event->getNumberOfSuperClusters();
00382 e.sc_isc=0;
00383
00384
00385
00386 for (int iisc=0; iisc<e.sc_nsc; iisc++) {
00387
00388 globalReset(&e);
00389
00390
00391 XANASuperCluster *xsc = (XANASuperCluster*)(*sc)[iisc];
00392
00393
00394
00395
00396 HepPoint3D point = xsc->getPosition();
00397 point.setMag(xsc->getEnergy());
00398 HepLorentzVector vecLor( point.x(), point.y(), point.z(), xsc->getEnergy());
00399 map<float, int> deltaR;
00400 deltaR.clear();
00401 TRefArray *gen = generatorEvent->getStableSignalParticles();
00402 TObjArray *genAll = generatorEvent->getParticles();
00403 for (int ip=0; ip<generatorEvent->getNumberOfStableSignalParticles(); ip++) {
00404 XANAGeneratorParticle *genpar = (XANAGeneratorParticle*)(*gen)[ip];
00405 if (abs(genpar->getPartId())==11
00406 || genpar->getPartId()==22 &&
00407 genpar->getMomentum().e()/xsc->getEnergy()>.1)
00408 deltaR[vecLor.deltaR(genpar->getMomentum())] = ip;
00409 }
00410 XANAGeneratorParticle *genpar = (XANAGeneratorParticle*)(*gen)[deltaR.begin()->second];
00411 XANAGeneratorParticle *genparMo1 = (XANAGeneratorParticle*)(*genAll)[genpar->getMo1()-1];
00412 XANAGeneratorParticle *genparMo2 = (XANAGeneratorParticle*)(*genAll)[genparMo1->getMo1()-1];
00413
00414 e.kn_deltaR = deltaR.begin()->first;
00415 e.kn_x[0] = genpar->getVertex().x();
00416 e.kn_x[1] = genpar->getVertex().y();
00417 e.kn_x[2] = genpar->getVertex().z();
00418 e.kn_x[3] = genpar->getVertex().e();
00419 e.kn_p[0] = genpar->getMomentum().x();
00420 e.kn_p[1] = genpar->getMomentum().y();
00421 e.kn_p[2] = genpar->getMomentum().z();
00422 e.kn_p[3] = genpar->getMomentum().e();
00423 e.kn_pMag = genpar->getMomentum().vect().mag();
00424 e.kn_E = genpar->getMomentum().e();
00425 e.kn_pT = genpar->getMomentum().vect().perp();
00426 e.kn_eta = genpar->getMomentum().eta();
00427 e.kn_phi = genpar->getMomentum().phi();
00428 e.kn_pT = genpar->getMomentum().vect().perp();
00429 e.kn_eta = genpar->getMomentum().eta();
00430 e.kn_phi = genpar->getMomentum().phi();
00431 e.kn_pid = genpar->getPartId();
00432 e.kn_pidMo1 = genparMo1->getPartId();
00433 e.kn_lineMo1 = genparMo1->getLine();
00434 e.kn_pidMo2 = genparMo2->getPartId();
00435 e.kn_lineMo2 = genparMo2->getLine();
00436
00437 ++e.sc_isc;
00438
00439 e.sc_Erec = xsc->getEnergy();
00440 e.sc_ErecS1 = xsc->getSum1();
00441 e.sc_ErecS4 = xsc->getSum4();
00442 e.sc_ErecS9 = xsc->getSum9();
00443 e.sc_ErecS25 = xsc->getSum25();
00444 e.sc_ErecScaleFactor = xsc->getEnergyScaleFactor();
00445 e.sc_ErecUncorr = e.sc_Erec/e.sc_ErecScaleFactor;
00446 e.sc_eta = xsc->getPosition().eta();
00447 e.sc_phi = xsc->getPosition().phi();
00448
00449 e.sc_seedE = xsc->getSeedCluster()->getEnergy();
00450 e.sc_seedEta = xsc->getSeedCluster()->getPosition().eta();
00451 e.sc_seedPhi = xsc->getSeedCluster()->getPosition().phi();
00452 e.sc_nBrems = xsc->getNumberOfBrems();
00453
00454 TRefArray *cbrems = xsc->getBremClusters();
00455 XANACluster *seed = xsc->getSeedCluster();
00456 e.sc_seedE = seed->getEnergy();
00457 e.sc_seedEta = seed->getPosition().eta();
00458 e.sc_seedPhi = seed->getPosition().phi();
00459 e.sc_seedNxtals = seed->getNumberOfRecHits();
00460 TRefArray *seedHits = seed->getEmClusterHits();
00461 float threshold;
00462 if (fabs(e.sc_eta)<1.48) threshold = 2.*SIGMA_EB;
00463 else threshold = 2.*SIGMA_EE;
00464 for (int irh=0; irh<e.sc_seedNxtals; irh++) {
00465 XANACaloRecHit *caloRhit = (XANACaloRecHit*)(*seedHits)[irh];
00466 if (caloRhit->getEnergy()>threshold) e.sc_seedNxtals2sigma++;
00467 }
00468
00469 for (int ibc=0; ibc<e.sc_nBrems; ibc++) {
00470 XANAEmCluster *bc = (XANAEmCluster*)(*cbrems)[ibc];
00471 e.sc_bremEtot += bc->getEnergy();
00472 e.sc_bremNxtals += bc->getNumberOfRecHits();
00473 int bremNxtalsTmp = bc->getNumberOfRecHits();
00474 TRefArray *bremHits = bc->getEmClusterHits();
00475 for (int irh=0; irh<bremNxtalsTmp; irh++)
00476 if (((XANACaloRecHit*)(*bremHits)[irh])->getEnergy()>2*threshold) {
00477 e.sc_bremNxtals2sigma++;
00478 }
00479 }
00480
00481
00482
00483
00484 deltaR.clear();
00485 TObjArray *gsfTks = event->getElectronTracks();
00486 for (int ip=0; ip<event->getNumberOfElectronTracks(); ip++) {
00487 XANAElectronGSFTrack *gsfTk = (XANAElectronGSFTrack*)(*gsfTks)[ip];
00488 float eta_tk = gsfTk->getGSFMomentumAtVertexMode().eta();
00489 float eta_sc = xsc->getPosition().eta();
00490 float phi_tk = gsfTk->getGSFMomentumAtVertexMode().phi();
00491 float phi_sc = xsc->getPosition().phi();
00492
00493 float deta = eta_sc - eta_tk;
00494 float dphi = phi_sc - phi_tk;
00495 if (fabs(dphi) > PI) dphi = dphi < 0 ? 2.*PI + dphi : dphi - 2.*PI;
00496
00497 deltaR[ sqrt(deta*deta+dphi*dphi) ] = ip;
00498 }
00499
00500 map<float, int>::iterator it;
00501 vector<int> tkList;
00502
00503 for (it=deltaR.begin(); it!=deltaR.end(); it++) {
00504 if (it->first<1.) {
00505
00506 tkList.push_back(it->second);
00507 ++e.tk_ntk;
00508 }
00509 }
00510
00511 for (unsigned int iitk=0; iitk<tkList.size(); iitk++) {
00512 ++e.tk_itk;
00513 XANAElectronGSFTrack *gsfTk = (XANAElectronGSFTrack*)(*gsfTks)[tkList[iitk]];
00514
00515 e.tk_xetaFP = gsfTk->getPositionAtFirstPoint().eta();
00516 e.tk_xphiFP = gsfTk->getPositionAtFirstPoint().phi();
00517 e.tk_xFP[0] = gsfTk->getPositionAtFirstPoint().x();
00518 e.tk_xFP[1] = gsfTk->getPositionAtFirstPoint().y();
00519 e.tk_xFP[2] = gsfTk->getPositionAtFirstPoint().z();
00520
00521 e.tk_xetaLP = gsfTk->getPositionAtLastPoint().eta();
00522 e.tk_xphiLP = gsfTk->getPositionAtLastPoint().phi();
00523 e.tk_xLP[0] = gsfTk->getPositionAtLastPoint().x();
00524 e.tk_xLP[1] = gsfTk->getPositionAtLastPoint().y();
00525 e.tk_xLP[2] = gsfTk->getPositionAtLastPoint().z();
00526
00527 e.tk_petaFPmode = gsfTk->getGSFMomentumAtFirstPointMode().eta();
00528 e.tk_pphiFPmode = gsfTk->getGSFMomentumAtFirstPointMode().phi();
00529 e.tk_pFPmode[0] = gsfTk->getGSFMomentumAtFirstPointMode().x();
00530 e.tk_pFPmode[1] = gsfTk->getGSFMomentumAtFirstPointMode().y();
00531 e.tk_pFPmode[2] = gsfTk->getGSFMomentumAtFirstPointMode().z();
00532 e.tk_pFPmodeMag = gsfTk->getGSFMomentumAtFirstPointMode().mag();
00533 e.tk_petaFPmean = gsfTk->getGSFMomentumAtFirstPointWeightedMean().eta();
00534 e.tk_pphiFPmean = gsfTk->getGSFMomentumAtFirstPointWeightedMean().phi();
00535 e.tk_pFPmean[0] = gsfTk->getGSFMomentumAtFirstPointWeightedMean().x();
00536 e.tk_pFPmean[1] = gsfTk->getGSFMomentumAtFirstPointWeightedMean().y();
00537 e.tk_pFPmean[2] = gsfTk->getGSFMomentumAtFirstPointWeightedMean().z();
00538 e.tk_pFPmeanMag = gsfTk->getGSFMomentumAtFirstPointWeightedMean().mag();
00539
00540 e.tk_petaLPmode = gsfTk->getGSFMomentumAtLastPointMode().eta();
00541 e.tk_pphiLPmode = gsfTk->getGSFMomentumAtLastPointMode().phi();
00542 e.tk_pLPmode[0] = gsfTk->getGSFMomentumAtLastPointMode().x();
00543 e.tk_pLPmode[1] = gsfTk->getGSFMomentumAtLastPointMode().y();
00544 e.tk_pLPmode[2] = gsfTk->getGSFMomentumAtLastPointMode().z();
00545 e.tk_pLPmodeMag = gsfTk->getGSFMomentumAtLastPointMode().mag();
00546 e.tk_petaLPmean = gsfTk->getGSFMomentumAtLastPointWeightedMean().eta();
00547 e.tk_pphiLPmean = gsfTk->getGSFMomentumAtLastPointWeightedMean().phi();
00548 e.tk_pLPmean[0] = gsfTk->getGSFMomentumAtLastPointWeightedMean().x();
00549 e.tk_pLPmean[1] = gsfTk->getGSFMomentumAtLastPointWeightedMean().y();
00550 e.tk_pLPmean[2] = gsfTk->getGSFMomentumAtLastPointWeightedMean().z();
00551 e.tk_pLPmeanMag = gsfTk->getGSFMomentumAtLastPointWeightedMean().mag();
00552
00553 e.tk_petaVXmode = gsfTk->getGSFMomentumAtVertexMode().eta();
00554 e.tk_pphiVXmode = gsfTk->getGSFMomentumAtVertexMode().phi();
00555 e.tk_pVXmode[0] = gsfTk->getGSFMomentumAtVertexMode().x();
00556 e.tk_pVXmode[1] = gsfTk->getGSFMomentumAtVertexMode().y();
00557 e.tk_pVXmode[2] = gsfTk->getGSFMomentumAtVertexMode().z();
00558 e.tk_pVXmodeMag = gsfTk->getGSFMomentumAtVertexMode().mag();
00559 e.tk_petaVXmean = gsfTk->getGSFMomentumAtVertexWeightedMean().eta();
00560 e.tk_pphiVXmean = gsfTk->getGSFMomentumAtVertexWeightedMean().phi();
00561 e.tk_pVXmean[0] = gsfTk->getGSFMomentumAtVertexWeightedMean().x();
00562 e.tk_pVXmean[1] = gsfTk->getGSFMomentumAtVertexWeightedMean().y();
00563 e.tk_pVXmean[2] = gsfTk->getGSFMomentumAtVertexWeightedMean().z();
00564 e.tk_pVXmeanMag = gsfTk->getGSFMomentumAtVertexWeightedMean().mag();
00565
00566 e.tk_chi2OverDof = gsfTk->getChi2OverDof();
00567 e.tk_chi2OverDofBack = gsfTk->getBackwardFitChi2OverDof();
00568 e.tk_tkNhits = gsfTk->getNumberOfHits();
00569 e.tk_tkNLostHits = gsfTk->getNumberOfLostHits();
00570 e.tk_chargeRec = gsfTk->getCharge();
00571 e.tk_highestWeightFP = gsfTk->getHighestWeightAtFirstPoint();
00572 e.tk_highestWeightLP = gsfTk->getHighestWeightAtLastPoint();
00573 e.tk_highestWeightVX = gsfTk->getHighestWeightAtVertex();
00574
00575
00576 XANAVertex *vtx = gsfTk->getVertex();
00577 if (vtx!=NULL) {
00578 e.vx_x[0] = vtx->getPosition().x();
00579 e.vx_x[1] = vtx->getPosition().y();
00580 e.vx_x[2] = vtx->getPosition().z();
00581 e.vx_xeDiag[0] = vtx->getPositionErrorDiag().x();
00582 e.vx_xeDiag[1] = vtx->getPositionErrorDiag().y();
00583 e.vx_xeDiag[2] = vtx->getPositionErrorDiag().z();
00584 e.vx_xeOffDiag[0] = vtx->getPositionErrorOffDiag().x();
00585 e.vx_xeOffDiag[1] = vtx->getPositionErrorOffDiag().y();
00586 e.vx_xeOffDiag[2] = vtx->getPositionErrorOffDiag().z();
00587 e.vx_chi2 = vtx->getChi2();
00588 e.vx_chi2Norm = vtx->getNormalizedChi2();
00589 e.vx_Ndof = (int)vtx->getNdof();
00590 e.vx_vectSum[0] = vtx->getVectorialSum().x();
00591 e.vx_vectSum[1] = vtx->getVectorialSum().y();
00592 e.vx_vectSum[2] = vtx->getVectorialSum().z();
00593
00594
00595 }
00596
00597
00598 e.el_iel=0;
00599 XANAElectronCandidate *elCan1 = xsc->getElectronCandidate();
00600 XANAElectronCandidate *elCan = gsfTk->getElectronCandidate();
00601 if (elCan1==elCan && elCan!=NULL) {
00602 e.el_iel = 1;
00603 e.el_phi = elCan->getPhi();
00604 e.el_eta = elCan->getEta();
00605 e.el_pVX[0] = elCan->getMomentumAtVertex().x();
00606 e.el_pVX[1] = elCan->getMomentumAtVertex().y();
00607 e.el_pVX[2] = elCan->getMomentumAtVertex().z();
00608 e.el_pVXMag = elCan->getMomentumAtVertex().mag();
00609 e.el_pTVX = elCan->getPt();
00610
00611 e.el_DetaScTkVx = elCan->getDeltaEtaSuperClusterAtVtx();
00612 e.el_DetaScTkEc = elCan->getDeltaEtaSuperClusterAtCalo();
00613 e.el_DetaSeedTkEc = elCan->getDeltaEtaSeedClusterAtCalo();
00614 e.el_DphiScTkVx = elCan->getDeltaEtaSuperClusterAtVtx();
00615 e.el_DphiScTkEc = elCan->getDeltaEtaSuperClusterAtCalo();
00616 e.el_DphiSeedTkEc = elCan->getDeltaEtaSeedClusterAtCalo();
00617
00618 e.el_tkIso = elCan->getTrackIsolation();
00619 e.el_ecIso = elCan->getCaloIsolation();
00620 e.el_EOP = elCan->getESuperClusterOverP();
00621 e.el_EseedOP = elCan->getESeedClusterOverP();
00622 e.el_HOE = elCan->getHadronicOverEm();
00623 }
00624
00625 tree->Fill();
00626 ifilled=1;
00627 k.ifirst=1;
00628
00629 }
00630 if (!tkList.size()) {
00631 tree->Fill();
00632 ifilled=1;
00633 k.ifirst=1;
00634
00635 }
00636 }
00637 if (!e.sc_nsc) {
00638 tree->Fill();
00639
00640 }
00641 event->clear();
00642 geantEvent->clear();
00643 generatorEvent->clear();
00644 }
00645 F->Close();
00646 delete F;
00647 }
00648
00649 cout << endl;
00650 cout << e.ievent << " events and " << nb << " bytes read " << endl;
00651 cout << endl;
00652
00653 fout->Write();
00654 fout->Close();
00655 return 0;
00656 }
00657
00658
00659
00660
00661
00662 void kineReset(simpleGeneratorEvent *k) {
00663 k->ifirst=k->nple=0;
00664 for (int i=0; i<10; i++) {
00665 for (int j=0; j<4; j++) {
00666 k->x[i][j]=0;
00667 k->p[i][j]=0;
00668 }
00669 k->pid[i]=k->mom[i]=0;
00670 k->pMag[i]=k->pT[i]=k->eta[i]=k->phi[i]=0;
00671 }
00672 k->mH=k->mZ1=k->mZ2=0;
00673 }
00674
00675
00676 void globalReset(simpleEvent *e)
00677 {
00678 e->kn_pid=0;
00679 e->kn_E=e->kn_pT=e->kn_eta=e->kn_phi=0;
00680 e->kn_pMag=0;
00681
00682 e->sc_Erec=e->sc_ErecUncorr=e->sc_ErecScaleFactor=0;
00683 e->sc_ErecS1=e->sc_ErecS4=e->sc_ErecS9=e->sc_ErecS25=0;
00684 e->sc_eta=e->sc_phi=0;
00685 e->sc_seedE=e->sc_seedEta=e->sc_seedPhi=0;
00686 e->sc_seedNxtals=e->sc_seedNxtals2sigma=0;
00687 e->sc_nBrems=0;
00688 e->sc_bremEtot=0;
00689 e->sc_bremNxtals=e->sc_bremNxtals2sigma=0;
00690 e->tk_ntk=e->tk_itk=0;
00691 e->tk_xetaFP=e->tk_xphiFP=0;
00692 e->tk_xetaLP=e->tk_xphiLP=0;
00693 e->tk_petaFPmode=e->tk_pphiFPmode=0;
00694 e->tk_pFPmodeMag=0;
00695 e->tk_petaFPmean=e->tk_pphiFPmean=0;
00696 e->tk_pFPmeanMag=0;
00697 e->tk_petaLPmode=e->tk_pphiLPmode=0;
00698 e->tk_pLPmodeMag=0;
00699 e->tk_petaLPmean=e->tk_pphiLPmean=0;
00700 e->tk_pLPmeanMag=0;
00701 e->tk_petaVXmode=e->tk_pphiVXmode=0;
00702 e->tk_pVXmodeMag=0;
00703 e->tk_petaVXmean=e->tk_pphiVXmean=0;
00704 e->tk_pVXmeanMag=0;
00705 e->tk_tkNhits=e->tk_tkNLostHits=0;
00706 e->tk_chi2OverDof=e->tk_chi2OverDofBack=0;
00707 e->tk_chargeRec=0;
00708 e->tk_highestWeightFP=e->tk_highestWeightLP=e->tk_highestWeightVX=0;
00709 e->vx_chi2=0;
00710 e->vx_chi2Norm=0;
00711 e->vx_Ndof=0;
00712 e->el_phi=0;
00713 e->el_eta=0;
00714 e->el_pVXMag=0;
00715 e->el_pTVX=0;
00716 e->el_pECMag=0;
00717 e->el_DetaScTkVx=0;
00718 e->el_DetaScTkEc=0;
00719 e->el_DetaSeedTkEc=0;
00720 e->el_DphiScTkVx=0;
00721 e->el_DphiScTkEc=0;
00722 e->el_DphiSeedTkEc=0;
00723 e->el_tkIso=0;
00724 e->el_ecIso=0;
00725 e->el_EOP=0;
00726 e->el_EseedOP=0;
00727 e->el_HOE=0;
00728
00729 e->kn_deltaR=0;
00730 for (int i=0; i<10; i++) {
00731 if (i<4) {
00732 e->kn_x[i]=0;
00733 e->kn_p[i]=0;
00734 }
00735 if (i<3) {
00736 e->tk_xFP[i]=0;
00737 e->tk_xLP[i]=0;
00738 e->tk_pFPmode[i]=0;
00739 e->tk_pFPmean[i]=0;
00740 e->tk_pLPmode[i]=0;
00741 e->tk_pLPmean[i]=0;
00742 e->tk_pVXmode[i]=0;
00743 e->tk_pVXmean[i]=0;
00744 e->el_pVX[i]=0;
00745 e->el_pEC[i]=0;
00746 e->vx_x[i]=0;
00747 e->vx_xeDiag[i]=0;
00748 e->vx_xeOffDiag[i]=0;
00749 e->vx_vectSum[i]=0;
00750 }
00751 }
00752 }
00753
00754
00755 int mom(int idx, TObjArray *tracks)
00756 {
00757 int pid=0;
00758 while (idx>0) {
00759 XANAGeantTrack *gtmp = (XANAGeantTrack*)(*tracks)[idx-1];
00760 idx = gtmp->getMotherIndex();
00761 pid = gtmp->getPartId();
00762 }
00763 return pid;
00764 }
00765
00766
00767
00768 XANAGeantVertex *getDaughterVertex(XANAGeantTrack *gtk, TObjArray *tracks)
00769 {
00770 for (int i=0; i<tracks->GetEntriesFast(); i++) {
00771 if (((XANAGeantTrack*)(*tracks)[i])->getMotherIndex()==gtk->getIndex()
00772 && ((XANAGeantTrack*)(*tracks)[i])->getVertex()->getNumberOfTracks()==2) {
00773 return ((XANAGeantTrack*)(*tracks)[i])->getVertex();
00774 }
00775 }
00776 return NULL;
00777 }
