1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 //_________________________________________________________________________
18 // Class for the analysis of particle - hadron correlations
19 // Particle (for example direct gamma) must be found in a previous analysis
20 //-- Author: Gustavo Conesa (LNF-INFN)
21 //////////////////////////////////////////////////////////////////////////////
24 // --- ROOT system ---
25 #include "Riostream.h"
28 //---- ANALYSIS system ----
30 #include "AliNeutralMesonSelection.h"
31 #include "AliAnaParticleHadronCorrelation.h"
32 #include "AliCaloTrackReader.h"
34 ClassImp(AliAnaParticleHadronCorrelation)
37 //____________________________________________________________________________
38 AliAnaParticleHadronCorrelation::AliAnaParticleHadronCorrelation() :
39 AliAnaPartCorrBaseClass(),
40 fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.),
41 fhPhiCharged(0), fhPhiNeutral(0), fhEtaCharged(0), fhEtaNeutral(0),
42 fhDeltaPhiCharged(0), fhDeltaPhiNeutral(0),
43 fhDeltaEtaCharged(0), fhDeltaEtaNeutral(0),
44 fhDeltaPhiChargedPt(0), fhDeltaPhiNeutralPt(0),
45 fhPtImbalanceNeutral(0), fhPtImbalanceCharged(0)
49 //Initialize parameters
53 //____________________________________________________________________________
54 AliAnaParticleHadronCorrelation::AliAnaParticleHadronCorrelation(const AliAnaParticleHadronCorrelation & g) :
55 AliAnaPartCorrBaseClass(g),
56 fDeltaPhiMaxCut(g.fDeltaPhiMaxCut), fDeltaPhiMinCut(g.fDeltaPhiMinCut),
57 fhPhiCharged(g.fhPhiCharged), fhPhiNeutral(g.fhPhiNeutral),
58 fhEtaCharged(g.fhEtaCharged), fhEtaNeutral(g.fhEtaNeutral),
59 fhDeltaPhiCharged(g.fhDeltaPhiCharged),
60 fhDeltaPhiNeutral(g.fhDeltaPhiNeutral),
61 fhDeltaEtaCharged(g.fhDeltaEtaCharged),
62 fhDeltaEtaNeutral(g.fhDeltaEtaNeutral),
63 fhDeltaPhiChargedPt(g.fhDeltaPhiChargedPt),
64 fhDeltaPhiNeutralPt(g.fhDeltaPhiNeutralPt),
65 fhPtImbalanceNeutral(g.fhPtImbalanceNeutral),
66 fhPtImbalanceCharged(g.fhPtImbalanceCharged)
72 //_________________________________________________________________________
73 AliAnaParticleHadronCorrelation & AliAnaParticleHadronCorrelation::operator = (const AliAnaParticleHadronCorrelation & source)
75 // assignment operator
77 if(this == &source)return *this;
78 ((AliAnaPartCorrBaseClass *)this)->operator=(source);
80 fDeltaPhiMaxCut = source.fDeltaPhiMaxCut ;
81 fDeltaPhiMinCut = source.fDeltaPhiMinCut ;
83 fhPhiCharged = source.fhPhiCharged ; fhPhiNeutral = source.fhPhiNeutral ;
84 fhEtaCharged = source.fhEtaCharged ; fhEtaNeutral = source.fhEtaNeutral ;
85 fhDeltaPhiCharged = source.fhDeltaPhiCharged ;
86 fhDeltaPhiNeutral = source.fhDeltaPhiNeutral ;
87 fhDeltaPhiNeutralPt = source.fhDeltaPhiNeutralPt ;
88 fhDeltaEtaCharged = source.fhDeltaEtaCharged ;
89 fhDeltaEtaNeutral = source.fhDeltaEtaNeutral ;
90 fhDeltaPhiChargedPt = source.fhDeltaPhiChargedPt ;
92 fhPtImbalanceNeutral = source.fhPtImbalanceNeutral ;
93 fhPtImbalanceCharged = source.fhPtImbalanceCharged ;
99 //________________________________________________________________________
100 TList * AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
103 // Create histograms to be saved in output file and
104 // store them in fOutputContainer
105 TList * outputContainer = new TList() ;
106 outputContainer->SetName("CorrelationHistos") ;
108 //Correlation with charged hadrons
109 if(GetReader()->IsCTSSwitchedOn()) {
110 fhPhiCharged = new TH2F
111 ("PhiCharged","#phi_{h^{#pm}} vs p_{T trigger}",
113 fhPhiCharged->SetYTitle("#phi_{h^{#pm}} (rad)");
114 fhPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
116 fhEtaCharged = new TH2F
117 ("EtaCharged","#eta_{h^{#pm}} vs p_{T trigger}",
119 fhEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");
120 fhEtaCharged->SetXTitle("p_{T trigger} (GeV/c)");
122 fhDeltaPhiCharged = new TH2F
123 ("DeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",
124 200,0,120,200,0,6.4);
125 fhDeltaPhiCharged->SetYTitle("#Delta #phi");
126 fhDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
128 fhDeltaPhiChargedPt = new TH2F
129 ("DeltaPhiChargedPt","#phi_{trigger} - #phi_{#p^{#pm}i} vs p_{T h^{#pm}}",
130 200,0,120,200,0,6.4);
131 fhDeltaPhiChargedPt->SetYTitle("#Delta #phi");
132 fhDeltaPhiChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
134 fhDeltaEtaCharged = new TH2F
135 ("DeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs p_{T trigger}",
137 fhDeltaEtaCharged->SetYTitle("#Delta #eta");
138 fhDeltaEtaCharged->SetXTitle("p_{T trigger} (GeV/c)");
140 fhPtImbalanceCharged =
141 new TH2F("CorrelationCharged","z_{trigger h^{#pm}} = p_{T h^{#pm}} / p_{T trigger}",
142 240,0.,120.,1000,0.,1.2);
143 fhPtImbalanceCharged->SetYTitle("z_{trigger h^{#pm}}");
144 fhPtImbalanceCharged->SetXTitle("p_{T trigger}");
146 outputContainer->Add(fhPhiCharged) ;
147 outputContainer->Add(fhEtaCharged) ;
148 outputContainer->Add(fhDeltaPhiCharged) ;
149 outputContainer->Add(fhDeltaEtaCharged) ;
150 outputContainer->Add(fhPtImbalanceCharged) ;
151 outputContainer->Add(fhDeltaPhiChargedPt) ;
152 } //Correlation with charged hadrons
154 //Correlation with neutral hadrons
155 if(GetReader()->IsEMCALSwitchedOn() || GetReader()->IsPHOSSwitchedOn()){
157 fhPhiNeutral = new TH2F
158 ("PhiNeutral","#phi_{#pi^{0}} vs p_{T trigger}",
160 fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
161 fhPhiNeutral->SetXTitle("p_{T trigger} (GeV/c)");
163 fhEtaNeutral = new TH2F
164 ("EtaNeutral","#eta_{#pi^{0}} vs p_{T trigger}",
166 fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
167 fhEtaNeutral->SetXTitle("p_{T trigger} (GeV/c)");
169 fhDeltaPhiNeutral = new TH2F
170 ("DeltaPhiNeutral","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T trigger}",
171 200,0,120,200,0,6.4);
172 fhDeltaPhiNeutral->SetYTitle("#Delta #phi");
173 fhDeltaPhiNeutral->SetXTitle("p_{T trigger} (GeV/c)");
175 fhDeltaPhiNeutralPt = new TH2F
176 ("DeltaPhiNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",
177 200,0,120,200,0,6.4);
178 fhDeltaPhiNeutralPt->SetYTitle("#Delta #phi");
179 fhDeltaPhiNeutralPt->SetXTitle("p_{T trigger} (GeV/c)");
181 fhDeltaEtaNeutral = new TH2F
182 ("DeltaEtaNeutral","#eta_{trigger} - #eta_{#pi^{0}} vs p_{T trigger}",
184 fhDeltaEtaNeutral->SetYTitle("#Delta #eta");
185 fhDeltaEtaNeutral->SetXTitle("p_{T trigger} (GeV/c)");
187 fhPtImbalanceNeutral =
188 new TH2F("CorrelationNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger}",
189 240,0.,120.,1000,0.,1.2);
190 fhPtImbalanceNeutral->SetYTitle("z_{trigger #pi^{0}}");
191 fhPtImbalanceNeutral->SetXTitle("p_{T trigger}");
193 outputContainer->Add(fhPhiNeutral) ;
194 outputContainer->Add(fhEtaNeutral) ;
195 outputContainer->Add(fhDeltaPhiNeutral) ;
196 outputContainer->Add(fhDeltaEtaNeutral) ;
197 outputContainer->Add(fhPtImbalanceNeutral) ;
200 //Keep neutral meson selection histograms if requiered
201 //Setting done in AliNeutralMesonSelection
202 if(GetNeutralMesonSelection()){
203 TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
204 cout<<"NMSHistos "<< nmsHistos<<endl;
205 if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
206 for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
209 }//Correlation with neutral hadrons
211 return outputContainer;
214 //____________________________________________________________________________
215 void AliAnaParticleHadronCorrelation::InitParameters()
218 //Initialize the parameters of the analysis.
220 SetPtCutRange(2,300);
221 fDeltaPhiMinCut = 1.5 ;
222 fDeltaPhiMaxCut = 4.5 ;
226 //__________________________________________________________________
227 void AliAnaParticleHadronCorrelation::Print(const Option_t * opt) const
230 //Print some relevant parameters set for the analysis
234 Info("*** Print *** ", "%s %s", GetName(), GetTitle() ) ;
235 printf("pT Hadron > %2.2f\n", GetMinPt()) ;
236 printf("pT Hadron < %2.2f\n", GetMaxPt()) ;
237 printf("Phi trigger particle-Hadron < %3.2f\n", fDeltaPhiMaxCut) ;
238 printf("Phi trigger particle-Hadron > %3.2f\n", fDeltaPhiMinCut) ;
243 //____________________________________________________________________________
244 void AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD()
246 //Particle-Hadron Correlation Analysis, fill AODs
248 printf("Begin hadron correlation analysis, fill AODs \n");
249 printf("In particle branch aod entries %d\n", GetAODBranch()->GetEntries());
250 printf("In CTS aod entries %d\n", GetAODCTS()->GetEntries());
251 printf("In EMCAL aod entries %d\n", GetAODEMCAL()->GetEntries());
252 printf("In PHOS aod entries %d\n", GetAODPHOS()->GetEntries());
255 //Loop on stored AOD particles, trigger
256 Int_t naod = GetAODBranch()->GetEntriesFast();
257 for(Int_t iaod = 0; iaod < naod ; iaod++){
258 AliAODParticleCorrelation* particle = dynamic_cast<AliAODParticleCorrelation*> (GetAODBranch()->At(iaod));
259 //Make correlation with charged hadrons
260 if(GetReader()->IsCTSSwitchedOn() )
261 MakeChargedCorrelation(particle, (TSeqCollection*)GetAODCTS(),kFALSE);
263 //Make correlation with neutral pions
264 //Trigger particle in PHOS, correlation with EMCAL
265 if(particle->GetDetector()=="PHOS" && GetReader()->IsEMCALSwitchedOn() && GetAODEMCAL()->GetEntries() > 0)
266 MakeNeutralCorrelation(particle,(TSeqCollection*)GetAODEMCAL(),kFALSE);
267 //Trigger particle in EMCAL, correlation with PHOS
268 else if(particle->GetDetector()=="EMCAL" && GetReader()->IsPHOSSwitchedOn() && GetAODPHOS()->GetEntries() > 0)
269 MakeNeutralCorrelation(particle,(TSeqCollection*)GetAODPHOS(),kFALSE);
270 //Trigger particle in CTS, correlation with PHOS, EMCAL and CTS
271 else if(particle->GetDetector()=="CTS" ){
272 if(GetReader()->IsPHOSSwitchedOn() && GetAODPHOS()->GetEntries() > 0)
273 MakeNeutralCorrelation(particle,(TSeqCollection*)GetAODPHOS(),kFALSE);
274 if(GetReader()->IsEMCALSwitchedOn() && GetAODEMCAL()->GetEntries() > 0)
275 MakeNeutralCorrelation(particle,(TSeqCollection*)GetAODEMCAL(),kFALSE);
280 if(GetDebug() > 1) printf("End hadron correlation analysis, fill AODs \n");
284 //____________________________________________________________________________
285 void AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
287 //Particle-Hadron Correlation Analysis, fill AODs
289 printf("Begin hadron correlation analysis, fill histograms \n");
290 printf("In particle branch aod entries %d\n", GetAODBranch()->GetEntries());
293 //Loop on stored AOD particles
294 Int_t naod = GetAODBranch()->GetEntriesFast();
295 for(Int_t iaod = 0; iaod < naod ; iaod++){
296 AliAODParticleCorrelation* particle = dynamic_cast<AliAODParticleCorrelation*> (GetAODBranch()->At(iaod));
299 printf("Particle %d, In Track Refs entries %d\n", iaod, (particle->GetRefTracks())->GetEntries());
300 printf("Particle %d, In Cluster Refs entries %d\n",iaod, (particle->GetRefClusters())->GetEntries());
303 //Make correlation with charged hadrons
304 if((particle->GetRefTracks())->GetEntries() > 0)
305 MakeChargedCorrelation(particle, (TSeqCollection*) (particle->GetRefTracks()),kTRUE);
307 //Make correlation with neutral pions
308 if((particle->GetRefClusters())->GetEntries() > 0)
309 MakeNeutralCorrelation(particle, (TSeqCollection*) (particle->GetRefClusters()), kTRUE);
313 if(GetDebug() > 1) printf("End hadron correlation analysis, fill histograms \n");
317 //____________________________________________________________________________
318 void AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODParticleCorrelation *aodParticle, TSeqCollection* pl, const Bool_t bFillHisto)
320 // Charged Hadron Correlation Analysis
321 if(GetDebug() > 1)printf("Make trigger particle - charged hadron correlation \n");
323 Double_t ptTrig = aodParticle->Pt();
324 Double_t phiTrig = aodParticle->Phi();
326 Double_t rat = -100.;
327 Double_t phi = -100. ;
328 Double_t eta = -100. ;
331 //Track loop, select tracks with good pt, phi and fill AODs or histograms
332 for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
333 AliAODTrack * track = dynamic_cast<AliAODTrack *>(pl->At(ipr)) ;
334 track->GetPxPyPz(p) ;
335 TLorentzVector mom(p[0],p[1],p[2],0);
339 if(phi<0) phi+=TMath::TwoPi();
342 if(IsFidutialCutOn()){
343 Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,"CTS") ;
347 //Select only hadrons in pt range
348 if(pt < GetMinPt() || pt > GetMaxPt()) continue ;
351 printf("charged hadron: pt %f, phi %f, phi trigger %f. Cuts: delta phi min %2.2f, max%2.2f, pT min %2.2f \n",
352 pt,phi,phiTrig,fDeltaPhiMinCut, fDeltaPhiMaxCut, GetMinPt());
356 fhEtaCharged->Fill(ptTrig,eta);
357 fhPhiCharged->Fill(ptTrig,phi);
358 fhDeltaEtaCharged->Fill(ptTrig,aodParticle->Eta()-eta);
359 fhDeltaPhiCharged->Fill(ptTrig,phiTrig-phi);
360 fhDeltaPhiChargedPt->Fill(pt,phiTrig-phi);
361 //Selection within angular range
362 if(((phiTrig-phi)> fDeltaPhiMinCut) && ((phiTrig-phi)<fDeltaPhiMaxCut) ){
363 if(GetDebug() > 2 ) printf("Selected charge for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f ",pt,phi,eta);
364 fhPtImbalanceCharged->Fill(ptTrig,rat);
369 aodParticle->AddTrack(track);
375 //____________________________________________________________________________
376 void AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODParticleCorrelation * aodParticle,TSeqCollection* pl, const Bool_t bFillHisto)
378 // Neutral Pion Correlation Analysis
379 if(GetDebug() > 1) printf("Make trigger particle - neutral hadron correlation \n");
382 Double_t rat = -100.;
383 Double_t phi = -100. ;
384 Double_t eta = -100. ;
385 Double_t ptTrig = aodParticle->Pt();
386 Double_t phiTrig = aodParticle->Phi();
387 Double_t etaTrig = aodParticle->Eta();
389 TLorentzVector gammai;
390 TLorentzVector gammaj;
392 Double_t vertex[] = {0,0,0};
394 if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(vertex);
396 //Cluster loop, select pairs with good pt, phi and fill AODs or histograms
397 for(Int_t iclus = 0;iclus < pl->GetEntries() ; iclus ++ ){
398 AliAODCaloCluster * calo = dynamic_cast< AliAODCaloCluster *>(pl->At(iclus)) ;
401 //Cluster selection, not charged, with photon or pi0 id and in fidutial cut
403 if(!SelectCluster(calo, vertex, gammai, pdg)) continue ;
406 printf("neutral cluster: pt %f, phi %f, phi trigger %f. Cuts: delta phi min %2.2f, max%2.2f, pT min %2.2f \n",
407 gammai.Pt(),gammai.Phi(),phiTrig,fDeltaPhiMinCut, fDeltaPhiMaxCut, GetMinPt());
409 //2 gamma overlapped, found with PID
410 if(pdg == AliCaloPID::kPi0){
412 //Select only hadrons in pt range
413 if(gammai.Pt() < GetMinPt() || gammai.Pt() > GetMaxPt()) continue ;
415 aodParticle->AddCluster(calo);
416 if(GetDebug() > 2) printf("Correlated with selected pi0 (pid): pt %f, phi %f",gammai.Pt(),gammai.Phi());
420 //Make invariant mass analysis
421 else if(pdg == AliCaloPID::kPhoton){
422 //Search the photon companion in case it comes from a Pi0 decay
423 //Apply several cuts to select the good pair;
424 for(Int_t jclus = iclus+1; jclus < pl->GetEntries() ; jclus ++ ){
425 AliAODCaloCluster * calo2 = dynamic_cast< AliAODCaloCluster *>(pl->At(jclus)) ;
427 //Cluster selection, not charged with photon or pi0 id and in fidutial cut
429 if(!SelectCluster(calo2,vertex, gammaj, pdgj)) continue ;
431 if(pdgj == AliCaloPID::kPhoton ){
433 if((gammai+gammaj).Pt() < GetMinPt() || (gammai+gammaj).Pt() > GetMaxPt()) continue ;
435 //Select good pair (aperture and invariant mass)
436 if(GetNeutralMesonSelection()->SelectPair(gammai, gammaj)){
438 if(GetDebug() > 2 ) printf("Neutral Hadron Correlation: AOD Selected gamma pair: pt %2.2f, phi %2.2f, eta %2.2f, M %2.3f",
439 (gammai+gammaj).Pt(),(gammai+gammaj).Phi(),(gammai+gammaj).Eta(), (gammai+gammaj).M());
440 Int_t labels[]={calo->GetLabel(0),calo2->GetLabel(0)};
441 Float_t pid[]={0,0,0,0,0,0,1,0,0,0,0,0,0};//Pi0 weight 1
442 Float_t pos[]={(gammai+gammaj).X(), (gammai+gammaj).Y(), (gammai+gammaj).Z()};
444 AliAODCaloCluster *caloCluster = new AliAODCaloCluster(0,2,labels,(gammai+gammaj).E(), pos, pid,calo->GetType(),0);
445 aodParticle->AddCluster(caloCluster);
451 else{ //Fill histograms
453 calo->GetMomentum(gammai,vertex);//Assume that come from vertex in straight line
456 if(pt < GetMinPt() || pt > GetMaxPt()) continue ;
462 if(GetDebug() > 2 ) printf("Neutral Hadron Correlation: Histograms selected gamma pair: pt %2.2f, phi %2.2f, eta %2.2f",pt,phi,eta);
464 fhEtaNeutral->Fill(ptTrig,eta);
465 fhPhiNeutral->Fill(ptTrig,phi);
466 fhDeltaEtaNeutral->Fill(ptTrig,etaTrig-eta);
467 fhDeltaPhiNeutral->Fill(ptTrig,phiTrig-phi);
468 fhDeltaPhiNeutralPt->Fill(pt,phiTrig-phi);
469 //Selection within angular range
470 if(((phiTrig-phi)> fDeltaPhiMinCut) && ((phiTrig-phi)<fDeltaPhiMaxCut) ){
471 if(GetDebug() > 2 ) printf("Selected neutral for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f ",pt,phi,eta);
472 fhPtImbalanceNeutral->Fill(ptTrig,rat);
479 //____________________________________________________________________________
480 Bool_t AliAnaParticleHadronCorrelation::SelectCluster(AliAODCaloCluster * calo, Double_t *vertex, TLorentzVector & mom, Int_t & pdg){
481 //Select cluster depending on its pid and acceptance selections
483 //Skip matched clusters with tracks
484 if(calo->GetNTracksMatched() > 0) {
488 calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
489 pdg = AliCaloPID::kPhoton;
491 //Get most probable PID, 2 options check PID weights (in MC this option is mandatory)
492 //or redo PID, recommended option for EMCal.
493 TString detector = "";
494 if(calo->IsPHOSCluster()) detector= "PHOS";
495 else if(calo->IsEMCALCluster()) detector= "EMCAL";
497 if(!IsCaloPIDRecalculationOn() || GetReader()->GetDataType() == AliCaloTrackReader::kMC )
498 pdg = GetCaloPID()->GetPdg(detector,calo->PID(),mom.E());//PID with weights
500 pdg = GetCaloPID()->GetPdg(detector,mom,calo->GetM02(),0,0,0,0);//PID recalculated
502 if(GetDebug() > 1) printf("PDG of identified particle %d\n",pdg);
504 //If it does not pass pid, skip
505 if(pdg != AliCaloPID::kPhoton || pdg != AliCaloPID::kPi0) {
510 //Check acceptance selection
511 if(IsFidutialCutOn()){
513 if(calo->IsPHOSCluster())
514 in = GetFidutialCut()->IsInFidutialCut(mom,"PHOS") ;
515 else if(calo->IsEMCALCluster())
516 in = GetFidutialCut()->IsInFidutialCut(mom,"EMCAL") ;
517 if(! in ) { return kFALSE ;}
520 if(GetDebug() > 1) printf("Correlation photon selection cuts passed: pT %3.2f, pdg %d\n",mom.Pt(), pdg);