Fork me on GitHub

Changes in / [54ec182:4564cba] in git


Ignore:
Files:
11 added
17 edited

Legend:

Unmodified
Added
Removed

  • Makefile

    r54ec182 r4564cba  
    464464        modules/ParticleDensity.h \
    465465        modules/TruthVertexFinder.h \
    466         modules/ExampleModule.h
     466        modules/ExampleModule.h \
     467        modules/LLPFilter.h \
     468        modules/CscClusterEfficiency.h \
     469        modules/CscClusterId.h
    467470tmp/modules/ModulesDict$(PcmSuf): \
    468471        tmp/modules/ModulesDict.$(SrcSuf)
    … …  
    518521        classes/DelphesFactory.h \
    519522        classes/SortableObject.h
     523tmp/classes/DelphesCscClusterFormula.$(ObjSuf): \
     524        classes/DelphesCscClusterFormula.$(SrcSuf) \
     525        classes/DelphesCscClusterFormula.h \
     526        classes/DelphesClasses.h
    520527tmp/classes/DelphesCylindricalFormula.$(ObjSuf): \
    521528        classes/DelphesCylindricalFormula.$(SrcSuf) \
    … …  
    742749        external/ExRootAnalysis/ExRootFilter.h \
    743750        external/ExRootAnalysis/ExRootResult.h
     751tmp/modules/CscClusterEfficiency.$(ObjSuf): \
     752        modules/CscClusterEfficiency.$(SrcSuf) \
     753        modules/CscClusterEfficiency.h \
     754        classes/DelphesClasses.h \
     755        classes/DelphesFactory.h \
     756        classes/DelphesCscClusterFormula.h \
     757        external/ExRootAnalysis/ExRootClassifier.h \
     758        external/ExRootAnalysis/ExRootFilter.h \
     759        external/ExRootAnalysis/ExRootResult.h
     760tmp/modules/CscClusterId.$(ObjSuf): \
     761        modules/CscClusterId.$(SrcSuf) \
     762        modules/CscClusterId.h \
     763        classes/DelphesClasses.h \
     764        classes/DelphesFactory.h \
     765        classes/DelphesCscClusterFormula.h \
     766        external/ExRootAnalysis/ExRootClassifier.h \
     767        external/ExRootAnalysis/ExRootFilter.h \
     768        external/ExRootAnalysis/ExRootResult.h
    744769tmp/modules/DecayFilter.$(ObjSuf): \
    745770        modules/DecayFilter.$(SrcSuf) \
    … …  
    882907        external/ExRootAnalysis/ExRootFilter.h \
    883908        external/ExRootAnalysis/ExRootResult.h
     909tmp/modules/LLPFilter.$(ObjSuf): \
     910        modules/LLPFilter.$(SrcSuf) \
     911        modules/LLPFilter.h \
     912        classes/DelphesClasses.h \
     913        classes/DelphesFactory.h \
     914        classes/DelphesFormula.h \
     915        external/ExRootAnalysis/ExRootClassifier.h \
     916        external/ExRootAnalysis/ExRootFilter.h \
     917        external/ExRootAnalysis/ExRootResult.h
    884918tmp/modules/LeptonDressing.$(ObjSuf): \
    885919        modules/LeptonDressing.$(SrcSuf) \
    … …  
    11631197DELPHES_OBJ +=  \
    11641198        tmp/classes/DelphesClasses.$(ObjSuf) \
     1199        tmp/classes/DelphesCscClusterFormula.$(ObjSuf) \
    11651200        tmp/classes/DelphesCylindricalFormula.$(ObjSuf) \
    11661201        tmp/classes/DelphesFactory.$(ObjSuf) \
    … …  
    12271262        tmp/modules/ClusterCounting.$(ObjSuf) \
    12281263        tmp/modules/ConstituentFilter.$(ObjSuf) \
     1264        tmp/modules/CscClusterEfficiency.$(ObjSuf) \
     1265        tmp/modules/CscClusterId.$(ObjSuf) \
    12291266        tmp/modules/DecayFilter.$(ObjSuf) \
    12301267        tmp/modules/Delphes.$(ObjSuf) \
    … …  
    12421279        tmp/modules/JetFlavorAssociation.$(ObjSuf) \
    12431280        tmp/modules/JetPileUpSubtractor.$(ObjSuf) \
     1281        tmp/modules/LLPFilter.$(ObjSuf) \
    12441282        tmp/modules/LeptonDressing.$(ObjSuf) \
    12451283        tmp/modules/Merger.$(ObjSuf) \
    … …  
    19281966        @touch $@
    19291967
     1968modules/LLPFilter.h: \
     1969        classes/DelphesModule.h
     1970        @touch $@
     1971
    19301972external/fastjet/internal/MinHeap.hh: \
    19311973        external/fastjet/internal/base.hh
    … …  
    22572299        @touch $@
    22582300
     2301modules/CscClusterEfficiency.h: \
     2302        classes/DelphesModule.h
     2303        @touch $@
     2304
    22592305external/fastjet/PseudoJetStructureBase.hh: \
    22602306        external/fastjet/internal/base.hh
    … …  
    23812427external/fastjet/config.h: \
    23822428        external/fastjet/config_win.h
     2429        @touch $@
     2430
     2431modules/CscClusterId.h: \
     2432        classes/DelphesModule.h
    23832433        @touch $@
    23842434

  • classes/ClassesLinkDef.h

    r54ec182 r4564cba  
    5757#pragma link C++ class Electron+;
    5858#pragma link C++ class Muon+;
     59#pragma link C++ class CscCluster+;
     60
    5961#pragma link C++ class Jet+;
    6062#pragma link C++ class Track+;
    … …  
    6668
    6769#endif
    68 

  • classes/DelphesClasses.cc

    r54ec182 r4564cba  
    3535CompBase *Electron::fgCompare = CompPT<Electron>::Instance();
    3636CompBase *Muon::fgCompare = CompPT<Muon>::Instance();
     37
    3738CompBase *Jet::fgCompare = CompPT<Jet>::Instance();
    3839CompBase *Track::fgCompare = CompPT<Track>::Instance();
    … …  
    4243CompBase *Vertex::fgCompare = CompSumPT2<Vertex>::Instance();
    4344CompBase *Candidate::fgCompare = CompMomentumPt<Candidate>::Instance();
     45CompBase *CscCluster::fgCompare = CompE<CscCluster>::Instance();
    4446
    4547//------------------------------------------------------------------------------
    … …  
    218220  Position(0.0, 0.0, 0.0, 0.0),
    219221  InitialPosition(0.0, 0.0, 0.0, 0.0),
     222  DecayPosition(0.0, 0.0, 0.0, 0.0),
    220223  PositionError(0.0, 0.0, 0.0, 0.0),
    221224  Area(0.0, 0.0, 0.0, 0.0),
    … …  
    230233  Phi(0), ErrorPhi(0),
    231234  Xd(0), Yd(0), Zd(0),
     235  XFirstHit(0), YFirstHit(0), ZFirstHit(0),
    232236  Nclusters(0.0),
    233237  dNdx(0.0),
    … …  
    388392  object.Position = Position;
    389393  object.InitialPosition = InitialPosition;
     394  object.DecayPosition = DecayPosition;
    390395  object.PositionError = PositionError;
    391396  object.Area = Area;
    … …  
    409414  object.Yd = Yd;
    410415  object.Zd = Zd;
     416  object.XFirstHit = XFirstHit;
     417  object.YFirstHit = YFirstHit;
     418  object.ZFirstHit = ZFirstHit;
    411419  object.Nclusters = Nclusters;
    412420  object.dNdx = dNdx;
    … …  
    523531  Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
    524532  InitialPosition.SetXYZT(0.0, 0.0, 0.0, 0.0);
     533  DecayPosition.SetXYZT(0.0, 0.0, 0.0, 0.0);
    525534  Area.SetXYZT(0.0, 0.0, 0.0, 0.0);
    526535  TrackCovariance.Zero();
    … …  
    544553  Yd = 0.0;
    545554  Zd = 0.0;
     555  XFirstHit = 0.0;
     556  YFirstHit = 0.0;
     557  ZFirstHit = 0.0;
    546558  Nclusters = 0.0;
    547559  dNdx = 0.0;

  • classes/DelphesClasses.h

    r54ec182 r4564cba  
    158158  Float_t Z; // particle vertex position (z component) | hepevt.vhep[number][2]
    159159
     160  Float_t decayX;
     161  Float_t decayY;
     162  Float_t decayZ;
     163  Float_t decayT;
     164
    160165  static CompBase *fgCompare; //!
    161166  const CompBase *GetCompare() const { return fgCompare; }
    … …  
    457462  Float_t Yd; // Y coordinate of point of closest approach to vertex
    458463  Float_t Zd; // Z coordinate of point of closest approach to vertex
     464
     465  Float_t XFirstHit; // X coordinate of point of closest approach to vertex
     466  Float_t YFirstHit; // Y coordinate of point of closest approach to vertex
     467  Float_t ZFirstHit; // Z coordinate of point of closest approach to vertex
    459468
    460469  Float_t L; // track path length
    … …  
    565574  Float_t Zd; // Z coordinate of point of closest approach to vertex
    566575
     576  Float_t XFirstHit; // X coordinate of point of closest approach to vertex
     577  Float_t YFirstHit; // Y coordinate of point of closest approach to vertex
     578  Float_t ZFirstHit; // Z coordinate of point of closest approach to vertex
     579
    567580  Float_t L; // track path length
    568581  Float_t D0; // track transverse impact parameter
    … …  
    638651  ClassDef(HectorHit, 1)
    639652};
     653//---------------------------------------------------------------------------
     654
     655class CscCluster: public SortableObject
     656{
     657public:
     658  Float_t Eta; // eta of LLP
     659  Float_t Phi; // phi of LLP
     660  Float_t PT; // pt of LLP
     661  Float_t Px;// px of LLP
     662  Float_t Py;// py of LLP
     663  Float_t Pz;// pz of LLP
     664  Float_t E; // E of LLP
     665  Float_t Ehad; // had energy of LLP
     666  Float_t Eem; // em energy of LLP
     667  Float_t pid; // LLP pid
     668  Float_t T; // LLP decay time-photon travel time
     669  Float_t X; // LLP decay x
     670  Float_t Y; //  LLP decay y
     671  Float_t Z; //  LLP decay z
     672  Float_t R; //  LLP decay z
     673  Float_t beta; // LLP beta
     674  Float_t ctau; //LLP ctau
     675
     676
     677  static CompBase *fgCompare; //!
     678  const CompBase *GetCompare() const { return fgCompare; }
     679
     680  ClassDef(CscCluster, 5)
     681};
    640682
    641683//---------------------------------------------------------------------------
    … …  
    682724  Float_t DeltaPhi;
    683725
    684   TLorentzVector Momentum, Position, InitialPosition, PositionError, Area;
     726  TLorentzVector Momentum, Position, InitialPosition, PositionError, DecayPosition, Area;
    685727
    686728  Float_t L; // path length
    … …  
    708750  Float_t Zd;
    709751
     752  Float_t XFirstHit;
     753  Float_t YFirstHit;
     754  Float_t ZFirstHit;
     755
    710756  // tracking resolution
    711757

  • classes/DelphesHepMC2Reader.cc

    r54ec182 r4564cba  
    438438{
    439439  Candidate *candidate;
     440  Candidate *candidateDaughter;
    440441  map<int, pair<int, int> >::iterator itMotherMap;
    441442  map<int, pair<int, int> >::iterator itDaughterMap;
    … …  
    446447    candidate = static_cast<Candidate *>(allParticleOutputArray->At(i));
    447448
     449
    448450    if(candidate->M1 > 0)
    449451    {
    … …  
    477479        candidate->D1 = -1;
    478480        candidate->D2 = -1;
    479       }
     481        const TLorentzVector &decayPosition = candidate->Position;
     482        candidate->DecayPosition.SetXYZT(decayPosition.X(), decayPosition.Y(), decayPosition.Z(), decayPosition.T());// decay position
     483     }
    480484      else
    481485      {
    482486        candidate->D1 = itDaughterMap->second.first;
    483487        candidate->D2 = itDaughterMap->second.second;
    484       }
    485     }
    486   }
    487 }
    488 
    489 //---------------------------------------------------------------------------
     488        candidateDaughter = static_cast<Candidate *>(allParticleOutputArray->At(candidate->D1));
     489        const TLorentzVector &decayPosition = candidateDaughter->Position;
     490        candidate->DecayPosition.SetXYZT(decayPosition.X(), decayPosition.Y(), decayPosition.Z(), decayPosition.T());// decay position
     491       
     492
     493      }
     494    }
     495  }
     496}
     497
     498//---------------------------------------------------------------------------

  • classes/DelphesHepMC3Reader.cc

    r54ec182 r4564cba  
    458458  TObjArray *array;
    459459  Candidate *candidate;
     460  Candidate *candidateDaughter;
    460461  TParticlePDG *pdgParticle;
    461462  int pdgCode;
    … …  
    573574        candidate->D1 = -1;
    574575        candidate->D2 = -1;
     576        const TLorentzVector &decayPosition = candidate->Position;
     577        candidate->DecayPosition.SetXYZT(decayPosition.X(), decayPosition.Y(), decayPosition.Z(), decayPosition.T());// decay position
    575578      }
    576579      else
    … …  
    578581        candidate->D1 = itDaughterMap->second.first;
    579582        candidate->D2 = itDaughterMap->second.second;
    580       }
    581     }
    582   }
    583 }
    584 
    585 //---------------------------------------------------------------------------
     583        candidateDaughter = static_cast<Candidate *>(allParticleOutputArray->At(candidate->D1));
     584        const TLorentzVector &decayPosition = candidateDaughter->Position;
     585        candidate->DecayPosition.SetXYZT(decayPosition.X(), decayPosition.Y(), decayPosition.Z(), decayPosition.T());// decay position
     586      }
     587    }
     588  }
     589}
     590
     591//---------------------------------------------------------------------------

  • examples/ExamplePVtxFind.C

    r54ec182 r4564cba  
    4545        TH1D* hTrFoundB = new TH1D("hTrFoundB", "Found BAD  primary tracks", 41, -0.5, 40.5);
    4646        TH1D* hTrDiff = new TH1D("hTrDiff", "Found - available tracks", 21, -10.5, 10.5);
     47        // Chi2 for cuts
     48        TH1D* hChi2SngP = new TH1D("hChi2SngP", "#chi^2 for single tracks", 200, 0., 50.);
     49        TH1D* hChi2MaxP = new TH1D("hChi2MaxP", "#chi^2 max contribution", 200, 0., 50.);
     50        TH1D* hChi2SngN = new TH1D("hChi2SngN", "#chi^2 for single tracks", 200, 0., 50.);
     51        TH1D* hChi2MaxN = new TH1D("hChi2MaxN", "#chi^2 max contribution", 200, 0., 50.);
     52       
     53       
    4754        //
    4855        // Loop over all events
    … …  
    5360                treeReader->ReadEntry(entry);
    5461                Int_t NtrG = branchTrack->GetEntries();
    55                 if(entry%500 ==0)std::cout << "Event "<<entry<<" opened containing " << NtrG << " tracks" << std::endl;
    56                 TVectorD** pr = new TVectorD * [NtrG];
    57                 TMatrixDSym** cv = new TMatrixDSym * [NtrG];
     62                TVectorD** pr = new TVectorD * [NtrG];          // Track Parameters
     63                TMatrixDSym** cv = new TMatrixDSym * [NtrG];    // Track covariances
     64                Bool_t *fprvx = new Bool_t[NtrG];               // Primary vertex flag
    5865                //
    5966                // test Particle branch
    … …  
    8895                                pr[it] = new TVectorD(obsPar);
    8996                                cv[it] = new TMatrixDSym(trk->CovarianceMatrix());
    90                                 // Check covariance
    91                                 /*
    92                                 TMatrixDSymEigen Eign(*cv[it]);
    93                                 TVectorD lambda = Eign.GetEigenValues();
    94                                 Bool_t good = kTRUE;
    95                                 for(Int_t i=0; i<5; i++) if(lambda(i)<0.0) good = kFALSE;
    96                                 if(!good){
    97                                         std::cout<<"Cov["<<it<<"]"<<" eigenvalues: "
    98                                         <<lambda(0)<<", "<<lambda(1)<<", "<<lambda(2)<<", "
    99                                         <<lambda(3)<<", "<<lambda(4)<<std::endl;
    100                                         TMatrixDSym Ctmp = *cv[it];
    101                                         std::cout<<"Error D (dir - cov) "<<trk->ErrorD0
    102                                         <<" - "<<TMath::Sqrt(Ctmp(0,0))<<std::endl;
    103                                 }
    104                                 */
    105                                 //std::cout << "Track loaded it= " << it << std::endl;
    10697                        //
    10798                        // Find true primary vertex
    … …  
    114105                                Double_t z = gp->Z;
    115106                                Bool_t prim = kTRUE;    // Is primary?
     107                                fprvx[it] = kFALSE;
    116108                                Int_t mp = gp->M1;      // Mother
    117109                                while (mp > 0) {
    … …  
    129121                                if (prim) {             // It's a primary track
    130122                                        Nprim++;
     123                                        fprvx[it] = kTRUE;
    131124                                        xpv = x;        // Store true primary
    132125                                        ypv = y;
    … …  
    136129                        }               // End loop on tracks
    137130                }
    138                 std::cout<<"PVtxFind true vertex: Nprim= "<<Nprim<<", x,y,z= "<<xpv<<", "<<ypv<<", "<<zpv<<std::endl;
     131                if(entry%500 ==0){
     132                  std::cout << "Event "<<entry<<" opened containing " << NtrG << " / "<< Nprim
     133                  << "   Total / primary tracks"<< std::endl;
     134                }
     135                //std::cout<<"PVtxFind true vertex: Nprim= "<<Nprim<<", x,y,z= "<<xpv<<", "<<ypv<<", "<<zpv<<std::endl;
    139136                //
    140137                // Find primary vertex
    … …  
    161158                TVectorD** PrSk = new TVectorD * [1];
    162159                TMatrixDSym** CvSk = new TMatrixDSym * [1];
    163                 Double_t MaxChi2 = 2.0+3*2.;
     160                Double_t MaxChi2 = 9.;
    164161                for (Int_t n = 0; n < NtrG; n++) {
    165162                        PrSk[0] = new TVectorD(*pr[n]);
    … …  
    169166                        Double_t Chi2One = Vskim->GetVtxChi2();
    170167                        //std::cout<<"Track "<<n<<", Chi2 = "<<Chi2One<<std::endl;
     168                        if(fprvx[n])hChi2SngP->Fill(Chi2One);
     169                        else hChi2SngN->Fill(Chi2One);
     170                        //
    171171                        if (Chi2One < MaxChi2) {
    172172                                nSkimmed[nSkim] = n;
    … …  
    184184                // Load tracks for primary fit
    185185                Int_t MinTrk = 1;       // Minumum # tracks for vertex fit
     186                std::vector<Int_t> trnum;
    186187                if (nSkim >= MinTrk) {
    187188                        TVectorD** PrFit = new TVectorD * [nSkim];
    … …  
    190191                                PrFit[n] = new TVectorD(*pr[nSkimmed[n]]);
    191192                                CvFit[n] = new TMatrixDSym(*cv[nSkimmed[n]]);
     193                                trnum.push_back(nSkimmed[n]);
    192194                        }
    193195                        delete[] nSkimmed;
    … …  
    200202                        //
    201203                        // Remove tracks with large chi2
    202                         Double_t MaxChi2Fit = 4.5;
     204                        Double_t MaxChi2Fit = 8.0;
    203205                        Bool_t Done = kFALSE;
    204206                        while (!Done) {
    … …  
    213215                                Double_t Chi2Mx = Chi2L[iMax];
    214216                                //std::cout << "Chi2Mx "<<Chi2Mx << std::endl;
     217                                if(fprvx[trnum[iMax]])hChi2MaxP->Fill(Chi2Mx);
     218                                else hChi2MaxN->Fill(Chi2Mx);
    215219                                if (Chi2Mx > MaxChi2Fit && Nfound > 1) {
    216220                                        //std::cout << "Before remove.  Nfound = "<<Nfound << std::endl;
    217221                                        Vtx->RemoveTrk(iMax);
     222                                        trnum.erase(trnum.begin() + iMax);
    218223                                        //std::cout << "After remove." << std::endl;
    219224                                        Nfound--;
    … …  
    230235                        if (Nfound >= Nmin) {
    231236                                TVectorD xvtx = Vtx->GetVtx();
    232                                 std::cout << "Found vertex " << xvtx(0)<<", "<<xvtx(1)<<", "<<xvtx(2) << std::endl;
     237                                //std::cout << "Found vertex " << xvtx(0)<<", "<<xvtx(1)<<", "<<xvtx(2) << std::endl;
    233238                                TMatrixDSym covX = Vtx->GetVtxCov();
    234239                                Double_t Chi2 = Vtx->GetVtxChi2();
    … …  
    256261                        delete[] PrFit;
    257262                        delete[] CvFit;
     263                        delete[] fprvx;
     264                        trnum.clear();
    258265                }
    259266
    … …  
    271278        TCanvas* Cnv = new TCanvas("Cnv", "Delphes primary vertex pulls", 50, 50, 900, 500);
    272279        Cnv->Divide(2, 2);
    273         Cnv->cd(1); gPad->SetLogy(1); gStyle->SetOptFit(1111);
     280        Cnv->cd(1); gPad->SetLogy(1); gStyle->SetOptFit(1111); gStyle->SetOptStat(111111);
    274281        hXpull->Fit("gaus"); hXpull->Draw();
    275         Cnv->cd(2); gPad->SetLogy(1); gStyle->SetOptFit(1111);
     282        Cnv->cd(2); gPad->SetLogy(1); gStyle->SetOptFit(1111); gStyle->SetOptStat(111111);
    276283        hYpull->Fit("gaus"); hYpull->Draw();
    277         Cnv->cd(3); gPad->SetLogy(1); gStyle->SetOptFit(1111);
     284        Cnv->cd(3); gPad->SetLogy(1); gStyle->SetOptFit(1111); gStyle->SetOptStat(111111);
    278285        hZpull->Fit("gaus"); hZpull->Draw();
    279286        Cnv->cd(4); hChi2->Draw();
    … …  
    288295        CnvN->cd(3);
    289296        hTrDiff->Draw();
     297        //
     298        TCanvas* CnvCh = new TCanvas("CnvCh", "#chi^2", 200, 200, 900, 500);
     299        CnvCh->Divide(2,1);
     300        CnvCh->cd(1);
     301        hChi2SngP->Draw();
     302        hChi2SngN->SetLineColor(kRed);
     303        hChi2SngN->Draw("SAME");
     304        CnvCh->cd(2);
     305        hChi2MaxP->Draw();
     306        hChi2MaxN->SetLineColor(kRed);
     307        hChi2MaxN->Draw("SAME");
    290308}

  • external/TrackCovariance/ObsTrk.cc

    r54ec182 r4564cba  
    7474        fCovILC.ResizeTo(5, 5);
    7575        fGenPar = XPtoPar(fGenX, fGenP, Q);
     76        fGenParMm = ParToMm(fGenPar);
    7677        fGenParACTS = ParToACTS(fGenPar);
    7778        fGenParILC = ParToILC(fGenPar);
    7879        //
    7980        fObsPar = GenToObsPar(fGenPar);
     81        fObsParMm = ParToMm(fObsPar);
    8082        fObsParACTS = ParToACTS(fObsPar);
    8183        fObsParILC = ParToILC(fObsPar);
    … …  
    8385        fObsP = ParToP(fObsPar);
    8486        fObsQ = ParToQ(fObsPar);
     87        fCovMm = CovToMm(fCov);
    8588        fCovACTS = CovToACTS(fObsPar, fCov);
    8689        fCovILC = CovToILC(fCov);
    … …  
    117120        //if (pt > maxPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is above grid range of " << maxPt << std::endl;
    118121        Double_t minAn = fGC->GetMinAng();
    119         //if (angd < minAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd 
     122        //if (angd < minAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd
    120123        //      << " is below grid range of " << minAn << std::endl;
    121124        Double_t maxAn = fGC->GetMaxAng();
    … …  
    131134        Double_t ZinNeg = fG->GetZminNeg();
    132135        Bool_t inside = TrkUtil::IsInside(fGenX, Rin, ZinNeg, ZinPos); // Check if in inner box
     136        SolTrack* trk = new SolTrack(fGenX, fGenP, fG);
     137        Double_t Xfirst, Yfirst, Zfirst;
     138        Int_t iLay = trk->FirstHit(Xfirst, Yfirst, Zfirst);
     139        fXfirst = TVector3(Xfirst, Yfirst, Zfirst);
     140  //std::cout<<"obs trk: "<<Xfirst<<","<<Yfirst<<","<<Zfirst<<std::endl;
     141
    133142        if (inside)
    134143        {
    … …  
    144153                //std::cout<<"ObsTrk:: outside: x= "<<fGenX(0)<<", y= "<<fGenX(1)
    145154                //                         <<", z= "<<fGenX(2)<<std::endl;
    146                 SolTrack* trk = new SolTrack(fGenX, fGenP, fG);
    147155                Bool_t Res = kTRUE; Bool_t MS = kTRUE;
    148156                trk->CovCalc(Res, MS);                                  // Calculate covariance matrix
    149                 Cov = trk->Cov();                                       // Track covariance
    150                 delete trk;
    151         }
     157                Cov = trk->Cov();
     158        }                                       // Track covariance
     159        delete trk;
    152160        //
    153161        fCov = Cov;

  • external/TrackCovariance/ObsTrk.h

    r54ec182 r4564cba  
    4646        TMatrixDSym fCovILC;                    // Covariance of track parameters in ILC format
    4747                                                                        // (d0, phi0, w, z0, tan(lambda))
     48        TVector3 fXfirst;                       // x,y,z of first track hit
    4849        //
    4950        // Service routines
    … …  
    8990        TMatrixDSym GetCovMm()  { return fCov; }        // in mm
    9091        TMatrixDSym GetCovACTS(){ return fCovACTS; }
    91         TMatrixDSym GetCovILC(){ return fCovILC; }
    92         //
     92        TMatrixDSym GetCovILC() { return fCovILC; }
     93        // First hit
     94        TVector3 GetFirstHit()  { return fXfirst; }
    9395};
    9496

  • external/TrackCovariance/SolTrack.cc

    r54ec182 r4564cba  
    176176{
    177177        //
    178         // Return lists of hits associated to a track including all scattering layers. 
     178        // Return lists of hits associated to a track including all scattering layers.
    179179        // Return value is the total number of measurement hits
    180180        // kmh = total number of measurement layers hit for given type
    … …  
    190190        {
    191191                Double_t R; Double_t phi; Double_t zz;
    192                 if (HitLayer(i, R, phi, zz)) 
     192                if (HitLayer(i, R, phi, zz))
    193193                {
    194194                        zhh[kh] = zz;
    … …  
    207207{
    208208        //
    209         // Return lists of hits associated to a track for all measurement layers. 
     209        // Return lists of hits associated to a track for all measurement layers.
    210210        // Return value is the total number of measurement hits
    211211        // kmh = total number of measurement layers hit for given type
    … …  
    234234}
    235235//
     236Int_t SolTrack::FirstHit(Double_t &Xfirst, Double_t &Yfirst, Double_t &Zfirst)
     237{
     238        Int_t iFirst = -1;
     239        Int_t iFirstLay = -1;   // Default return with no hits
     240        Xfirst = 0.;
     241        Yfirst = 0.;
     242        Zfirst = 0.;
     243        Int_t Nmh = nmHit();    // # measurement hits
     244        if(Nmh > 0){
     245                Int_t    *ih = new Int_t   [Nmh];
     246                Double_t *Xh = new Double_t[Nmh];
     247                Double_t *Yh = new Double_t[Nmh];
     248                Double_t *Zh = new Double_t[Nmh];
     249                Double_t *dh = new Double_t[Nmh];
     250                //
     251                Int_t n = HitListXYZ(ih, Xh, Yh, Zh);
     252                //
     253                for(Int_t i=0; i<Nmh; i++){
     254                        Double_t rr = TMath::Sqrt(Xh[i]*Xh[i]+Yh[i]*Yh[i]);     // Hit radius
     255                        dh[i] = TMath::ASin(C() * TMath::Sqrt((rr * rr - D() * D()) / (1. + 2 * C() * D()))) / C();     // Arc length traveled
     256                }
     257                //
     258                Int_t *hord = new Int_t[Nmh];                   // hit order by increasing arc length
     259                TMath::Sort(Nmh, dh, hord, kFALSE);             // Order by increasing arc length
     260                iFirst = hord[0];                               // First hit pointer
     261                Xfirst = Xh[iFirst];
     262                Yfirst = Yh[iFirst];
     263                Zfirst = Zh[iFirst];
     264                iFirstLay = ih[iFirst];
     265                //
     266                // Clean
     267                delete [] ih;
     268                delete [] Xh;
     269                delete [] Yh;
     270                delete [] Zh;
     271                delete [] dh;
     272                delete [] hord;
     273        }
     274//
     275        return  iFirstLay;      // Return first hit layer number
     276}
     277//
    236278// Track plot
    237279//
    … …  
    285327        //
    286328        //
    287         // Input flags: 
     329        // Input flags:
    288330        //                              Res = .TRUE. turn on resolution effects/Use standard resolutions
    289331        //                                        .FALSE. set all resolutions to 0
    … …  
    368410                        nz = 1.0;
    369411                }
    370                 Double_t corr = TMath::Abs(pxi*nx + pyi * ny + pzi * nz) / p(); 
     412                Double_t corr = TMath::Abs(pxi*nx + pyi * ny + pzi * nz) / p();
    371413                Double_t Rlf = fG->lTh(i) / (corr*fG->lX0(i));                                  // Rad. length fraction
    372414                thms[ii] = 0.0136*TMath::Sqrt(Rlf)*(1.0 + 0.038*TMath::Log(Rlf)) / p();         // MS angle
    … …  
    396438                Int_t ityp  = fG->lTyp(i);                      // Layer type Barrel or Z
    397439                Int_t nmeai = fG->lND(i);                       // # measurements in layer
    398                
     440
    399441                if (fG->isMeasure(i))
    400442                {
    … …  
    428470                                                Rm(im, 2) = csa * dRphi(2) - ssa * dRz(2);      // C derivative
    429471                                                Rm(im, 3) = csa * dRphi(3) - ssa * dRz(3);      // z0 derivative
    430                                                 Rm(im, 4) = csa * dRphi(4) - ssa * dRz(4);      // cot(theta) derivative       
     472                                                Rm(im, 4) = csa * dRphi(4) - ssa * dRz(4);      // cot(theta) derivative
    431473                                        }
    432474                                        if (ityp == 2)          // Z type layer (Measure R-phi at const. Z)
    … …  
    489531                                                {
    490532                                                        Double_t strk = 0;
    491                                                         if (nmk + 1 == 1) strk = fG->lStU(k);   // Stereo angle upper 
     533                                                        if (nmk + 1 == 1) strk = fG->lStU(k);   // Stereo angle upper
    492534                                                        if (nmk + 1 == 2) strk = fG->lStL(k);   // Stereo angle lower
    493535                                                        //if (im == km && Res) Sm(im, km) += sig*sig;   // Detector resolution on diagonal
    … …  
    500542                                                        //
    501543                                                        // Loop on all layers below for MS contributions
    502                                                         for (Int_t jj = 0; jj < kk; jj++)               
     544                                                        for (Int_t jj = 0; jj < kk; jj++)
    503545                                                        {
    504546                                                                Double_t di = dik(ii, jj);

  • external/TrackCovariance/SolTrack.h

    r54ec182 r4564cba  
    7373        Int_t HitList(Int_t *&ihh, Double_t *&rhh, Double_t *&zhh);
    7474        Int_t HitListXYZ(Int_t *&ihh, Double_t *&Xh, Double_t *&Yh, Double_t *&Zh);
     75        Int_t FirstHit(Double_t &Xfirst, Double_t &Yfirst, Double_t &Zfirst);   // First hit position
    7576        //
    7677        // Track graph

  • external/TrackCovariance/VertexFit.cc

    r54ec182 r4564cba  
    1717{
    1818        fNtr = 0;
    19         fRold = -1.0;
     19        fRstart = -1.0;
    2020        fVtxDone = kFALSE;
    2121        fVtxCst = kFALSE;
    … …  
    3131{
    3232        fNtr = Ntr;
    33         fRold = -1.0;
     33        fRstart = -1.0;
    3434        fVtxDone = kFALSE;
    3535        fVtxCst = kFALSE;
    … …  
    5555{
    5656        fNtr = Ntr;
    57         fRold = -1.0;
     57        fRstart = -1.0;
    5858        fVtxDone = kFALSE;
    5959        fVtxCst = kFALSE;
    … …  
    220220        std::vector<TVectorD*> x0i;                             // Tracks at ma
    221221        std::vector<TVectorD*> ni;                              // Track derivative wrt phase
    222         std::vector<TMatrixDSym*> Ci;                   // Position error matrix at fixed phase
     222        std::vector<TMatrixDSym*> Ci;                           // Position error matrix at fixed phase
    223223        std::vector<TVectorD*> wi;                              // Ci*ni
     224        std::vector<Double_t> s_in;                             // Starting phase
     225        //
     226        //
    224227        //
    225228        // Track loop
    226229        for (Int_t i = 0; i < fNtr; i++)
    227230        {
    228                 Double_t s = 0.;
    229231                TVectorD par = *fPar[i];
    230232                TMatrixDSym Cov = *fCov[i];
    231                 x0i.push_back(new TVectorD(Fill_x0(par)));
     233                Double_t s = 0.;
     234                // Case when starting radius is provided
     235                if(fRstart > TMath::Abs(par(0))){
     236                        s = 2.*TMath::ASin(par(2)*TMath::Sqrt((fRstart*fRstart-par(0)*par(0))/(1.+2.*par(2)*par(0))));
     237                }
     238                //
     239                x0i.push_back(new TVectorD(Fill_x(par, s)));
    232240                ni.push_back(new TVectorD(derXds(par, s)));
    233241                TMatrixD A = derXdPar(par, s);
    … …  
    235243                TMatrixDSym Cinv = RegInv(*Ci[i]);
    236244                wi.push_back(new TVectorD(Cinv * (*ni[i])));
     245                s_in.push_back(s);
    237246        }
    238247        //std::cout << "Vtx init completed. fNtr = "<<fNtr << std::endl;
    … …  
    262271        for (Int_t i = 0; i < fNtr; i++){
    263272                Double_t si = Dot(*wi[i], fXv - (*x0i[i])) / Ci[i]->Similarity(*wi[i]);
    264                 ffi.push_back(si);
     273                ffi.push_back(si+s_in[i]);
    265274                //TVectorD xvi = Fill_x(*fPar[i],si);
    266275                //std::cout << "Fast vertex "<<i<<": xvi = "<<xvi(0)<<", "<<xvi(1)<<", "<<xvi(2)
    … …  
    280289        Ci.clear();
    281290        wi.clear();
     291        s_in.clear();
    282292}
    283293//
    … …  
    402412        //
    403413        fVtxDone = kTRUE;               // Set fit completion flag
    404         fRold = TMath::Sqrt(fXv(0)*fXv(0) + fXv(1)*fXv(1));     // Store fit
     414        //fRstart = TMath::Sqrt(fXv(0)*fXv(0) + fXv(1)*fXv(1)); // Store fit
    405415        //std::cout << "Found vertex " << fXv(0) << ", " << fXv(1) << ", " << fXv(2)
    406416        //      << ", after "<<Ntry<<" iterations"<<std::endl;

  • external/TrackCovariance/VertexFit.h

    r54ec182 r4564cba  
    3636        // Results
    3737        Bool_t fVtxDone;                        // Flag vertex fit completed
    38         Double_t fRold;                         // Current value of vertex radius
     38        Double_t fRstart;                       // Starting value of vertex radius (0 = none)
    3939        TVectorD fXv;                           // Found vertex
    4040        TMatrixDSym fcovXv;                     // Vertex covariance
    … …  
    8383        // Handle tracks/constraints
    8484        void AddVtxConstraint(TVectorD xv, TMatrixDSym cov);    // Add gaussian vertex constraint
    85         void AddTrk(TVectorD *par, TMatrixDSym *Cov);                           // Add track to input list
    86         void RemoveTrk(Int_t iTrk);                                                             // Remove iTrk track
     85        void AddTrk(TVectorD *par, TMatrixDSym *Cov);           // Add track to input list
     86        void RemoveTrk(Int_t iTrk);                             // Remove iTrk track
     87        void SetStartR(Double_t R) { fRstart = R; };            // Set starting radius
    8788        //
    8889};

  • modules/ModulesLinkDef.h

    r54ec182 r4564cba  
    7979#include "modules/TruthVertexFinder.h"
    8080#include "modules/ExampleModule.h"
     81#include "modules/LLPFilter.h"
     82#include "modules/CscClusterEfficiency.h"
     83#include "modules/CscClusterId.h"
    8184
    8285#ifdef __CINT__
    … …  
    139142#pragma link C++ class TruthVertexFinder+;
    140143#pragma link C++ class ExampleModule+;
     144#pragma link C++ class LLPFilter+;
     145#pragma link C++ class CscClusterEfficiency+;
     146#pragma link C++ class CscClusterId+;
    141147
    142148#endif

  • modules/TrackCovariance.cc

    r54ec182 r4564cba  
    121121    const TLorentzVector &candidateMomentum = particle->Momentum;
    122122
     123
    123124    Bool_t inside = TrkUtil::IsInside(candidatePosition.Vect(), Rin, ZinNeg, ZinPos); // Check if in inner box
    124125    Bool_t Accept = kTRUE;
    … …  
    130131
    131132    ObsTrk track(candidatePosition.Vect(), candidateMomentum.Vect(), candidate->Charge, fCovariance, fGeometry);
     133
    132134
    133135    mother    = candidate;
    … …  
    150152    candidate->Yd = track.GetObsX().Y()*1e03;
    151153    candidate->Zd = track.GetObsX().Z()*1e03;
     154
     155    candidate->XFirstHit = track.GetFirstHit().X()*1e03;
     156    candidate->YFirstHit = track.GetFirstHit().Y()*1e03;
     157    candidate->ZFirstHit = track.GetFirstHit().Z()*1e03;
    152158
    153159    candidate->D0       = track.GetObsPar()[0]*1e03;

  • modules/TreeWriter.cc

    r54ec182 r4564cba  
    4545#include "TString.h"
    4646
     47#include <set>
    4748#include <algorithm>
    4849#include <iostream>
    … …  
    7677  fClassMap[Electron::Class()] = &TreeWriter::ProcessElectrons;
    7778  fClassMap[Muon::Class()] = &TreeWriter::ProcessMuons;
     79  fClassMap[CscCluster::Class()] = &TreeWriter::ProcessCscCluster;
    7880  fClassMap[Jet::Class()] = &TreeWriter::ProcessJets;
    7981  fClassMap[MissingET::Class()] = &TreeWriter::ProcessMissingET;
    … …  
    149151{
    150152  TIter it1(candidate->GetCandidates());
     153  set<Candidate *> s;
     154  set<Candidate *>::iterator it3;
    151155  it1.Reset();
     156  s.clear();
    152157  array->Clear();
    153158
    … …  
    159164    if(candidate->GetCandidates()->GetEntriesFast() == 0)
    160165    {
    161       array->Add(candidate);
     166      s.insert(candidate);
    162167      continue;
    163168    }
    … …  
    167172    if(candidate->GetCandidates()->GetEntriesFast() == 0)
    168173    {
    169       array->Add(candidate);
     174      s.insert(candidate);
    170175      continue;
    171176    }
    … …  
    175180    while((candidate = static_cast<Candidate *>(it2.Next())))
    176181    {
    177       array->Add(candidate->GetCandidates()->At(0));
     182      candidate = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
     183      if(candidate->GetCandidates()->GetEntriesFast() == 0)
     184      {
     185        s.insert(candidate);
     186      }
    178187    }
     188  }
     189
     190  for(it3 = s.begin(); it3 != s.end(); ++it3)
     191  {
     192    array->Add(*it3);
    179193  }
    180194}
    … …  
    238252    entry->Z = position.Z();
    239253    entry->T = position.T() * 1.0E-3 / c_light;
     254
    240255  }
    241256}
    … …  
    384399    entry->Zd = candidate->Zd;
    385400
     401    entry->XFirstHit = candidate->XFirstHit;
     402    entry->YFirstHit = candidate->YFirstHit;
     403    entry->ZFirstHit = candidate->ZFirstHit;
     404
    386405    const TLorentzVector &momentum = candidate->Momentum;
    387406
    … …  
    544563    entry->Zd = candidate->Zd;
    545564
     565    entry->XFirstHit = candidate->XFirstHit;
     566    entry->YFirstHit = candidate->YFirstHit;
     567    entry->ZFirstHit = candidate->ZFirstHit;
     568
    546569    const TLorentzVector &momentum = candidate->Momentum;
    547570
    … …  
    886909  }
    887910}
     911//------------------------------------------------------------------------------
     912
     913void TreeWriter::ProcessCscCluster(ExRootTreeBranch *branch, TObjArray *array)
     914{
     915  TIter iterator(array);
     916  Candidate *candidate = 0;
     917  CscCluster *entry = 0;
     918  Double_t pt, signPz, cosTheta, eta, rapidity;
     919
     920  const Double_t c_light = 2.99792458E8; // in unit of m/s
     921
     922  array->Sort();
     923
     924
     925  // loop over all clusters
     926  iterator.Reset();
     927  while((candidate = static_cast<Candidate *>(iterator.Next())))
     928  {
     929    const TLorentzVector &momentum = candidate->Momentum;
     930    const TLorentzVector &position = candidate->DecayPosition;
     931
     932    pt = momentum.Pt();
     933    cosTheta = TMath::Abs(momentum.CosTheta());
     934    signPz = (momentum.Pz() >= 0.0) ? 1.0 : -1.0;
     935    eta = (cosTheta == 1.0 ? signPz * 999.9 : momentum.Eta());
     936
     937    entry = static_cast<CscCluster *>(branch->NewEntry());
     938
     939    entry->SetBit(kIsReferenced);
     940    entry->SetUniqueID(candidate->GetUniqueID());
     941
     942    entry->Eta = eta;
     943    entry->Phi = momentum.Phi();
     944
     945    entry->PT = momentum.Pt(); // pt of LLP
     946    entry->Px = momentum.Px();// px of LLP
     947    entry->Py = momentum.Py();// py of LLP
     948    entry->Pz = momentum.Pz();// pz of LLP
     949    entry->E = momentum.E(); // E of LLP
     950    entry->pid = candidate->PID; // LLP pid
     951    entry->Eem = candidate->Eem; // LLP Eem
     952    entry->Ehad = candidate->Ehad; // LLP Ehad
     953    Double_t beta = momentum.P()/momentum.E();
     954    Double_t gamma = 1.0/sqrt(1-beta*beta);
     955    Double_t decayDistance = sqrt(pow(position.X(),2)+pow(position.Y(),2)+pow(position.Z(),2)); // mm
     956    entry->beta = beta; // LLP pid
     957    entry->ctau = decayDistance/(beta * gamma); // LLP travel time in its rest frame
     958    entry->T = decayDistance*(1./beta-1)* 1.0E-3/c_light*1e9; // ns
     959    entry->X = position.X(); // LLP decay x
     960    entry->Y = position.Y(); //  LLP decay y
     961    entry->Z = position.Z(); //  LLP decay z
     962  }
     963}
    888964
    889965//------------------------------------------------------------------------------

  • modules/TreeWriter.h

    r54ec182 r4564cba  
    6060  void ProcessElectrons(ExRootTreeBranch *branch, TObjArray *array);
    6161  void ProcessMuons(ExRootTreeBranch *branch, TObjArray *array);
     62  void ProcessCscCluster(ExRootTreeBranch *branch, TObjArray *array);
    6263  void ProcessTauJets(ExRootTreeBranch *branch, TObjArray *array);
    6364  void ProcessJets(ExRootTreeBranch *branch, TObjArray *array);
    … …  
    6768  void ProcessWeight(ExRootTreeBranch *branch, TObjArray *array);
    6869  void ProcessHectorHit(ExRootTreeBranch *branch, TObjArray *array);
     70
    6971
    7072#if !defined(__CINT__) && !defined(__CLING__)
Note: See TracChangeset for help on using the changeset viewer.