AliPhysics  97344c9 (97344c9)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskEmcalJetShapesMC.cxx
Go to the documentation of this file.
1 //
2 // Jet QG tagging analysis task.
3 //
4 // Author: D. Caffarri, L. Cunqueiro
5 
6 #include <TClonesArray.h>
7 #include <TH1F.h>
8 #include <TH2F.h>
9 #include <TH3F.h>
10 #include <THnSparse.h>
11 #include <TTree.h>
12 #include <TList.h>
13 #include <TLorentzVector.h>
14 #include <TProfile.h>
15 #include <TChain.h>
16 #include <TSystem.h>
17 #include <TFile.h>
18 #include <TKey.h>
19 #include <AliAnalysisDataSlot.h>
20 #include <AliAnalysisDataContainer.h>
21 #include "TMatrixD.h"
22 #include "TMatrixDSym.h"
23 #include "TMatrixDSymEigen.h"
24 #include "TVector3.h"
25 #include "TVector2.h"
26 #include "AliVCluster.h"
27 #include "AliVTrack.h"
28 #include "AliEmcalJet.h"
29 #include "AliRhoParameter.h"
30 #include "AliLog.h"
31 #include "AliEmcalParticle.h"
32 #include "AliMCEvent.h"
33 #include "AliGenPythiaEventHeader.h"
34 #include "AliAODMCHeader.h"
35 #include "AliMCEvent.h"
36 #include "AliAnalysisManager.h"
37 #include "AliJetContainer.h"
38 #include "AliParticleContainer.h"
39 #include "AliEmcalPythiaInfo.h"
40 #include "TRandom3.h"
41 #include "AliEmcalJetFinder.h"
42 #include "AliAODEvent.h"
44 
45 using std::cout;
46 using std::endl;
47 
49 
50 //________________________________________________________________________
53  fContainer(0),
54  fMinFractionShared(0),
55  fJetShapeType(kGenShapes),
56  fJetShapeSub(kNoSub),
57  fJetSelection(kInclusive),
58  fPtThreshold(-9999.),
59  fRMatching(0.2),
60  fJetRadius(0.4),
61  fSubjetRadius(0.2),
62  fSelectedShapes(0),
63  fminpTTrig(20.),
64  fmaxpTTrig(50.),
65  fangWindowRecoil(0.6),
66  fSemigoodCorrect(0),
67  fHolePos(0),
68  fHoleWidth(0),
69  fCentSelectOn(kTRUE),
70  fCentMin(0),
71  fCentMax(10),
72  fOneConstSelectOn(kFALSE),
73  fDerivSubtrOrder(0),
74  fPhiJetCorr6(0x0),
75  fPhiJetCorr7(0x0),
76  fEtaJetCorr6(0x0),
77  fEtaJetCorr7(0x0),
78  fPtJetCorr(0x0),
79  fPtJet(0x0),
80  fhpTjetpT(0x0),
81  fhPt(0x0),
82  fhPhi(0x0),
83  fNbOfConstvspT(0x0),
84  fTreeObservableTagging(0x0)
85 
86 {
87  for(Int_t i=0;i<23;i++){
88  fShapesVar[i]=0;}
89  SetMakeGeneralHistograms(kTRUE);
90 }
91 
92 //________________________________________________________________________
94  AliAnalysisTaskEmcalJet(name, kTRUE),
95  fContainer(0),
96  fMinFractionShared(0),
97  fJetShapeType(kGenShapes),
98  fJetShapeSub(kNoSub),
99  fJetSelection(kInclusive),
100  fPtThreshold(-9999.),
101  fRMatching(0.2),
102  fSelectedShapes(0),
103  fminpTTrig(20.),
104  fmaxpTTrig(50.),
105  fangWindowRecoil(0.6),
106  fSemigoodCorrect(0),
107  fHolePos(0),
108  fHoleWidth(0),
109  fCentSelectOn(kTRUE),
110  fCentMin(0),
111  fCentMax(10),
112  fOneConstSelectOn(kFALSE),
113  fDerivSubtrOrder(0),
114  fPhiJetCorr6(0x0),
115  fPhiJetCorr7(0x0),
116  fEtaJetCorr6(0x0),
117  fEtaJetCorr7(0x0),
118  fPtJetCorr(0x0),
119  fPtJet(0x0),
120  fhpTjetpT(0x0),
121  fhPt(0x0),
122  fhPhi(0x0),
123  fNbOfConstvspT(0x0),
124  fTreeObservableTagging(0x0)
125 
126 {
127  // Standard constructor.
128 
129 
130  for(Int_t i=0;i<23;i++){
131  fShapesVar[i]=0;}
132 
134 
135  DefineOutput(1, TList::Class());
136  DefineOutput(2, TTree::Class());
137 
138 
139 }
140 
141 //________________________________________________________________________
143 {
145  delete fTreeObservableTagging;
147  }
148 
149 }
150 
151 //________________________________________________________________________
153 {
154  // Create user output.
156 
157  Bool_t oldStatus = TH1::AddDirectoryStatus();
158  TH1::AddDirectory(oldStatus);
159 
160  //fTreeObservableTagging = new TTree("fTreeJetShape", "fTreeJetShape");
161 
162  //TH1::AddDirectory(oldStatus);
163 
164  const char* nameoutput = GetOutputSlot(2)->GetContainer()->GetName();
165  fTreeObservableTagging = new TTree(nameoutput, nameoutput);
166 
167  const Int_t nVar = 23;
168 
169  TString *fShapesVarNames = new TString [nVar];
170 
171  fShapesVarNames[0] = "partonCode";
172  fShapesVarNames[1] = "ptJet";
173  fShapesVarNames[2] = "ptDJet";
174  fShapesVarNames[3] = "mJet";
175  fShapesVarNames[4] = "nbOfConst";
176  fShapesVarNames[5] = "angularity";
177  fShapesVarNames[6] = "circularity";
178  fShapesVarNames[7] = "lesub";
179  fShapesVarNames[8] = "CoreFraction";
180  fShapesVarNames[9] = "Nsubjet1";
181  fShapesVarNames[10] = "Nsubjet2";
182  fShapesVarNames[11] = "DeltaR";
183  fShapesVarNames[12] = "OpenAngle";
184  fShapesVarNames[13] = "weightPythia";
185 
186  fShapesVarNames[14] = "NT70";
187  fShapesVarNames[15] = "nConstNT70";
188  fShapesVarNames[16] = "NT80";
189  fShapesVarNames[17] = "nConstNT80";
190  fShapesVarNames[18] = "NT90";
191  fShapesVarNames[19] = "nConstNT90";
192  fShapesVarNames[20] = "NT95";
193  fShapesVarNames[21] = "nConstNT95";
194 
195  fShapesVarNames[22] = "SubjetFraction";
196 
197 
198  for(Int_t ivar=0; ivar < nVar; ivar++){
199  cout<<"looping over variables"<<endl;
200  fTreeObservableTagging->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/F", fShapesVarNames[ivar].Data()));}
201 
202 
203 
204  fPhiJetCorr6= new TH2F("fPhiJetCorr6", "fPhiJetCorr6", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
205  fOutput->Add(fPhiJetCorr6);
206  fEtaJetCorr6= new TH2F("fEtaJetCorr6", "fEtaJetCorr6", 50, -1.5, 1.5, 50, -1.5, 1.5);
207  fOutput->Add(fEtaJetCorr6);
208 
209  fPhiJetCorr7= new TH2F("fPhiJetCorr7", "fPhiJetCorr7", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
210  fOutput->Add(fPhiJetCorr7);
211  fEtaJetCorr7= new TH2F("fEtaJetCorr7", "fEtaJetCorr7", 50, -1.5, 1.5, 50, -1.5, 1.5);
212  fOutput->Add(fEtaJetCorr7);
213 
214  fPtJetCorr= new TH2F("fPtJetCorr", "fPtJetCorr", 100, 0, 200, 100, 0, 200);
215  fOutput->Add(fPtJetCorr);
216  fPtJet= new TH1F("fPtJet", "fPtJet", 100, 0, 200);
217  fOutput->Add(fPtJet);
218 
219  fhpTjetpT= new TH2F("fhpTjetpT", "fhpTjetpT", 200, 0, 200, 200, 0, 200);
220  fOutput->Add(fhpTjetpT);
221  fhPt= new TH1F("fhPt", "fhPt", 200, 0, 200);
222  fOutput->Add(fhPt);
223  fhPhi= new TH1F("fhPhi", "fhPhi", 100, -TMath::Pi(), TMath::Pi());
224  fOutput->Add(fhPhi);
225 
226  fNbOfConstvspT=new TH2F("fNbOfConstvspT", "fNbOfConstvspT", 100, 0, 100, 200, 0, 200);
227  fOutput->Add(fNbOfConstvspT);
228 
229  //fOutput->Add(fTreeObservableTagging);
230 
231 
232  PostData(1, fOutput); // Post data for ALL output slots > 0 here
233  PostData(2, fTreeObservableTagging);
234 
235  delete [] fShapesVarNames;
236 
237 }
238 
239 //________________________________________________________________________
241 {
242  // Run analysis code here, if needed. It will be executed before FillHistograms().
243 
244  return kTRUE;
245 }
246 
247 //________________________________________________________________________
249 {
250  // Fill histograms.
251  //cout<<"IntoFillHistograms"<<endl;
252  AliEmcalJet* jet1 = NULL;
253  AliJetContainer *jetCont = GetJetContainer(0);
254 
255  Float_t kWeight=1;
256  if (fCentSelectOn)
257  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
258 
259  AliAODTrack *triggerHadron = 0x0;
260 
261  if (fJetSelection == kRecoil) {
262  //Printf("Recoil jets!!!, fminpTTrig = %f, fmaxpTTrig = %f", fminpTTrig, fmaxpTTrig);
263  Int_t triggerHadronLabel = SelectTrigger(fminpTTrig, fmaxpTTrig);
264 
265 
266  if (triggerHadronLabel==-99999) {
267  //Printf ("Trigger Hadron not found, return");
268  return 0;}
269 
270 
272  TClonesArray *trackArrayAn = partContAn->GetArray();
273  triggerHadron = static_cast<AliAODTrack*>(trackArrayAn->At(triggerHadronLabel));
274 
275  if (!triggerHadron) {
276  //Printf("No Trigger hadron with the found label!!");
277  return 0;
278  }
279 
280  if(fSemigoodCorrect){
281  Double_t disthole=RelativePhi(triggerHadron->Phi(),fHolePos);
282  if(TMath::Abs(disthole)+fHoleWidth>TMath::Pi()-fangWindowRecoil){
283  return 0;}
284  }
285 
286  fhPt->Fill(triggerHadron->Pt());
287 
288  }
289 
290  if(jetCont) {
291  jetCont->ResetCurrentID();
292  while((jet1 = jetCont->GetNextAcceptJet())) {
293  //Printf("jet1=%p", jet1);
294  if (!jet1) continue;
295  AliEmcalJet* jet2 = 0x0;
296  AliEmcalJet* jet3 = 0x0;
297  fPtJet->Fill(jet1->Pt());
298  AliEmcalJet *jetUS = NULL;
299  Int_t ifound=0, jfound=0;
300  Int_t ilab=-1, jlab=-1;
301 
303  Double_t disthole=RelativePhi(jet1->Phi(),fHolePos);
304  if(TMath::Abs(disthole)<fHoleWidth){
305  continue;
306  }
307  }
308 
309  Float_t dphiRecoil = 0.;
310  if (fJetSelection == kRecoil){
311  dphiRecoil = RelativePhi(triggerHadron->Phi(), jet1->Phi());
312  if (TMath::Abs(dphiRecoil) < (TMath::Pi() - fangWindowRecoil)) {
313  // Printf("Recoil jets back to back not found! continuing");
314  continue;
315  }
316 
317  fhpTjetpT->Fill(triggerHadron->Pt(), jet1->Pt());
318  //Printf(" ************ FILLING HISTOS****** shapeSub = %d, triggerHadron = %f, jet1 = %f", fJetShapeSub, triggerHadron->Pt(), jet1->Pt());
319  fhPhi->Fill(RelativePhi(triggerHadron->Phi(), jet1->Phi()));
320 
321  }
322 
323 
324  fShapesVar[0] = 0.;
325 
326  if (fJetShapeType == kGenShapes){
327  const AliEmcalPythiaInfo *partonsInfo = 0x0;
328  partonsInfo = GetPythiaInfo();
329  //Printf("partonsInfo=%p", partonsInfo);
330  Double_t jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi6());
331  Double_t detap1=(jet1->Eta())-(partonsInfo->GetPartonEta6());
332  kWeight=partonsInfo->GetPythiaEventWeight();
333  //Printf("kWeight=%f", kWeight);
334  fShapesVar[13] = kWeight;
335 
336  Float_t dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
337  fEtaJetCorr6->Fill(jet1->Eta(), partonsInfo->GetPartonEta6());
338  fPhiJetCorr6->Fill(jet1->Phi(), partonsInfo->GetPartonPhi6());
339  if(dRp1 < fRMatching) {
340  fShapesVar[0] = partonsInfo->GetPartonFlag6();
341  fPtJetCorr ->Fill(partonsInfo->GetPartonPt6(), jet1->Pt());
342  }
343  else {
344  jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi7());
345  detap1=(jet1->Eta())-(partonsInfo->GetPartonEta7());
346  dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
347  fEtaJetCorr7->Fill(jet1->Eta(), partonsInfo->GetPartonEta7());
348  fPhiJetCorr7->Fill(jet1->Phi(), partonsInfo->GetPartonPhi7());
349  if(dRp1 < fRMatching) {
350  fShapesVar[0] = partonsInfo->GetPartonFlag7();
351  fPtJetCorr->Fill(partonsInfo->GetPartonPt7(), jet1->Pt());
352  }
353  else fShapesVar[0]=0;
354  }
355  }
356 
357  Double_t ptSubtracted = 0;
358  ptSubtracted= jet1->Pt();
359  //Printf("ptSubtracted=%f", ptSubtracted);
360 
361 
362  if (ptSubtracted < fPtThreshold) continue;
363 
364  if ((fCentSelectOn == kFALSE) && (jet1->GetNumberOfTracks() <= 1)) continue;
365 
366  AliEmcalJetFinder *Reclusterer1; //Object containg Subjets from Subtracted Hybrid Jets
367  Reclusterer1 = Recluster(jet1, 0, fJetRadius, fSubjetRadius, 1, 0, "SubJetFinder_1");
368 
369  fShapesVar[1] = ptSubtracted;
370  fShapesVar[2] = GetJetpTD(jet1,0);
371  fShapesVar[3] = GetJetMass(jet1,0);
372  fShapesVar[4] = 1.*GetJetNumberOfConstituents(jet1,0);
373  fShapesVar[5] = GetJetAngularity(jet1,0);
374  fShapesVar[6] = GetJetCircularity(jet1,0);
375  fShapesVar[7] = GetJetLeSub(jet1,0);
376  fShapesVar[8] = GetJetCoreFrac(jet1,0);
377  fShapesVar[9] = fjNSubJettiness(jet1,0,1,0,1,0);
378  fShapesVar[10]= fjNSubJettiness(jet1,0,2,0,1,0);
379  fShapesVar[11]= fjNSubJettiness(jet1,0,2,0,1,2);
380  fShapesVar[12]= fjNSubJettiness(jet1,0,2,0,1,1);
381 
382  Float_t nTFractions[8]={0.,0.,0.,0.,0.,0.,0.,0.};
383  NTValues(jet1, 0, nTFractions);
384  //shape 13 is pythia weight!
385  for (Int_t ishape=14; ishape<22; ishape++) fShapesVar[ishape] = nTFractions[ishape-14];
386 
387  fShapesVar[22]= GetSubjetFraction(jet1,0,fJetRadius,Reclusterer1);
388 
389  fTreeObservableTagging->Fill();
390 
391  }
392 
393  }
394 
395  return kTRUE;
396 }
397 
398 //________________________________________________________________________
400  //calc subtracted jet mass
401  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
403  else return jet->GetShapeProperties()->GetSecondOrderSubtracted();
404  else
405  return jet->M();
406 }
407 
408 //________________________________________________________________________
410 
411  AliJetContainer *jetCont = GetJetContainer(jetContNb);
412  if (!jet->GetNumberOfTracks())
413  return 0;
414  Double_t den=0.;
415  Double_t num = 0.;
416  AliVParticle *vp1 = 0x0;
417  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
418  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
419 
420  if (!vp1){
421  Printf("AliVParticle associated to constituent not found");
422  continue;
423  }
424 
425  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
426  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
427  Double_t dr = TMath::Sqrt(dr2);
428  num=num+vp1->Pt()*dr;
429  den=den+vp1->Pt();
430  }
431  return num/den;
432 }
433 
434 //________________________________________________________________________
436 
437  if((fJetShapeSub==kDerivSub) && (jetContNb==0))
440  else
441  return Angularity(jet, jetContNb);
442 
443 }
444 
445 
446 //________________________________________________________________________
448 
449  AliJetContainer *jetCont = GetJetContainer(jetContNb);
450  if (!jet->GetNumberOfTracks())
451  return 0;
452  Double_t den=0.;
453  Double_t num = 0.;
454  AliVParticle *vp1 = 0x0;
455  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
456  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
457 
458  if (!vp1){
459  Printf("AliVParticle associated to constituent not found");
460  continue;
461  }
462 
463  num=num+vp1->Pt()*vp1->Pt();
464  den=den+vp1->Pt();
465  }
466  return TMath::Sqrt(num)/den;
467 }
468 
469 //________________________________________________________________________
471  //calc subtracted jet mass
472  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
474  else return jet->GetShapeProperties()->GetSecondOrderSubtractedpTD();
475  else
476  return PTD(jet, jetContNb);
477 
478 }
479 
480 //_____________________________________________________________________________
482 
483  AliJetContainer *jetCont = GetJetContainer(jetContNb);
484  if (!jet->GetNumberOfTracks())
485  return 0;
486  Double_t mxx = 0.;
487  Double_t myy = 0.;
488  Double_t mxy = 0.;
489  int nc = 0;
490  Double_t sump2 = 0.;
491  Double_t pxjet=jet->Px();
492  Double_t pyjet=jet->Py();
493  Double_t pzjet=jet->Pz();
494 
495 
496  //2 general normalized vectors perpendicular to the jet
497  TVector3 ppJ1(pxjet, pyjet, pzjet);
498  TVector3 ppJ3(- pxjet* pzjet, - pyjet * pzjet, pxjet * pxjet + pyjet * pyjet);
499  ppJ3.SetMag(1.);
500  TVector3 ppJ2(-pyjet, pxjet, 0);
501  ppJ2.SetMag(1.);
502  AliVParticle *vp1 = 0x0;
503  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
504  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
505 
506  if (!vp1){
507  Printf("AliVParticle associated to constituent not found");
508  continue;
509  }
510 
511  TVector3 pp(vp1->Px(), vp1->Py(), vp1->Pz());
512 
513  //local frame
514  TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
515  TVector3 pPerp = pp - pLong;
516  //projection onto the two perpendicular vectors defined above
517 
518  Float_t ppjX = pPerp.Dot(ppJ2);
519  Float_t ppjY = pPerp.Dot(ppJ3);
520  Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
521  if(ppjT<=0) return 0;
522 
523  mxx += (ppjX * ppjX / ppjT);
524  myy += (ppjY * ppjY / ppjT);
525  mxy += (ppjX * ppjY / ppjT);
526  nc++;
527  sump2 += ppjT;}
528 
529  if(nc<2) return 0;
530  if(sump2==0) return 0;
531  // Sphericity Matrix
532  Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
533  TMatrixDSym m0(2,ele);
534 
535  // Find eigenvectors
536  TMatrixDSymEigen m(m0);
537  TVectorD eval(2);
538  TMatrixD evecm = m.GetEigenVectors();
539  eval = m.GetEigenValues();
540  // Largest eigenvector
541  int jev = 0;
542  // cout<<eval[0]<<" "<<eval[1]<<endl;
543  if (eval[0] < eval[1]) jev = 1;
544  TVectorD evec0(2);
545  // Principle axis
546  evec0 = TMatrixDColumn(evecm, jev);
547  Double_t compx=evec0[0];
548  Double_t compy=evec0[1];
549  TVector2 evec(compx, compy);
550  Double_t circ=0;
551  if(jev==1) circ=2*eval[0];
552  if(jev==0) circ=2*eval[1];
553 
554  return circ;
555 
556 
557 
558 }
559 
560 
561 
562 
563 //________________________________________________________________________
565  //calc subtracted jet mass
566 
567  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
570  else
571  return Circularity(jet, jetContNb);
572 
573 }
574 
575 //________________________________________________________________________
577 
578  AliJetContainer *jetCont = GetJetContainer(jetContNb);
579  if (!jet->GetNumberOfTracks())
580  return 0;
581  Double_t den=0.;
582  Double_t num = 0.;
583  AliVParticle *vp1 = 0x0;
584  AliVParticle *vp2 = 0x0;
585  std::vector<int> ordindex;
586  ordindex=jet->GetPtSortedTrackConstituentIndexes(jetCont->GetParticleContainer()->GetArray());
587  //Printf("Nbof const = %d", jet->GetNumberOfTracks());
588  //Printf("ordindex[0] = %d, ordindex[1] = %d", ordindex[0], ordindex[1]);
589 
590  if(ordindex.size()<2) return -1;
591 
592  vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[0], jetCont->GetParticleContainer()->GetArray()));
593  if (!vp1){
594  Printf("AliVParticle associated to Leading constituent not found");
595  return -1;
596  }
597 
598  vp2 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[1], jetCont->GetParticleContainer()->GetArray()));
599  if (!vp2){
600  Printf("AliVParticle associated to Subleading constituent not found");
601  return -1;
602  }
603 
604 
605  num=vp1->Pt();
606  den=vp2->Pt();
607  //Printf("vp1->Pt() =%f, vp2->Pt() =%f", vp1->Pt(), vp2->Pt());
608 
609 return num-den;
610 }
611 
612 //________________________________________________________________________
614  //calc subtracted jet mass
615 
616  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
619  else
620  return LeSub(jet, jetContNb);
621 
622 }
623 
624 //________________________________________________________________________
626  //calc subtracted jet mass
627 
628  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
631  else
632  return jet->GetNumberOfTracks();
633 
634  }
635 
636 
637 //________________________________________________________________________
638 AliEmcalJetFinder *AliAnalysisTaskEmcalJetShapesMC::Recluster(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, Double_t SubJetRadius, Double_t SubJetMinPt, Int_t Algorithm, const char* Name){
639 
640  AliJetContainer *JetCont = GetJetContainer(JetContNb);
641  AliEmcalJetFinder *Reclusterer = new AliEmcalJetFinder(Name); //JetFinder Object for reclustered jets
642  Reclusterer->SetRadius(SubJetRadius);
643  Reclusterer->SetJetMinPt(SubJetMinPt);
644  Reclusterer->SetJetAlgorithm(Algorithm); //0 for anti-kt 1 for kt
645  Reclusterer->SetJetMaxEta(0.9-JetRadius);
646  Reclusterer->SetRecombSheme(0);
647  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
648 
649  //Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
650  Double_t dVtx[3]={0.,0.,0.};
651  if(Reclusterer->AliEmcalJetFinder::Filter(Jet, JetCont, dVtx)){;} //reclustering jet1 using the jetfinderobject Reclusterer
652  return Reclusterer;
653 }
654 
655 
656 
657 //________________________________________________________________________
658 /*Double_t AliAnalysisTaskEmcalJetShapesMC::NSubJettiness(AliEmcalJet *Jet, Int_t JetContNb, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t A, Int_t B){
659  AliJetContainer *JetCont = GetJetContainer(JetContNb);
660  AliEmcalJet *SubJet=NULL;
661  Double_t DeltaR1=0;
662  Double_t DeltaR2=0;
663  AliVParticle *JetParticle=0x0;
664  Double_t SubJetiness_Numerator = 0;
665  Double_t SubJetiness_Denominator = 0;
666  Double_t Index=-2;
667  Bool_t Error=kFALSE;
668  if (!Jet->GetNumberOfTracks()) return -2;
669  if (Reclusterer->GetNumberOfJets() < N) return -2;
670  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
671  JetParticle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
672  for (Int_t j=1; j<=N; j++){
673  Index=SubJetOrdering(Jet,Reclusterer,j,0,kTRUE);
674  if(Index==-999){
675  Error = kTRUE;
676  i=Jet->GetNumberOfTracks();
677  break;
678  }
679  if(j==1){
680  DeltaR1=TMath::Power((Reclusterer->GetJet(Index)->Pt()),A)*TMath::Power((TMath::Sqrt((((JetParticle->Eta())-(Reclusterer->GetJet(Index)->Eta()))*((JetParticle->Eta())- (Reclusterer->GetJet(Index)->Eta())))+((RelativePhi((Reclusterer->GetJet(Index)->Phi()),JetParticle->Phi()))*(RelativePhi((Reclusterer->GetJet(Index)->Phi()),JetParticle->Phi()))))),B);
681  }
682  else{
683  DeltaR2=TMath::Power((Reclusterer->GetJet(Index)->Pt()),A)*TMath::Power((TMath::Sqrt((((JetParticle->Eta())-(Reclusterer->GetJet(Index)->Eta()))*((JetParticle->Eta())- (Reclusterer->GetJet(Index)->Eta())))+((RelativePhi((Reclusterer->GetJet(Index)->Phi()),JetParticle->Phi()))*(RelativePhi((Reclusterer->GetJet(Index)->Phi()),JetParticle->Phi()))))),B);
684  if (DeltaR2<DeltaR1) DeltaR1=DeltaR2;
685  }
686  }
687  SubJetiness_Numerator=SubJetiness_Numerator+(JetParticle->Pt()*DeltaR1);
688  if (A>=0) SubJetiness_Denominator=SubJetiness_Denominator+(TMath::Power((Reclusterer->GetJet(SubJetOrdering(Jet,Reclusterer,1,0,kTRUE))->Pt()),A)*JetParticle->Pt()*TMath::Power(fJetRadius,B));
689  else SubJetiness_Denominator=SubJetiness_Denominator+(TMath::Power((Reclusterer->GetJet(SubJetOrdering(Jet,Reclusterer,N,0,kTRUE))->Pt()),A)*JetParticle->Pt()*TMath::Power(fJetRadius,B));
690  }
691  if (SubJetiness_Denominator!=0 && !Error) return SubJetiness_Numerator/SubJetiness_Denominator;
692  else return -2;
693 
694 }
695 */
696 
697 //----------------------------------------------------------------------
699  AliJetContainer *JetCont = GetJetContainer(JetContNb);
700  AliEmcalJet *SubJet=NULL;
701  Double_t DeltaR1=0;
702  Double_t DeltaR2=0;
703  AliVParticle *JetParticle=0x0;
704  Double_t SubJetiness_Numerator = 0;
705  Double_t SubJetiness_Denominator = 0;
706  Double_t Index=-2;
707  if (Reclusterer->GetNumberOfJets() < 1) return -2;
708  Index=SubJetOrdering(Jet,Reclusterer,1,0,kTRUE);
709  if(Index==-999) return -2;
710  SubJetiness_Numerator=(Reclusterer->GetJet(Index)->Pt());
711  SubJetiness_Denominator=Jet->Pt();
712  return SubJetiness_Numerator/SubJetiness_Denominator;
713 
714 
715 }
716 //__________________________________________________________________________________
718 
719  AliJetContainer *jetCont = GetJetContainer(jetContNb);
720  if (!jet->GetNumberOfTracks())
721  return 0;
722  Double_t den=0.;
723  Double_t num = 0.;
724  AliVParticle *vp1 = 0x0;
725  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
726  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
727 
728  if (!vp1){
729  Printf("AliVParticle associated to constituent not found");
730  continue;
731  }
732 
733  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
734  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
735  Double_t dr = TMath::Sqrt(dr2);
736  if(dr<=fSubjetRadius) num=num+vp1->Pt();
737 
738  }
739  return num/jet->Pt();
740 }
741 
742 
743 
744 
745 //________________________________________________________________________
747  //calc subtracted jet mass
748 
749  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
752  else
753  return CoreFrac(jet, jetContNb);
754 
755 }
756 
757 
758 
759 
760 //----------------------------------------------------------------------
762  AliEmcalJet *SubJet=NULL;
763  Double_t SortingVariable;
764  Int_t ArraySize =N+1;
765  TArrayD *JetSorter = new TArrayD(ArraySize);
766  TArrayD *JetIndexSorter = new TArrayD(ArraySize);
767  for (Int_t i=0; i<ArraySize; i++){
768  JetSorter->SetAt(0,i);
769  }
770  for (Int_t i=0; i<ArraySize; i++){
771  JetIndexSorter->SetAt(0,i);
772  }
773  if(Reclusterer->GetNumberOfJets()<N) return -999;
774  for (Int_t i=0; i<Reclusterer->GetNumberOfJets(); i++){
775  SubJet=Reclusterer->GetJet(i);
776  if (Type==0) SortingVariable=SubJet->Pt();
777  else if (Type==1) SortingVariable=SubJet->E();
778  else if (Type==2) SortingVariable=SubJet->M();
779  for (Int_t j=0; j<N; j++){
780  if (SortingVariable>JetSorter->GetAt(j)){
781  for (Int_t k=N-1; k>=j; k--){
782  JetSorter->SetAt(JetSorter->GetAt(k),k+1);
783  JetIndexSorter->SetAt(JetIndexSorter->GetAt(k),k+1);
784  }
785  JetSorter->SetAt(SortingVariable,j);
786  JetIndexSorter->SetAt(i,j);
787  break;
788  }
789  }
790  }
791  if (!Index) return JetSorter->GetAt(N-1);
792  else return JetIndexSorter->GetAt(N-1);
793 }
794 
795 
796 
797 //returns -1 if the Nth hardest jet is requested where N>number of available jets
798 //type: 0=Pt 1=E 2=M
799 //Index TRUE=returns index FALSE=returns value of quantatiy in question
800 
801 
802 
803 
804 
805 //______________________________________________________________________________
807 
808  AliJetContainer *jetCont = GetJetContainer(jetContNb);
809  if (!jet->GetNumberOfTracks())
810  return 0;
811  Double_t mxx = 0.;
812  Double_t myy = 0.;
813  Double_t mxy = 0.;
814  int nc = 0;
815  Double_t sump2 = 0.;
816 
817  AliVParticle *vp1 = 0x0;
818  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
819  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
820 
821  if (!vp1){
822  Printf("AliVParticle associated to constituent not found");
823  continue;
824  }
825 
826  Double_t ppt=vp1->Pt();
827  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
828 
829  Double_t deta = vp1->Eta()-jet->Eta();
830  mxx += ppt*ppt*deta*deta;
831  myy += ppt*ppt*dphi*dphi;
832  mxy -= ppt*ppt*deta*TMath::Abs(dphi);
833  nc++;
834  sump2 += ppt*ppt;
835 
836  }
837  if(nc<2) return 0;
838  if(sump2==0) return 0;
839  // Sphericity Matrix
840  Double_t ele[4] = {mxx , mxy , mxy , myy };
841  TMatrixDSym m0(2,ele);
842 
843  // Find eigenvectors
844  TMatrixDSymEigen m(m0);
845  TVectorD eval(2);
846  TMatrixD evecm = m.GetEigenVectors();
847  eval = m.GetEigenValues();
848  // Largest eigenvector
849  int jev = 0;
850  // cout<<eval[0]<<" "<<eval[1]<<endl;
851  if (eval[0] < eval[1]) jev = 1;
852  TVectorD evec0(2);
853  // Principle axis
854  evec0 = TMatrixDColumn(evecm, jev);
855  Double_t compx=evec0[0];
856  Double_t compy=evec0[1];
857  TVector2 evec(compx, compy);
858  Double_t sig=0;
859  if(jev==1) sig=TMath::Sqrt(TMath::Abs(eval[0])/sump2);
860  if(jev==0) sig=TMath::Sqrt(TMath::Abs(eval[1])/sump2);
861 
862  return sig;
863 
864 }
865 
866 //________________________________________________________________________
868  //calc subtracted jet mass
869 
870  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
873  else
874  return Sigma2(jet, jetContNb);
875 
876 }
877 
878 
879 //_________________________________________________________________________
881 
882  AliJetContainer *jetCont = GetJetContainer(jetContNb);
883  if (!jet->GetNumberOfTracks())
884  return;
885 
886  Double_t ptJet = jet->Pt();
887 
888  AliVParticle *vp1 = 0x0;
889  std::vector<int> ordindex;
890  ordindex=jet->GetPtSortedTrackConstituentIndexes(jetCont->GetParticleContainer()->GetArray());
891  //Printf("Nbof const = %d", jet->GetNumberOfTracks());
892  //Printf("ordindex[0] = %d, ordindex[1] = %d", ordindex[0], ordindex[1]);
893  //if(ordindex.size()<2) return -1;
894 
895  for(Int_t iconst =0; iconst<jet->GetNumberOfTracks(); iconst++){
896 
897  vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[iconst], jetCont->GetParticleContainer()->GetArray()));
898  if (!vp1){
899  Printf("AliVParticle associated to Leading constituent not found");
900  return;
901  }
902 
903  if (nTFractions[0] <= 0.7*ptJet){
904  nTFractions[0] += vp1->Pt();
905  nTFractions[1] +=1;
906  }
907 
908  if (nTFractions[2] <= 0.8*ptJet){
909  nTFractions[2] += vp1->Pt();
910  nTFractions[3] +=1;
911  }
912 
913  if (nTFractions[4] <= 0.9*ptJet){
914  nTFractions[4] += vp1->Pt();
915  nTFractions[5] +=1;
916  }
917 
918  if (nTFractions[6] <= 0.95*ptJet){
919  nTFractions[6] += vp1->Pt();
920  nTFractions[7] +=1;
921  }
922  }
923 }
924 //_________________________________________________________________________________________________
926 
927  //WARNING!!! Only works for parent jets that are clustered with Anti-Kt! To change go to AliEmcalJetFinder.cxx and look at the Nsubjettiness() function
928 
929  //Algorithm==0 -> kt_axes;
930  // Algorithm==1 -> ca_axes;
931  //Algorithm==2 -> antikt_0p2_axes;
932  //Algorithm==3 -> wta_kt_axes;
933  //Algorithm==4 -> wta_ca_axes;
934  //Algorithm==5 -> onepass_kt_axes;
935  //Algorithm==6 -> onepass_ca_axes;
936  //Algorithm==7 -> onepass_antikt_0p2_axes;
937  //Algorithm==8 -> onepass_wta_kt_axes;
938  //Algorithm==9 -> onepass_wta_ca_axes;
939  //Algorithm==10 -> min_axes;
940 
941 
942  //Option==0 returns Nsubjettiness Value
943  //Option==1 && N==2 returns opening angle between two subjet axes(Delta R?)
944  //Option==2 && N==2 returns Delta R
945 
946  if (Jet->GetNumberOfTracks()>=N){
947  AliJetContainer *JetCont = GetJetContainer(JetContNb);
948  AliEmcalJetFinder *JetFinder=new AliEmcalJetFinder("Nsubjettiness");
949  //Printf("jetfinder =%p, JetCont=%p", JetFinder, JetCont);
950  JetFinder->SetJetMaxEta(0.9-fJetRadius);
951  JetFinder->SetRadius(fJetRadius);
952  JetFinder->SetJetAlgorithm(Algorithm); //0 for anti-kt 1 for kt //this is for the JET!!!!!!!!!! Not the SubJets
953  JetFinder->SetRecombSheme(0);
954  JetFinder->SetJetMinPt(Jet->Pt());
955  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
956  //Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
957  Double_t dVtx[3]={0,0,0};
958  //Printf("JetFinder->Nsubjettiness =%f", JetFinder->Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubjetRadius,Beta,Option));
959  return JetFinder->Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubjetRadius,Beta,Option);
960 
961  }
962  else return -2;
963 }
964 
965 
966 
967 //________________________________________________________________________
969 
971  TClonesArray *tracksArray = partCont->GetArray();
972 
973  if(!partCont || !tracksArray) return -99999;
974  AliAODTrack *track = 0x0;
975  AliEmcalParticle *emcPart = 0x0;
976 
977 
978  TList *trackList = new TList();
979  Int_t triggers[100];
980  for (Int_t iTrigger=0; iTrigger<100; iTrigger++) triggers[iTrigger] = 0;
981  Int_t iTT = 0;
982 
983  for(Int_t iTrack=0; iTrack <= tracksArray->GetEntriesFast(); iTrack++){
984 
985 
986  if (fJetShapeSub == kConstSub){
987  emcPart = static_cast<AliEmcalParticle*>(tracksArray->At(iTrack));
988  if (!emcPart) continue;
989  if(TMath::Abs(emcPart->Eta())>0.9) continue;
990  if (emcPart->Pt()<0.15) continue;
991 
992  if ((emcPart->Pt() >= minpT) && (emcPart->Pt()< maxpT)) {
993  trackList->Add(emcPart);
994  triggers[iTT] = iTrack;
995  iTT++;
996  }
997  }
998  else{
999  track = static_cast<AliAODTrack*>(tracksArray->At(iTrack));
1000  if (!track) continue;
1001  if(TMath::Abs(track->Eta())>0.9) continue;
1002  if (track->Pt()<0.15) continue;
1003  if (!(track->TestFilterBit(768))) continue;
1004 
1005  if ((track->Pt() >= minpT) && (track->Pt()< maxpT)) {
1006  trackList->Add(track);
1007  triggers[iTT] = iTrack;
1008  iTT++;
1009 
1010  }
1011  }
1012  }
1013 
1014  if (iTT == 0) return -99999;
1015  Int_t nbRn = 0, index = 0 ;
1016  TRandom3* random = new TRandom3(0);
1017  nbRn = random->Integer(iTT);
1018 
1019  index = triggers[nbRn];
1020  //Printf("iTT Total= %d, nbRn = %d, Index = %d",iTT, nbRn, index );
1021  return index;
1022 
1023 }
1024 
1025 //__________________________________________________________________________________
1027 
1028  if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1029  else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1030  if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1031  else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1032  double dphi = mphi-vphi;
1033  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1034  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1035  return dphi;//dphi in [-Pi, Pi]
1036 }
1037 
1038 
1039 //________________________________________________________________________
1041  //
1042  // retrieve event objects
1043  //
1045  return kFALSE;
1046 
1047  return kTRUE;
1048 }
1049 
1050 //_______________________________________________________________________
1052 {
1053  // Called once at the end of the analysis.
1054 
1055  AliInfo("Terminate");
1056  AliAnalysisTaskSE::Terminate();
1057 
1058  fOutput = dynamic_cast<AliEmcalList*> (GetOutputData(1));
1059  if (!fOutput) {
1060  AliError("fOutput not available");
1061  return;
1062  }
1063 
1064  fTreeObservableTagging = dynamic_cast<TTree*>(GetOutputData(2));
1065  if (!fTreeObservableTagging){
1066  Printf("ERROR: fTreeObservableTagging not available");
1067  return;
1068  }
1069 
1070 }
1071 
void SetRadius(Double_t val)
Float_t GetSigma2(AliEmcalJet *jet, Int_t jetContNb=0)
Double_t GetFirstOrderSubtractedAngularity() const
Float_t GetJetMass(AliEmcalJet *jet, Int_t jetContNb=0)
double Double_t
Definition: External.C:58
void SetJetMinPt(Double_t val)
Double_t GetSecondOrderSubtractedSigma2() const
Definition: External.C:236
Float_t GetPartonEta6() const
AliJetContainer * GetJetContainer(Int_t i=0) const
Float_t GetPartonEta7() const
Double_t GetSecondOrderSubtractedConstituent() const
Double_t Eta() const
Definition: AliEmcalJet.h:114
ClassImp(AliAnalysisTaskEmcalJetShapesMC) AliAnalysisTaskEmcalJetShapesMC
Double_t Py() const
Definition: AliEmcalJet.h:100
Double_t Phi() const
Definition: AliEmcalJet.h:110
Float_t GetPythiaEventWeight() const
Float_t GetPartonPhi7() const
Double_t Nsubjettiness(AliEmcalJet *pJet, AliJetContainer *pContJets, Double_t dVtx[3], Int_t N, Int_t Algorithm, Double_t Radius, Double_t Beta, Int_t Option=0, Int_t Measure=0, Double_t Beta_SD=0, Double_t ZCut=0.1)
Float_t CoreFrac(AliEmcalJet *jet, Int_t jetContNb=0)
Float_t GetJetCircularity(AliEmcalJet *jet, Int_t jetContNb=0)
Double_t E() const
Definition: AliEmcalJet.h:112
Float_t GetJetpTD(AliEmcalJet *jet, Int_t jetContNb=0)
Int_t GetPartonFlag7() const
Container for particles within the EMCAL framework.
Float_t GetPartonPhi6() const
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:132
Double_t Px() const
Definition: AliEmcalJet.h:99
Float_t GetJetCoreFrac(AliEmcalJet *jet, Int_t jetContNb=0)
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Double_t fjNSubJettiness(AliEmcalJet *Jet, Int_t JetContNb, Int_t N, Int_t Algorithm, Double_t Beta, Int_t Option)
Int_t GetPartonFlag6() const
AliParticleContainer * GetParticleContainer() const
Double_t GetFirstOrderSubtractedConstituent() const
void SetRecombSheme(Int_t val)
void SetJetAlgorithm(Int_t val)
int Int_t
Definition: External.C:63
Double_t Eta() const
Double_t Pt() const
unsigned int UInt_t
Definition: External.C:33
Float_t GetJetNumberOfConstituents(AliEmcalJet *jet, Int_t jetContNb=0)
float Float_t
Definition: External.C:68
Task to store and correlate the MC shapes.
TTree * fTreeObservableTagging
! Tree with tagging variables subtracted MC or true MC or raw
Double_t RelativePhi(Double_t mphi, Double_t vphi)
Double_t GetSecondOrderSubtractedAngularity() const
Double_t GetSubjetFraction(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, AliEmcalJetFinder *Reclusterer)
Definition: Option.C:68
AliEmcalJetFinder * Recluster(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, Double_t SubJetRadius, Double_t SubJetMinPt, Int_t Algorithm, const char *Name)
Double_t fCent
!event centrality
Float_t PTD(AliEmcalJet *jet, Int_t jetContNb=0)
void SetJetMaxEta(Double_t val)
AliEmcalJet * GetNextAcceptJet()
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Float_t GetJetAngularity(AliEmcalJet *jet, Int_t jetContNb=0)
Int_t SelectTrigger(Float_t minpT, Float_t maxpT)
AliEmcalJet * GetJet(Int_t index)
Double_t Pt() const
Definition: AliEmcalJet.h:102
Enhanced TList-derived class that implements correct merging for pt_hard binned production.
Definition: AliEmcalList.h:25
Float_t Angularity(AliEmcalJet *jet, Int_t jetContNb=0)
Double_t GetFirstOrderSubtractedCircularity() const
AliEmcalList * fOutput
!output list
Double_t SubJetOrdering(AliEmcalJet *Jet, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t Type, Bool_t Index)
Short_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:153
void NTValues(AliEmcalJet *jet, Int_t jetContNb, Float_t *nTFractions)
const Int_t nVar
Double_t GetSecondOrderSubtractedCircularity() const
Float_t GetJetLeSub(AliEmcalJet *jet, Int_t jetContNb=0)
Store some informaion about a Pythia eventThis class is used to store some information about a Pythia...
void SetMakeGeneralHistograms(Bool_t g)
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
const AliEmcalPythiaInfo * GetPythiaInfo() const
Double_t Pz() const
Definition: AliEmcalJet.h:101
Float_t Circularity(AliEmcalJet *jet, Int_t jetContNb=0)
Declaration of class AliEmcalPythiaInfo.
const char Option_t
Definition: External.C:48
Double_t GetFirstOrderSubtractedSigma2() const
Float_t Sigma2(AliEmcalJet *jet, Int_t jetContNb=0)
bool Bool_t
Definition: External.C:53
AliEmcalJetShapeProperties * GetShapeProperties() const
Definition: AliEmcalJet.h:253
Float_t GetPartonPt7() const
Float_t LeSub(AliEmcalJet *jet, Int_t jetContNb=0)
Double_t M() const
Definition: AliEmcalJet.h:113
Container for jet within the EMCAL jet framework.
std::vector< int > GetPtSortedTrackConstituentIndexes(TClonesArray *tracks) const
Float_t GetPartonPt6() const