AliPhysics  8bb951a (8bb951a)
 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  fShapesVar(0),
59  fPtThreshold(-9999.),
60  fRMatching(0.2),
61  fJetRadius(0.4),
62  fSubjetRadius(0.2),
63  fSelectedShapes(0),
64  fminpTTrig(20.),
65  fmaxpTTrig(50.),
66  fangWindowRecoil(0.6),
67  fSemigoodCorrect(0),
68  fHolePos(0),
69  fHoleWidth(0),
70  fCentSelectOn(kTRUE),
71  fCentMin(0),
72  fCentMax(10),
73  fOneConstSelectOn(kFALSE),
74  fDerivSubtrOrder(0),
75  fPhiJetCorr6(0x0),
76  fPhiJetCorr7(0x0),
77  fEtaJetCorr6(0x0),
78  fEtaJetCorr7(0x0),
79  fPtJetCorr(0x0),
80  fPtJet(0x0),
81  fhpTjetpT(0x0),
82  fhPt(0x0),
83  fhPhi(0x0),
84  fNbOfConstvspT(0x0),
85  fTreeObservableTagging(0)
86 
87 {
88  SetMakeGeneralHistograms(kTRUE);
89 }
90 
91 //________________________________________________________________________
93  AliAnalysisTaskEmcalJet(name, kTRUE),
94  fContainer(0),
95  fMinFractionShared(0),
96  fJetShapeType(kGenShapes),
97  fJetShapeSub(kNoSub),
98  fJetSelection(kInclusive),
99  fShapesVar(0),
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(0)
125 
126 {
127  // Standard constructor.
128 
130 
131  DefineOutput(1, TList::Class());
132  DefineOutput(2, TTree::Class());
133 
134 
135 }
136 
137 //________________________________________________________________________
139 {
140  // Destructor.
141 }
142 
143 //________________________________________________________________________
145 {
146  // Create user output.
147 
149 
150  Bool_t oldStatus = TH1::AddDirectoryStatus();
151  TH1::AddDirectory(kFALSE);
152 
153  //fTreeObservableTagging = new TTree("fTreeJetShape", "fTreeJetShape");
154 
155  const char* nameoutput = GetOutputSlot(2)->GetContainer()->GetName();
156  fTreeObservableTagging = new TTree(nameoutput, nameoutput);
157 
158  const Int_t nVar = 21;
159 
160  fShapesVar = new Float_t [nVar];
161  TString *fShapesVarNames = new TString [nVar];
162 
163  fShapesVarNames[0] = "partonCode";
164  fShapesVarNames[1] = "ptJet";
165  fShapesVarNames[2] = "ptDJet";
166  fShapesVarNames[3] = "mJet";
167  fShapesVarNames[4] = "nbOfConst";
168  fShapesVarNames[5] = "angularity";
169  fShapesVarNames[6] = "circularity";
170  fShapesVarNames[7] = "lesub";
171  fShapesVarNames[8] = "CoreFraction";
172  fShapesVarNames[9] = "Nsubjet1";
173  fShapesVarNames[10] = "Nsubjet2";
174  fShapesVarNames[11] = "SubjetFraction";
175  fShapesVarNames[12] = "weightPythia";
176 
177  fShapesVarNames[13] = "NT70";
178  fShapesVarNames[14] = "nConstNT70";
179  fShapesVarNames[15] = "NT80";
180  fShapesVarNames[16] = "nConstNT80";
181  fShapesVarNames[17] = "NT90";
182  fShapesVarNames[18] = "nConstNT90";
183  fShapesVarNames[19] = "NT95";
184  fShapesVarNames[20] = "nConstNT95";
185 
186 
187  for(Int_t ivar=0; ivar < nVar; ivar++){
188  cout<<"looping over variables"<<endl;
189  fTreeObservableTagging->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/F", fShapesVarNames[ivar].Data()));}
190 
191 
192 
193  fPhiJetCorr6= new TH2F("fPhiJetCorr6", "fPhiJetCorr6", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
194  fOutput->Add(fPhiJetCorr6);
195  fEtaJetCorr6= new TH2F("fEtaJetCorr6", "fEtaJetCorr6", 50, -1.5, 1.5, 50, -1.5, 1.5);
196  fOutput->Add(fEtaJetCorr6);
197 
198  fPhiJetCorr7= new TH2F("fPhiJetCorr7", "fPhiJetCorr7", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
199  fOutput->Add(fPhiJetCorr7);
200  fEtaJetCorr7= new TH2F("fEtaJetCorr7", "fEtaJetCorr7", 50, -1.5, 1.5, 50, -1.5, 1.5);
201  fOutput->Add(fEtaJetCorr7);
202 
203  fPtJetCorr= new TH2F("fPtJetCorr", "fPtJetCorr", 100, 0, 200, 100, 0, 200);
204  fOutput->Add(fPtJetCorr);
205  fPtJet= new TH1F("fPtJet", "fPtJet", 100, 0, 200);
206  fOutput->Add(fPtJet);
207 
208  fhpTjetpT= new TH2F("fhpTjetpT", "fhpTjetpT", 200, 0, 200, 200, 0, 200);
209  fOutput->Add(fhpTjetpT);
210  fhPt= new TH1F("fhPt", "fhPt", 200, 0, 200);
211  fOutput->Add(fhPt);
212  fhPhi= new TH1F("fhPhi", "fhPhi", 100, -TMath::Pi(), TMath::Pi());
213  fOutput->Add(fhPhi);
214 
215  fNbOfConstvspT=new TH2F("fNbOfConstvspT", "fNbOfConstvspT", 100, 0, 100, 200, 0, 200);
216  fOutput->Add(fNbOfConstvspT);
217 
218  //fOutput->Add(fTreeObservableTagging);
219 
220  TH1::AddDirectory(oldStatus);
221  PostData(1, fOutput); // Post data for ALL output slots > 0 here
222  PostData(2, fTreeObservableTagging);
223 
224  delete [] fShapesVarNames;
225 
226 }
227 
228 //________________________________________________________________________
230 {
231  // Run analysis code here, if needed. It will be executed before FillHistograms().
232 
233  return kTRUE;
234 }
235 
236 //________________________________________________________________________
238 {
239  // Fill histograms.
240  //cout<<"base container"<<endl;
241  AliEmcalJet* jet1 = NULL;
242  AliJetContainer *jetCont = GetJetContainer(0);
243 
244  Float_t kWeight=1;
245  if (fCentSelectOn)
246  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
247 
248  AliAODTrack *triggerHadron = 0x0;
249 
250  if (fJetSelection == kRecoil) {
251  //Printf("Recoil jets!!!, fminpTTrig = %f, fmaxpTTrig = %f", fminpTTrig, fmaxpTTrig);
252  Int_t triggerHadronLabel = SelectTrigger(fminpTTrig, fmaxpTTrig);
253 
254 
255  if (triggerHadronLabel==-99999) {
256  //Printf ("Trigger Hadron not found, return");
257  return 0;}
258 
259 
261  TClonesArray *trackArrayAn = partContAn->GetArray();
262  triggerHadron = static_cast<AliAODTrack*>(trackArrayAn->At(triggerHadronLabel));
263 
264  if (!triggerHadron) {
265  //Printf("No Trigger hadron with the found label!!");
266  return 0;
267  }
268 
269  if(fSemigoodCorrect){
270  Double_t disthole=RelativePhi(triggerHadron->Phi(),fHolePos);
271  if(TMath::Abs(disthole)+fHoleWidth>TMath::Pi()-fangWindowRecoil){
272  return 0;}
273  }
274 
275  fhPt->Fill(triggerHadron->Pt());
276 
277  }
278 
279  if(jetCont) {
280  jetCont->ResetCurrentID();
281  while((jet1 = jetCont->GetNextAcceptJet())) {
282  //Printf("jet1=%p", jet1);
283  if (!jet1) continue;
284  AliEmcalJet* jet2 = 0x0;
285  AliEmcalJet* jet3 = 0x0;
286  fPtJet->Fill(jet1->Pt());
287  AliEmcalJet *jetUS = NULL;
288  Int_t ifound=0, jfound=0;
289  Int_t ilab=-1, jlab=-1;
290 
292  Double_t disthole=RelativePhi(jet1->Phi(),fHolePos);
293  if(TMath::Abs(disthole)<fHoleWidth){
294  continue;
295  }
296  }
297 
298  Float_t dphiRecoil = 0.;
299  if (fJetSelection == kRecoil){
300  dphiRecoil = RelativePhi(triggerHadron->Phi(), jet1->Phi());
301  if (TMath::Abs(dphiRecoil) < (TMath::Pi() - fangWindowRecoil)) {
302  // Printf("Recoil jets back to back not found! continuing");
303  continue;
304  }
305 
306  fhpTjetpT->Fill(triggerHadron->Pt(), jet1->Pt());
307  //Printf(" ************ FILLING HISTOS****** shapeSub = %d, triggerHadron = %f, jet1 = %f", fJetShapeSub, triggerHadron->Pt(), jet1->Pt());
308  fhPhi->Fill(RelativePhi(triggerHadron->Phi(), jet1->Phi()));
309 
310  }
311 
312 
313  fShapesVar[0] = 0.;
314 
315  if (fJetShapeType == kGenShapes){
316  const AliEmcalPythiaInfo *partonsInfo = 0x0;
317  partonsInfo = GetPythiaInfo();
318  //Printf("partonsInfo=%p", partonsInfo);
319  Double_t jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi6());
320  Double_t detap1=(jet1->Eta())-(partonsInfo->GetPartonEta6());
321  kWeight=partonsInfo->GetPythiaEventWeight();
322  //Printf("kWeight=%f", kWeight);
323  fShapesVar[12] = kWeight;
324 
325  Float_t dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
326  fEtaJetCorr6->Fill(jet1->Eta(), partonsInfo->GetPartonEta6());
327  fPhiJetCorr6->Fill(jet1->Phi(), partonsInfo->GetPartonPhi6());
328  if(dRp1 < fRMatching) {
329  fShapesVar[0] = partonsInfo->GetPartonFlag6();
330  fPtJetCorr ->Fill(partonsInfo->GetPartonPt6(), jet1->Pt());
331  }
332  else {
333  jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi7());
334  detap1=(jet1->Eta())-(partonsInfo->GetPartonEta7());
335  dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
336  fEtaJetCorr7->Fill(jet1->Eta(), partonsInfo->GetPartonEta7());
337  fPhiJetCorr7->Fill(jet1->Phi(), partonsInfo->GetPartonPhi7());
338  if(dRp1 < fRMatching) {
339  fShapesVar[0] = partonsInfo->GetPartonFlag7();
340  fPtJetCorr->Fill(partonsInfo->GetPartonPt7(), jet1->Pt());
341  }
342  else fShapesVar[0]=0;
343  }
344  }
345 
346  Double_t ptSubtracted = 0;
347  ptSubtracted= jet1->Pt();
348  //Printf("ptSubtracted=%f", ptSubtracted);
349 
350 
351  if (ptSubtracted < fPtThreshold) continue;
352 
353  if ((fCentSelectOn == kFALSE) && (jet1->GetNumberOfTracks() <= 1)) continue;
354 
355  AliEmcalJetFinder *Reclusterer1; //Object containg Subjets from Subtracted Hybrid Jets
356  Reclusterer1 = Recluster(jet1, 0, fSubjetRadius, 0, 0, "SubJetFinder_1");
357 
358  fShapesVar[1] = ptSubtracted;
359  fShapesVar[2] = GetJetpTD(jet1,0);
360  fShapesVar[3] = GetJetMass(jet1,0);
361  fShapesVar[4] = 1.*GetJetNumberOfConstituents(jet1,0);
362  fShapesVar[5] = GetJetAngularity(jet1,0);
363  fShapesVar[6] = GetJetCircularity(jet1,0);
364  fShapesVar[7] = GetJetLeSub(jet1,0);
365  fShapesVar[8] = GetJetCoreFrac(jet1,0);
366  fShapesVar[9] = NSubJettiness(jet1, 0, fJetRadius, Reclusterer1, 1, 0, 1);
367  fShapesVar[10]= NSubJettiness(jet1, 0, fJetRadius, Reclusterer1, 2, 0, 1);
368  fShapesVar[11]= GetSubjetFraction(jet1,0,fJetRadius,Reclusterer1);
369 
370  Float_t nTFractions[8]={0.,0.,0.,0.,0.,0.,0.,0.};
371  NTValues(jet1, 0, nTFractions);
372  for (Int_t ishape=13; ishape<21; ishape++) fShapesVar[ishape] = nTFractions[ishape-13];
373 
374 
375  fTreeObservableTagging->Fill();
376 
377  }
378 
379  }
380 
381  return kTRUE;
382 }
383 
384 //________________________________________________________________________
386  //calc subtracted jet mass
387  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
388  if (fDerivSubtrOrder == 1) return jet->GetFirstOrderSubtracted();
389  else return jet->GetSecondOrderSubtracted();
390  else
391  return jet->M();
392 }
393 
394 //________________________________________________________________________
396 
397  AliJetContainer *jetCont = GetJetContainer(jetContNb);
398  if (!jet->GetNumberOfTracks())
399  return 0;
400  Double_t den=0.;
401  Double_t num = 0.;
402  AliVParticle *vp1 = 0x0;
403  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
404  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
405 
406  if (!vp1){
407  Printf("AliVParticle associated to constituent not found");
408  continue;
409  }
410 
411  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
412  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
413  Double_t dr = TMath::Sqrt(dr2);
414  num=num+vp1->Pt()*dr;
415  den=den+vp1->Pt();
416  }
417  return num/den;
418 }
419 
420 //________________________________________________________________________
422 
423  if((fJetShapeSub==kDerivSub) && (jetContNb==0))
424  if (fDerivSubtrOrder == 1) return jet->GetFirstOrderSubtractedAngularity();
425  else return jet->GetSecondOrderSubtractedAngularity();
426  else
427  return Angularity(jet, jetContNb);
428 
429 }
430 
431 
432 //________________________________________________________________________
433 Float_t AliAnalysisTaskEmcalJetShapesMC::PTD(AliEmcalJet *jet, Int_t jetContNb){
434 
435  AliJetContainer *jetCont = GetJetContainer(jetContNb);
436  if (!jet->GetNumberOfTracks())
437  return 0;
438  Double_t den=0.;
439  Double_t num = 0.;
440  AliVParticle *vp1 = 0x0;
441  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
442  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
443 
444  if (!vp1){
445  Printf("AliVParticle associated to constituent not found");
446  continue;
447  }
448 
449  num=num+vp1->Pt()*vp1->Pt();
450  den=den+vp1->Pt();
451  }
452  return TMath::Sqrt(num)/den;
453 }
454 
455 //________________________________________________________________________
457  //calc subtracted jet mass
458  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
459  if (fDerivSubtrOrder == 1) return jet->GetFirstOrderSubtractedpTD();
460  else return jet->GetSecondOrderSubtractedpTD();
461  else
462  return PTD(jet, jetContNb);
463 
464 }
465 
466 //_____________________________________________________________________________
468 
469  AliJetContainer *jetCont = GetJetContainer(jetContNb);
470  if (!jet->GetNumberOfTracks())
471  return 0;
472  Double_t mxx = 0.;
473  Double_t myy = 0.;
474  Double_t mxy = 0.;
475  int nc = 0;
476  Double_t sump2 = 0.;
477  Double_t pxjet=jet->Px();
478  Double_t pyjet=jet->Py();
479  Double_t pzjet=jet->Pz();
480 
481 
482  //2 general normalized vectors perpendicular to the jet
483  TVector3 ppJ1(pxjet, pyjet, pzjet);
484  TVector3 ppJ3(- pxjet* pzjet, - pyjet * pzjet, pxjet * pxjet + pyjet * pyjet);
485  ppJ3.SetMag(1.);
486  TVector3 ppJ2(-pyjet, pxjet, 0);
487  ppJ2.SetMag(1.);
488  AliVParticle *vp1 = 0x0;
489  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
490  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
491 
492  if (!vp1){
493  Printf("AliVParticle associated to constituent not found");
494  continue;
495  }
496 
497  TVector3 pp(vp1->Px(), vp1->Py(), vp1->Pz());
498 
499  //local frame
500  TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
501  TVector3 pPerp = pp - pLong;
502  //projection onto the two perpendicular vectors defined above
503 
504  Float_t ppjX = pPerp.Dot(ppJ2);
505  Float_t ppjY = pPerp.Dot(ppJ3);
506  Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
507  if(ppjT<=0) return 0;
508 
509  mxx += (ppjX * ppjX / ppjT);
510  myy += (ppjY * ppjY / ppjT);
511  mxy += (ppjX * ppjY / ppjT);
512  nc++;
513  sump2 += ppjT;}
514 
515  if(nc<2) return 0;
516  if(sump2==0) return 0;
517  // Sphericity Matrix
518  Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
519  TMatrixDSym m0(2,ele);
520 
521  // Find eigenvectors
522  TMatrixDSymEigen m(m0);
523  TVectorD eval(2);
524  TMatrixD evecm = m.GetEigenVectors();
525  eval = m.GetEigenValues();
526  // Largest eigenvector
527  int jev = 0;
528  // cout<<eval[0]<<" "<<eval[1]<<endl;
529  if (eval[0] < eval[1]) jev = 1;
530  TVectorD evec0(2);
531  // Principle axis
532  evec0 = TMatrixDColumn(evecm, jev);
533  Double_t compx=evec0[0];
534  Double_t compy=evec0[1];
535  TVector2 evec(compx, compy);
536  Double_t circ=0;
537  if(jev==1) circ=2*eval[0];
538  if(jev==0) circ=2*eval[1];
539 
540  return circ;
541 
542 
543 
544 }
545 
546 
547 
548 
549 //________________________________________________________________________
551  //calc subtracted jet mass
552 
553  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
555  else return jet->GetSecondOrderSubtractedCircularity();
556  else
557  return Circularity(jet, jetContNb);
558 
559 }
560 
561 //________________________________________________________________________
563 
564  AliJetContainer *jetCont = GetJetContainer(jetContNb);
565  if (!jet->GetNumberOfTracks())
566  return 0;
567  Double_t den=0.;
568  Double_t num = 0.;
569  AliVParticle *vp1 = 0x0;
570  AliVParticle *vp2 = 0x0;
571  std::vector<int> ordindex;
572  ordindex=jet->SortConstituentsPt(jetCont->GetParticleContainer()->GetArray());
573  //Printf("Nbof const = %d", jet->GetNumberOfTracks());
574  //Printf("ordindex[0] = %d, ordindex[1] = %d", ordindex[0], ordindex[1]);
575 
576  if(ordindex.size()<2) return -1;
577 
578  vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[0], jetCont->GetParticleContainer()->GetArray()));
579  if (!vp1){
580  Printf("AliVParticle associated to Leading constituent not found");
581  return -1;
582  }
583 
584  vp2 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[1], jetCont->GetParticleContainer()->GetArray()));
585  if (!vp2){
586  Printf("AliVParticle associated to Subleading constituent not found");
587  return -1;
588  }
589 
590 
591  num=vp1->Pt();
592  den=vp2->Pt();
593  //Printf("vp1->Pt() =%f, vp2->Pt() =%f", vp1->Pt(), vp2->Pt());
594 
595 return num-den;
596 }
597 
598 //________________________________________________________________________
600  //calc subtracted jet mass
601 
602  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
603  if (fDerivSubtrOrder == 1) return jet->GetFirstOrderSubtractedLeSub();
604  else return jet->GetSecondOrderSubtractedLeSub();
605  else
606  return LeSub(jet, jetContNb);
607 
608 }
609 
610 //________________________________________________________________________
612  //calc subtracted jet mass
613 
614  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
616  else return jet->GetSecondOrderSubtractedConstituent();
617  else
618  return jet->GetNumberOfTracks();
619 
620  }
621 
622 
623 //________________________________________________________________________
624 AliEmcalJetFinder *AliAnalysisTaskEmcalJetShapesMC::Recluster(AliEmcalJet *Jet, Int_t JetContNb, Double_t SubJetRadius, Double_t SubJetMinPt, Int_t Algorithm, const char* Name){
625 
626  AliJetContainer *JetCont = GetJetContainer(JetContNb);
627  AliEmcalJetFinder *Reclusterer = new AliEmcalJetFinder(Name); //JetFinder Object for reclustered jets
628  Reclusterer->SetRadius(SubJetRadius);
629  Reclusterer->SetJetMinPt(SubJetMinPt);
630  Reclusterer->SetJetAlgorithm(Algorithm); //0 for anti-kt 1 for kt
631  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
632  Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
633  if(Reclusterer->AliEmcalJetFinder::Filter(Jet, JetCont, dVtx)){;} //reclustering jet1 using the jetfinderobject Reclusterer
634  return Reclusterer;
635 }
636 
637 
638 
639 //________________________________________________________________________
640 Double_t AliAnalysisTaskEmcalJetShapesMC::NSubJettiness(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t A, Int_t B){
641  AliJetContainer *JetCont = GetJetContainer(JetContNb);
642  AliEmcalJet *SubJet=NULL;
643  Double_t DeltaR1=0;
644  Double_t DeltaR2=0;
645  AliVParticle *JetParticle=0x0;
646  Double_t SubJetiness_Numerator = 0;
647  Double_t SubJetiness_Denominator = 0;
648  Double_t Index=-2;
649  Bool_t Error=kFALSE;
650  if (!Jet->GetNumberOfTracks()) return -2;
651  if (Reclusterer->GetNumberOfJets() < N) return -2;
652  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
653  JetParticle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
654  for (Int_t j=1; j<=N; j++){
655  Index=SubJetOrdering(Jet,Reclusterer,j,0,kTRUE);
656  if(Index==-999){
657  Error = kTRUE;
658  i=Jet->GetNumberOfTracks();
659  break;
660  }
661  if(j==1){
662  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);
663  }
664  else{
665  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);
666  if (DeltaR2<DeltaR1) DeltaR1=DeltaR2;
667  }
668  }
669  SubJetiness_Numerator=SubJetiness_Numerator+(JetParticle->Pt()*DeltaR1);
670  if (A>=0) SubJetiness_Denominator=SubJetiness_Denominator+(TMath::Power((Reclusterer->GetJet(SubJetOrdering(Jet,Reclusterer,1,0,kTRUE))->Pt()),A)*JetParticle->Pt()*TMath::Power(JetRadius,B));
671  else SubJetiness_Denominator=SubJetiness_Denominator+(TMath::Power((Reclusterer->GetJet(SubJetOrdering(Jet,Reclusterer,N,0,kTRUE))->Pt()),A)*JetParticle->Pt()*TMath::Power(JetRadius,B));
672  }
673  if (SubJetiness_Denominator!=0 && !Error) return SubJetiness_Numerator/SubJetiness_Denominator;
674  else return -2;
675 
676 }
677 
678 //----------------------------------------------------------------------
679 Double_t AliAnalysisTaskEmcalJetShapesMC::GetSubjetFraction(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, AliEmcalJetFinder *Reclusterer){
680  AliJetContainer *JetCont = GetJetContainer(JetContNb);
681  AliEmcalJet *SubJet=NULL;
682  Double_t DeltaR1=0;
683  Double_t DeltaR2=0;
684  AliVParticle *JetParticle=0x0;
685  Double_t SubJetiness_Numerator = 0;
686  Double_t SubJetiness_Denominator = 0;
687  Double_t Index=-2;
688  if (Reclusterer->GetNumberOfJets() < 1) return -2;
689  Index=SubJetOrdering(Jet,Reclusterer,1,0,kTRUE);
690  if(Index==-999) return -2;
691  SubJetiness_Numerator=(Reclusterer->GetJet(Index)->Pt());
692  SubJetiness_Denominator=Jet->Pt();
693  return SubJetiness_Numerator/SubJetiness_Denominator;
694 
695 
696 }
697 //__________________________________________________________________________________
699 
700  AliJetContainer *jetCont = GetJetContainer(jetContNb);
701  if (!jet->GetNumberOfTracks())
702  return 0;
703  Double_t den=0.;
704  Double_t num = 0.;
705  AliVParticle *vp1 = 0x0;
706  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
707  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
708 
709  if (!vp1){
710  Printf("AliVParticle associated to constituent not found");
711  continue;
712  }
713 
714  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
715  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
716  Double_t dr = TMath::Sqrt(dr2);
717  if(dr<=fSubjetRadius) num=num+vp1->Pt();
718 
719  }
720  return num/jet->Pt();
721 }
722 
723 
724 
725 
726 //________________________________________________________________________
728  //calc subtracted jet mass
729 
730  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
731  if (fDerivSubtrOrder == 1) return jet->GetFirstOrderSubtractedLeSub();
732  else return jet->GetSecondOrderSubtractedLeSub();
733  else
734  return CoreFrac(jet, jetContNb);
735 
736 }
737 
738 
739 
740 
741 //----------------------------------------------------------------------
742 Double_t AliAnalysisTaskEmcalJetShapesMC::SubJetOrdering(AliEmcalJet *Jet, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t Type, Bool_t Index){
743  AliEmcalJet *SubJet=NULL;
744  Double_t SortingVariable;
745  Int_t ArraySize =N+1;
746  TArrayD *JetSorter = new TArrayD(ArraySize);
747  TArrayD *JetIndexSorter = new TArrayD(ArraySize);
748  for (Int_t i=0; i<ArraySize; i++){
749  JetSorter->SetAt(0,i);
750  }
751  for (Int_t i=0; i<ArraySize; i++){
752  JetIndexSorter->SetAt(0,i);
753  }
754  if(Reclusterer->GetNumberOfJets()<N) return -999;
755  for (Int_t i=0; i<Reclusterer->GetNumberOfJets(); i++){
756  SubJet=Reclusterer->GetJet(i);
757  if (Type==0) SortingVariable=SubJet->Pt();
758  else if (Type==1) SortingVariable=SubJet->E();
759  else if (Type==2) SortingVariable=SubJet->M();
760  for (Int_t j=0; j<N; j++){
761  if (SortingVariable>JetSorter->GetAt(j)){
762  for (Int_t k=N-1; k>=j; k--){
763  JetSorter->SetAt(JetSorter->GetAt(k),k+1);
764  JetIndexSorter->SetAt(JetIndexSorter->GetAt(k),k+1);
765  }
766  JetSorter->SetAt(SortingVariable,j);
767  JetIndexSorter->SetAt(i,j);
768  break;
769  }
770  }
771  }
772  if (!Index) return JetSorter->GetAt(N-1);
773  else return JetIndexSorter->GetAt(N-1);
774 }
775 
776 
777 
778 //returns -1 if the Nth hardest jet is requested where N>number of available jets
779 //type: 0=Pt 1=E 2=M
780 //Index TRUE=returns index FALSE=returns value of quantatiy in question
781 
782 
783 
784 
785 
786 //______________________________________________________________________________
788 
789  AliJetContainer *jetCont = GetJetContainer(jetContNb);
790  if (!jet->GetNumberOfTracks())
791  return 0;
792  Double_t mxx = 0.;
793  Double_t myy = 0.;
794  Double_t mxy = 0.;
795  int nc = 0;
796  Double_t sump2 = 0.;
797 
798  AliVParticle *vp1 = 0x0;
799  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
800  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
801 
802  if (!vp1){
803  Printf("AliVParticle associated to constituent not found");
804  continue;
805  }
806 
807  Double_t ppt=vp1->Pt();
808  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
809 
810  Double_t deta = vp1->Eta()-jet->Eta();
811  mxx += ppt*ppt*deta*deta;
812  myy += ppt*ppt*dphi*dphi;
813  mxy -= ppt*ppt*deta*TMath::Abs(dphi);
814  nc++;
815  sump2 += ppt*ppt;
816 
817  }
818  if(nc<2) return 0;
819  if(sump2==0) return 0;
820  // Sphericity Matrix
821  Double_t ele[4] = {mxx , mxy , mxy , myy };
822  TMatrixDSym m0(2,ele);
823 
824  // Find eigenvectors
825  TMatrixDSymEigen m(m0);
826  TVectorD eval(2);
827  TMatrixD evecm = m.GetEigenVectors();
828  eval = m.GetEigenValues();
829  // Largest eigenvector
830  int jev = 0;
831  // cout<<eval[0]<<" "<<eval[1]<<endl;
832  if (eval[0] < eval[1]) jev = 1;
833  TVectorD evec0(2);
834  // Principle axis
835  evec0 = TMatrixDColumn(evecm, jev);
836  Double_t compx=evec0[0];
837  Double_t compy=evec0[1];
838  TVector2 evec(compx, compy);
839  Double_t sig=0;
840  if(jev==1) sig=TMath::Sqrt(TMath::Abs(eval[0])/sump2);
841  if(jev==0) sig=TMath::Sqrt(TMath::Abs(eval[1])/sump2);
842 
843  return sig;
844 
845 }
846 
847 //________________________________________________________________________
849  //calc subtracted jet mass
850 
851  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
852  if (fDerivSubtrOrder == 1) return jet->GetFirstOrderSubtractedSigma2();
853  else return jet->GetSecondOrderSubtractedSigma2();
854  else
855  return Sigma2(jet, jetContNb);
856 
857 }
858 
859 
860 //_________________________________________________________________________
861 void AliAnalysisTaskEmcalJetShapesMC::NTValues(AliEmcalJet *jet, Int_t jetContNb, Float_t* nTFractions){
862 
863  AliJetContainer *jetCont = GetJetContainer(jetContNb);
864  if (!jet->GetNumberOfTracks())
865  return;
866 
867  Double_t ptJet = jet->Pt();
868 
869  AliVParticle *vp1 = 0x0;
870  std::vector<int> ordindex;
871  ordindex=jet->SortConstituentsPt(jetCont->GetParticleContainer()->GetArray());
872  //Printf("Nbof const = %d", jet->GetNumberOfTracks());
873  //Printf("ordindex[0] = %d, ordindex[1] = %d", ordindex[0], ordindex[1]);
874  //if(ordindex.size()<2) return -1;
875 
876  for(Int_t iconst =0; iconst<jet->GetNumberOfTracks(); iconst++){
877 
878  vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[iconst], jetCont->GetParticleContainer()->GetArray()));
879  if (!vp1){
880  Printf("AliVParticle associated to Leading constituent not found");
881  return;
882  }
883 
884  if (nTFractions[0] <= 0.7*ptJet){
885  nTFractions[0] += vp1->Pt();
886  nTFractions[1] +=1;
887  }
888 
889  if (nTFractions[2] <= 0.8*ptJet){
890  nTFractions[2] += vp1->Pt();
891  nTFractions[3] +=1;
892  }
893 
894  if (nTFractions[4] <= 0.9*ptJet){
895  nTFractions[4] += vp1->Pt();
896  nTFractions[5] +=1;
897  }
898 
899  if (nTFractions[6] <= 0.95*ptJet){
900  nTFractions[6] += vp1->Pt();
901  nTFractions[7] +=1;
902  }
903  }
904 }
905 
906 
907 //________________________________________________________________________
908 Int_t AliAnalysisTaskEmcalJetShapesMC::SelectTrigger(Float_t minpT, Float_t maxpT){
909 
911  TClonesArray *tracksArray = partCont->GetArray();
912 
913  if(!partCont || !tracksArray) return -99999;
914  AliAODTrack *track = 0x0;
915  AliEmcalParticle *emcPart = 0x0;
916 
917 
918  TList *trackList = new TList();
919  Int_t triggers[100];
920  for (Int_t iTrigger=0; iTrigger<100; iTrigger++) triggers[iTrigger] = 0;
921  Int_t iTT = 0;
922 
923  for(Int_t iTrack=0; iTrack <= tracksArray->GetEntriesFast(); iTrack++){
924 
925 
926  if (fJetShapeSub == kConstSub){
927  emcPart = static_cast<AliEmcalParticle*>(tracksArray->At(iTrack));
928  if (!emcPart) continue;
929  if(TMath::Abs(emcPart->Eta())>0.9) continue;
930  if (emcPart->Pt()<0.15) continue;
931 
932  if ((emcPart->Pt() >= minpT) && (emcPart->Pt()< maxpT)) {
933  trackList->Add(emcPart);
934  triggers[iTT] = iTrack;
935  iTT++;
936  }
937  }
938  else{
939  track = static_cast<AliAODTrack*>(tracksArray->At(iTrack));
940  if (!track) continue;
941  if(TMath::Abs(track->Eta())>0.9) continue;
942  if (track->Pt()<0.15) continue;
943  if (!(track->TestFilterBit(768))) continue;
944 
945  if ((track->Pt() >= minpT) && (track->Pt()< maxpT)) {
946  trackList->Add(track);
947  triggers[iTT] = iTrack;
948  iTT++;
949 
950  }
951  }
952  }
953 
954  if (iTT == 0) return -99999;
955  Int_t nbRn = 0, index = 0 ;
956  TRandom3* random = new TRandom3(0);
957  nbRn = random->Integer(iTT);
958 
959  index = triggers[nbRn];
960  //Printf("iTT Total= %d, nbRn = %d, Index = %d",iTT, nbRn, index );
961  return index;
962 
963 }
964 
965 //__________________________________________________________________________________
966 Double_t AliAnalysisTaskEmcalJetShapesMC::RelativePhi(Double_t mphi,Double_t vphi){
967 
968  if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
969  else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
970  if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
971  else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
972  double dphi = mphi-vphi;
973  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
974  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
975  return dphi;//dphi in [-Pi, Pi]
976 }
977 
978 
979 //________________________________________________________________________
981  //
982  // retrieve event objects
983  //
985  return kFALSE;
986 
987  return kTRUE;
988 }
989 
990 //_______________________________________________________________________
992 {
993  // Called once at the end of the analysis.
994 
995  // fTreeObservableTagging = dynamic_cast<TTree*>(GetOutputData(1));
996  // if (!fTreeObservableTagging){
997  // Printf("ERROR: fTreeObservableTagging not available");
998  // return;
999  // }
1000 
1001 }
1002 
void SetRadius(Double_t val)
Float_t GetSigma2(AliEmcalJet *jet, Int_t jetContNb=0)
Float_t GetJetMass(AliEmcalJet *jet, Int_t jetContNb=0)
void SetJetMinPt(Double_t val)
Float_t GetPartonEta6() const
AliJetContainer * GetJetContainer(Int_t i=0) const
Float_t GetPartonEta7() const
Double_t Eta() const
Definition: AliEmcalJet.h:60
ClassImp(AliAnalysisTaskEmcalJetShapesMC) AliAnalysisTaskEmcalJetShapesMC
Double_t Py() const
Definition: AliEmcalJet.h:45
Double_t Phi() const
Definition: AliEmcalJet.h:55
Float_t GetPythiaEventWeight() const
Double_t GetSecondOrderSubtractedLeSub() const
Definition: AliEmcalJet.h:256
Double_t GetFirstOrderSubtractedConstituent() const
Definition: AliEmcalJet.h:245
Float_t GetPartonPhi7() const
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:58
TList * fOutput
!output list
Double_t GetFirstOrderSubtracted() const
Definition: AliEmcalJet.h:178
Float_t GetJetpTD(AliEmcalJet *jet, Int_t jetContNb=0)
Int_t GetPartonFlag7() const
Container for particles within the EMCAL framework.
Float_t GetPartonPhi6() const
Double_t GetFirstOrderSubtractedLeSub() const
Definition: AliEmcalJet.h:255
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:83
Double_t Px() const
Definition: AliEmcalJet.h:44
Float_t GetJetCoreFrac(AliEmcalJet *jet, Int_t jetContNb=0)
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Int_t GetPartonFlag6() const
AliParticleContainer * GetParticleContainer() const
void SetJetAlgorithm(Int_t val)
Double_t Eta() const
Double_t Pt() const
Float_t GetJetNumberOfConstituents(AliEmcalJet *jet, Int_t jetContNb=0)
Double_t GetSecondOrderSubtractedSigma2() const
Definition: AliEmcalJet.h:235
std::vector< int > SortConstituentsPt(TClonesArray *tracks) const
Double_t RelativePhi(Double_t mphi, Double_t vphi)
Double_t GetSubjetFraction(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, AliEmcalJetFinder *Reclusterer)
Double_t fCent
!event centrality
TClonesArray * GetArray() const
Float_t PTD(AliEmcalJet *jet, Int_t jetContNb=0)
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 GetSecondOrderSubtracted() const
Definition: AliEmcalJet.h:179
Double_t Pt() const
Definition: AliEmcalJet.h:47
Double_t GetSecondOrderSubtractedpTD() const
Definition: AliEmcalJet.h:215
Float_t Angularity(AliEmcalJet *jet, Int_t jetContNb=0)
Double_t GetFirstOrderSubtractedpTD() const
Definition: AliEmcalJet.h:214
Double_t SubJetOrdering(AliEmcalJet *Jet, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t Type, Bool_t Index)
Double_t GetFirstOrderSubtractedAngularity() const
Definition: AliEmcalJet.h:204
Short_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:107
void NTValues(AliEmcalJet *jet, Int_t jetContNb, Float_t *nTFractions)
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)
Double_t GetFirstOrderSubtractedCircularity() const
Definition: AliEmcalJet.h:224
const AliEmcalPythiaInfo * GetPythiaInfo() const
Double_t Pz() const
Definition: AliEmcalJet.h:46
Float_t Circularity(AliEmcalJet *jet, Int_t jetContNb=0)
AliEmcalJetFinder * Recluster(AliEmcalJet *Jet, Int_t JetContNb, Double_t SubJetRadius, Double_t SubJetMinPt, Int_t Algorithm, const char *Name)
Declaration of class AliEmcalPythiaInfo.
Double_t GetSecondOrderSubtractedAngularity() const
Definition: AliEmcalJet.h:205
Float_t Sigma2(AliEmcalJet *jet, Int_t jetContNb=0)
Double_t GetSecondOrderSubtractedCircularity() const
Definition: AliEmcalJet.h:225
Double_t NSubJettiness(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t A, Int_t B)
Float_t GetPartonPt7() const
Float_t LeSub(AliEmcalJet *jet, Int_t jetContNb=0)
Double_t M() const
Definition: AliEmcalJet.h:59
void ResetCurrentID(Int_t i=-1)
Double_t GetSecondOrderSubtractedConstituent() const
Definition: AliEmcalJet.h:246
Double_t GetFirstOrderSubtractedSigma2() const
Definition: AliEmcalJet.h:234
Float_t GetPartonPt6() const