200 likes | 395 Views
Physics Analysis with ROOT. Yajun Mao School of Physics, Peking University. July 27, 2005. Outline. Physics Analysis in General An New Approach Examples as in HERMES Summary. Physics Analysis in General. 4-mom, 4-vertex, pid, tof,…. Global tracks. PIDs. More detailed info for pid.
E N D
Physics Analysis with ROOT Yajun Mao School of Physics, Peking University July 27, 2005 BESIII Software Meeting,IHEP
Outline Physics Analysis in General An New Approach Examples as in HERMES Summary BESIII Software Meeting,IHEP
Physics Analysis in General 4-mom, 4-vertex, pid, tof,… Global tracks PIDs More detailed info for pid Tracks in Sub-detectors Partial tracks in sub-det… DST or mDST Hits in sub-detectors Hits in sub-detecors … … … Beam status, lum, DAQ status,… Beam、DAQ… Run# Mother particle info … … Event# NTUPLE Grand-Mother particle info … … Physics Analysis Particle#1 info Particle#2 info … … Particle#N info BESIII Software Meeting,IHEP
What’s the Problem? Straightforward between DST & NTUPLE Only good for simple case, what if we need some additional information? what if we want to check other accompany particles? No good means to find some deep hiding bugs 1 6 2 3 7 4 5 Broken data structure Easy to make mistake Hard to share analysis code, low efficiency… BESIII Software Meeting,IHEP
A Possible Solution Split the analysis to several steps Take the analysis as event filtering Keep the whole event in each steps, not only those selected tracks Use a consistent data structure Code re-usable to other analysis Less bugs for more people’s use Higher efficiency, could concentrate in physics, rather than programming BESIII Software Meeting,IHEP
A track is a measured particle such as particle 1, 2 ,3 ,4 ,5 but a particle doesn’t to be a track for example, particle 6, 7 a photon is not a track too 1 6 2 3 7 4 5 From Track To Particle? Define a class for particle Cut away the connection to a special experiment More general sharing BESIII Software Meeting,IHEP
A Class for Particles: HParticle TObject TAttLine TAtt3D TParticle HParticle BESIII Software Meeting,IHEP
The TParticle Class(1) Int_t fPdgCode; // PDG code of the particle Int_t fStatusCode; // generation status code Int_t fMother[2]; // Indices of the mother particles Int_t fDaughter[2]; // Indices of the daughter particles Float_t fWeight; // particle weight Double_t fCalcMass; // Calculated mass Double_t fPx; // x component of momentum Double_t fPy; // y component of momentum Double_t fPz; // z component of momentum Double_t fE; // Energy Double_t fVx; // x of production vertex Double_t fVy; // y of production vertex Double_t fVz; // z of production vertex Double_t fVt; // t of production vertex Double_t fPolarTheta; // Polar angle of polarisation Double_t fPolarPhi; // azymutal angle of polarisation TParticlePDG* fParticlePDG; //! reference to the particle record in PDG database BESIII Software Meeting,IHEP
The TParticle Class(2) Int_t GetStatusCode () const { return fStatusCode; } Int_t GetPdgCode () const { return fPdgCode; } Int_t GetFirstMother () const { return fMother[0]; } Int_t GetMother (Int_t i) const { return fMother[i]; } Int_t GetSecondMother () const { return fMother[1]; } Int_t GetFirstDaughter() const { return fDaughter[0]; } Int_t GetDaughter (Int_t i) const { return fDaughter[i]; } Int_t GetLastDaughter () const { return fDaughter[1]; } Double_t GetCalcMass () const { return fCalcMass; } Double_t GetMass () { return GetPDG()->Mass();} Int_t GetNDaughters () const { return fDaughter[1]>0 ? fDaughter[1]-fDaughter[0]+1 : 0;} void GetPolarisation(TVector3 &v); Float_t GetWeight () const { return fWeight;} TParticlePDG* GetPDG (Int_t mode = 0); Int_t Beauty () { return GetPDG()->Beauty(); } Int_t Charm () { return GetPDG()->Charm(); } Int_t Strangeness () { return GetPDG()->Strangeness();} void Momentum(TLorentzVector &v) { v.SetPxPyPzE(fPx,fPy,fPz,fE);} void ProductionVertex(TLorentzVector &v) { v.SetXYZT(fVx,fVy,fVz,fVt);} // ****** redefine several most oftenly used // methods of LORENTZ_VECTOR Double_t Vx () const { return fVx;} Double_t Vy () const { return fVy;} Double_t Vz () const { return fVz;} Double_t T () const { return fT;} … … BESIII Software Meeting,IHEP
The HParticle Class(1) Float_t fLifeTime; // Life time of this particle TVector3 fDecayVertex; // the decay vertex of the particle Float_t fDaughterDCA; // the closest approach of 2 daughter tracks Int_t fTrackIndex; // index of track associated to the particle Int_t fIndex; // index of this particle Int_t fAcceptanceOK; // if the track is in acceptance; Float_t fLengthOfFlight; // The length of flight before it decays Int_t fDaughterIndex[3]; // indices of Daughters; Float_t fRecoMassWin[2]; // Reconstructed mass windows(in GeV) Int_t fDaughterAllIndex[15]; // indices including grand daughters; Double_t fUserVariable[10]; // User Defined Variables; BESIII Software Meeting,IHEP
The HParticle Class(2) //!-------------- reconstruction related ---------------- // The particle is reconstructed from its daugthers d1, d2, void Reconstruction(HParticle* d1, HParticle* d2); // The particle is reconstructed from its daugthers d1, d2, d3 void Reconstruction(HParticle* d1, HParticle* d2, HParticle* d3); // Set a window for reconstructed mass void SetRecoMassWin(Float_t xmin, Float_t xmax) {fRecoMassWin[0]=xmin;fRecoMassWin[1]=xmax; } // Get the window for reconstructed mass Float_t GetRecoMassWin(Int_t i) { return fRecoMassWin[i]; } // find the distance of closest approach with other particle Float_t DCA(HParticle *d, TVector3 &pos); // Get index of i-th all daughters of this particle // GetDaugtherAllIndex(0) returns the number of all daughters Int_t GetDaughterAllIndex(Int_t i) { return fDaughterAllIndex[i]; } void SetDaughterAllIndex(Int_t i, Int_t idx){fDaughterAllIndex[i] = idx; } // Check if particle d exists in this particle Bool_t Exist(HParticle *d); //!-------------- decay related ---------------- // Simulate the particle decays to its daughters d1, d2 void Decay(HParticle *d1, HParticle *d2); …… BESIII Software Meeting,IHEP
How To Use HParticle In case of simulation //define lambda TVector3 p(px,py,pz); TVector3 v(vx,vy,vz); HParticle *Lambda = new HParticle(3212,p,v); // define a proton HParticle *Proton = new HParticle(2212); // define a pion- HParticle *PiMinus = new HParticle(-211); // simulate decay Lambda->Decay(Proton,PiMinus) In case of reconstruction //define proton TVector3 ProtonP(px1,py1,pz1); TVector3 ProtonV(vx1,vy1,vz1); HParticle *Proton = new HParticle(2212,ProtonP,ProtonV); //define pion- TVector3 PiMinusP(px2,py2,pz2); TVector3 PiMinusV(vx2,vy2,vz2); HParticle *PiMinus = new HParticle(-211,PiMinusP,PiMinusV); //define lambda HParticle *Lambda = new HParticle(3212); // reconstruct lambda from proton and pion- Lambda->Reconstruction(Proton,PiMinus) BESIII Software Meeting,IHEP
Physical Event Reconstruction Apply Rough Cut ROOT TTree DST or mDST DST Structure HParticle Container Header Part. 1 Part. N DST Structure Only one HParticle Object In TClonesArray Drawable Drawable Physics Analysis with HParticle ~200GB ~2GB ~20MB BESIII Software Meeting,IHEP
HERMES mDST(ADAMO Tables) g1Track g1Beam g1DAQ smRICH … … g1DAQ smRICH g1Target ROOT Ttree Burst g1Beam g1Target Event … …. … … g1Track smTrack … …. HERMES Data Flow ROOT Converter Event Reconstruction For a physical process TTree Event g1Track smRICH smTrack Particles … … BESIII Software Meeting,IHEP
PPbar Selection // create an instance of ReconHERMES ReconHERMES *process = new ReconHERMES(); // create dummy particles you have interest to study HParticle PPbar(9003); HParticle Proton(2212); HParticle Pbar(-2212); PPbar.SetRecoMassWin(1.8,3.0); //Set Mass window for PPbar // Loop over all events in this burst for(Int_t j=0; j<Nevt; j++){ process->SetEventNum(j); // set to current event process->TrackToHParticle(Particles); // convert tracks to HParticle objects // only process those events which contain all necessary final states if(process->EventOK(Particles,Nfinal,finalPDG)){ TOrdCollection PPbars; // reconstruct PPbar from Particles and put result in PPbars process->Reconstruction(Particles,&PPbars,&PPbar,&Proton,&Pbar); // need at least 1 PPbar for this analysis if(PPbars.GetSize()>0){ Particles->AddAll(&PPbars); //save Particles here } } BESIII Software Meeting,IHEP
Structure After PPbar Selection BESIII Software Meeting,IHEP
Rough Cut and NTuples // create a dummy particle you have interest to study HParticle PPbar(9003); // create an instance of ReconHERMES ReconHERMES *process = new ReconHERMES(); for(Int_t i=0; i<Tot1; i++){ T->GetEntry(i); Int_t NParticle = Particles->GetEntries(); Int_t NTrack = g1Tracks->GetEntries(); // Find associate information to be saved for (Int_t j=NTrack; j<NParticle; j++){ HParticle *fPPbar = (HParticle *)Particles->At(j); if( fPPbar->GetPdgCode() == fakepart.GetPdgCode()){ // find one candidate Int_t index_proton = fPPbar->GetDaughterIndex(0); Int_t index_pbar= fPPbar->GetDaughterIndex(1); HParticle *fProton = (HParticle *)Particles->At(index_proton); HParticle *fPbar = (HParticle *)Particles->At(index_pbar); // put cut here // save data } } BESIII Software Meeting,IHEP
Structure After Rough Cut BESIII Software Meeting,IHEP
Highlights of HParticle 1) Inherit from ROOT TParticle class 2) Use PDG data base for particle property 3) Use PDG numbering scheme 4) Include its vertex in 4-D phase space 5) Include its momentum in 4-D phase space 6) Could handle a new particle 7) Could be reconstructed from its 2/3-body daughters 8) Could decay to its 2/3-body daughters 9) In case of reconstruction, provide daughter DCA 10) In case of decay, provide track intersected with z=z0 plane BESIII Software Meeting,IHEP
Summary A new analysis procedure was developed based on HParticle class Could handle very complicated physical process in a rather clear/simply way Same data structure for different processes Could check data in each step with ROOT browser, easy to find bugs Easy to add additional info in the events Most of the codes can be shared Simple/clear naming scheme for variables BESIII Software Meeting,IHEP