AliPhysics  c923f52 (c923f52)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskFakeJets.cxx
Go to the documentation of this file.
1 //
2 // Trying to see if fake jets can be identified by their shape, collection of selected shapes
3 //
4 // Author: L. Cunqueiro
5 // n-subjetiness code copied from Nima Zhardosti's task AliAnalysisTaskSubetFraction
6 
7 #include <TClonesArray.h>
8 #include <TH1F.h>
9 #include <TH2F.h>
10 #include <TH3F.h>
11 #include <THnSparse.h>
12 #include <TTree.h>
13 #include <TList.h>
14 #include <TLorentzVector.h>
15 #include <TProfile.h>
16 #include <TChain.h>
17 #include <TSystem.h>
18 #include <TFile.h>
19 #include <TKey.h>
20 #include "TMatrixD.h"
21 #include "TMatrixDSym.h"
22 #include "TMatrixDSymEigen.h"
23 #include "TVector3.h"
24 #include "TVector2.h"
25 #include "AliVCluster.h"
26 #include "AliVTrack.h"
27 #include "AliEmcalJet.h"
28 #include "AliRhoParameter.h"
29 #include "AliLog.h"
30 #include "AliEmcalParticle.h"
31 #include "AliMCEvent.h"
32 #include "AliGenPythiaEventHeader.h"
33 #include "AliAODMCHeader.h"
34 #include "AliMCEvent.h"
35 #include "AliAnalysisManager.h"
36 #include "AliJetContainer.h"
37 #include "AliParticleContainer.h"
38 #include "TRandom3.h"
39 #include "AliEmcalJetFinder.h"
40 #include "AliAODEvent.h"
42 
43 using std::cout;
44 using std::endl;
45 
47 
48 //________________________________________________________________________
51  fContainer(0),
52  fMinFractionShared(0),
53  fJetShapeType(kData),
54  fJetShapeSub(kNoSub),
55  fJetSelection(kInclusive),
56  fShapesVar(0),
57  fPtThreshold(-9999.),
58  fRMatching(0.2),
59  fminpTTrig(20.),
60  fmaxpTTrig(50.),
61  fangWindowRecoil(0.6),
62  fSemigoodCorrect(0),
63  fHolePos(0),
64  fHoleWidth(0),
65  fCentSelectOn(kTRUE),
66  fCentMin(0),
67  fCentMax(10),
68  fOneConstSelectOn(kFALSE),
69  fDerivSubtrOrder(0),
70  fSubjetRadius(0.2),
71  fJetRadius(0.4),
72  fh2ResponseUW(0x0),
73  fh2ResponseW(0x0),
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  fTreeFakeJets(0)
85 
86 {
87  SetMakeGeneralHistograms(kTRUE);
88 }
89 
90 //________________________________________________________________________
92  AliAnalysisTaskEmcalJet(name, kTRUE),
93  fContainer(0),
94  fMinFractionShared(0),
95  fJetShapeType(kData),
96  fJetShapeSub(kNoSub),
97  fJetSelection(kInclusive),
98  fShapesVar(0),
99  fPtThreshold(-9999.),
100  fRMatching(0.2),
101  fminpTTrig(20.),
102  fmaxpTTrig(50.),
103  fangWindowRecoil(0.6),
104  fSemigoodCorrect(0),
105  fHolePos(0),
106  fHoleWidth(0),
107  fCentSelectOn(kTRUE),
108  fCentMin(0),
109  fCentMax(10),
110  fOneConstSelectOn(kFALSE),
111  fDerivSubtrOrder(0),
112  fSubjetRadius(0.2),
113  fJetRadius(0.4),
114  fh2ResponseUW(0x0),
115  fh2ResponseW(0x0),
116  fPhiJetCorr6(0x0),
117  fPhiJetCorr7(0x0),
118  fEtaJetCorr6(0x0),
119  fEtaJetCorr7(0x0),
120  fPtJetCorr(0x0),
121  fPtJet(0x0),
122  fhpTjetpT(0x0),
123  fhPt(0x0),
124  fhPhi(0x0),
125  fNbOfConstvspT(0x0),
126  fTreeFakeJets(0)
127 
128 {
129  // Standard constructor.
130 
132 
133  DefineOutput(1, TTree::Class());
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  fTreeFakeJets = new TTree("fTreeFakeJets", "fTreeFakeJets");
154  const Int_t nVar = 11;
155  fShapesVar = new Float_t [nVar];
156  TString *fShapesVarNames = new TString [nVar];
157 
158 
159  fShapesVarNames[0] = "ptJet";
160  fShapesVarNames[1] = "ptDJet";
161  fShapesVarNames[2] = "angularity";
162  fShapesVarNames[3] = "angularitysquared";
163  fShapesVarNames[4]= "hardtrack";
164  fShapesVarNames[5]="hard2track";
165  fShapesVarNames[6]="corefrac";
166  fShapesVarNames[7]="nsubjet1";
167  fShapesVarNames[8]="nsubjet2";
168  fShapesVarNames[9]="subjetfrac";
169  fShapesVarNames[10]="mass";
170  for(Int_t ivar=0; ivar < nVar; ivar++){
171  cout<<"looping over variables"<<endl;
172  fTreeFakeJets->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/F", fShapesVarNames[ivar].Data()));
173 
174  //if( ivar == 4 ) fTreeObservableTagging->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/I", fShapesVarNames[ivar].Data()));
175 
176  }
177 
178  fh2ResponseUW= new TH2F("fh2ResponseUW", "fh2ResponseUW", 100, 0, 200, 100, 0, 200);
179  fOutput->Add(fh2ResponseUW);
180  fh2ResponseW= new TH2F("fh2ResponseW", "fh2ResponseW", 100, 0, 200, 100, 0, 200);
181  fOutput->Add(fh2ResponseW);
182  fPhiJetCorr6= new TH2F("fPhiJetCorr6", "fPhiJetCorr6", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
183  fOutput->Add(fPhiJetCorr6);
184  fEtaJetCorr6= new TH2F("fEtaJetCorr6", "fEtaJetCorr6", 50, -1.5, 1.5, 50, -1.5, 1.5);
185  fOutput->Add(fEtaJetCorr6);
186 
187  fPhiJetCorr7= new TH2F("fPhiJetCorr7", "fPhiJetCorr7", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
188  fOutput->Add(fPhiJetCorr7);
189  fEtaJetCorr7= new TH2F("fEtaJetCorr7", "fEtaJetCorr7", 50, -1.5, 1.5, 50, -1.5, 1.5);
190  fOutput->Add(fEtaJetCorr7);
191 
192  fPtJetCorr= new TH2F("fPtJetCorr", "fPtJetCorr", 100, 0, 200, 100, 0, 200);
193  fOutput->Add(fPtJetCorr);
194  fPtJet= new TH1F("fPtJet", "fPtJet", 100, 0, 200);
195  fOutput->Add(fPtJet);
196 
197  fhpTjetpT= new TH2F("fhpTjetpT", "fhpTjetpT", 200, 0, 200, 200, 0, 200);
198  fOutput->Add(fhpTjetpT);
199  fhPt= new TH1F("fhPt", "fhPt", 200, 0, 200);
200  fOutput->Add(fhPt);
201  fhPhi= new TH1F("fhPhi", "fhPhi", 100, -TMath::Pi(), TMath::Pi());
202  fOutput->Add(fhPhi);
203 
204  fNbOfConstvspT=new TH2F("fNbOfConstvspT", "fNbOfConstvspT", 100, 0, 100, 200, 0, 200);
205  fOutput->Add(fNbOfConstvspT);
206 
207  fOutput->Add(fTreeFakeJets);
208  TH1::AddDirectory(oldStatus);
209  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
210 
211 }
212 
213 //________________________________________________________________________
215 {
216  // Run analysis code here, if needed. It will be executed before FillHistograms().
217 
218  return kTRUE;
219 }
220 
221 //________________________________________________________________________
223 {
224  // Fill histograms.
225  //cout<<"base container"<<endl;
226  AliEmcalJet* jet1 = NULL;
227  AliJetContainer *jetCont = GetJetContainer(0);
228  Float_t kWeight=1;
229  if (fCentSelectOn)
230  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
231 
232  AliAODTrack *triggerHadron = 0x0;
233 
234  if (fJetSelection == kRecoil) {
235  //Printf("Recoil jets!!!, fminpTTrig = %f, fmaxpTTrig = %f", fminpTTrig, fmaxpTTrig);
236  Int_t triggerHadronLabel = SelectTrigger(fminpTTrig, fmaxpTTrig);
237 
238 
239  if (triggerHadronLabel==-99999) {
240  //Printf ("Trigger Hadron not found, return");
241  return 0;}
242 
243 
245  TClonesArray *trackArrayAn = partContAn->GetArray();
246  triggerHadron = static_cast<AliAODTrack*>(trackArrayAn->At(triggerHadronLabel));
247 
248  if (!triggerHadron) {
249  //Printf("No Trigger hadron with the found label!!");
250  return 0;
251  }
252 
253  if(fSemigoodCorrect){
254  Double_t disthole=RelativePhi(triggerHadron->Phi(),fHolePos);
255  if(TMath::Abs(disthole)+fHoleWidth>TMath::Pi()-fangWindowRecoil){
256  return 0;}
257  }
258 
259  fhPt->Fill(triggerHadron->Pt());
260 
261  }
262 
263  if(jetCont) {
264  jetCont->ResetCurrentID();
265  while((jet1 = jetCont->GetNextAcceptJet())) {
266  if (!jet1) continue;
267 
268  fPtJet->Fill(jet1->Pt());
269  AliEmcalJet *jetUS = NULL;
270  Int_t ifound=0;
271  Int_t ilab=-1;
272 
274  Double_t disthole=RelativePhi(jet1->Phi(),fHolePos);
275  if(TMath::Abs(disthole)<fHoleWidth){
276  continue;}
277  }
278 
279 
280 
281 
282  Double_t ptSubtracted = jet1->Pt();
283  if (ptSubtracted < fPtThreshold) continue;
284 
285  if (fOneConstSelectOn == kTRUE) fNbOfConstvspT->Fill(GetJetNumberOfConstituents(jet1,0), ptSubtracted);
286 
287  if ((fCentSelectOn == kFALSE) && (jet1->GetNumberOfTracks() <= 1)) continue;
288 
289  AliEmcalJetFinder *Reclusterer1; //Object containg Subjets from Subtracted Hybrid Jets
290  Reclusterer1 = Recluster(jet1, 0, fSubjetRadius, 0, 0, "SubJetFinder_1");
291 
292 
293  fShapesVar[0] = ptSubtracted;
294  fShapesVar[1] = GetJetpTD(jet1,0);
295  fShapesVar[2] = GetJetAngularity(jet1,0);
296  fShapesVar[3] = GetJetAngularitySquared(jet1,0);
297  fShapesVar[4]= GetJetHardTrack(jet1,0);
298  fShapesVar[5]=GetJetSecHardTrack(jet1,0);
299  fShapesVar[6]=GetJetCoreFrac(jet1,0);
300  fShapesVar[7]=NSubJettiness(jet1, 0, fJetRadius, Reclusterer1, 1, 0, 1);
301  fShapesVar[8]=NSubJettiness(jet1, 0, fJetRadius, Reclusterer1, 2, 0, 1);
302  fShapesVar[9]=GetSubjetFraction(jet1,0,fJetRadius,Reclusterer1);
303 
304  fShapesVar[10]=GetJetMass(jet1,0);
305  fTreeFakeJets->Fill();
306 
307  }}
308 
309 
310 
311  return kTRUE;
312  }
313 
314 //________________________________________________________________________
316  //calc subtracted jet mass
317  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
319  else return jet->GetShapeProperties()->GetSecondOrderSubtracted();
320  else
321  return jet->M();
322 }
323 
324 //________________________________________________________________________
326 
327  AliJetContainer *jetCont = GetJetContainer(jetContNb);
328  if (!jet->GetNumberOfTracks())
329  return 0;
330  Double_t den=0.;
331  Double_t num = 0.;
332  AliVParticle *vp1 = 0x0;
333  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
334  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
335 
336  if (!vp1){
337  Printf("AliVParticle associated to constituent not found");
338  continue;
339  }
340 
341  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
342  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
343  Double_t dr = TMath::Sqrt(dr2);
344  num=num+vp1->Pt()*dr;
345  den=den+vp1->Pt();
346  }
347  return num/den;
348 }
349 
350 //________________________________________________________________________
352 
353  if((fJetShapeSub==kDerivSub) && (jetContNb==0))
356  else
357  return Angularity(jet, jetContNb);
358 
359 }
360 //________________________________________________________________________
362 
363  AliJetContainer *jetCont = GetJetContainer(jetContNb);
364  if (!jet->GetNumberOfTracks())
365  return 0;
366  Double_t den=0.;
367  Double_t num = 0.;
368  AliVParticle *vp1 = 0x0;
369  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
370  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
371 
372  if (!vp1){
373  Printf("AliVParticle associated to constituent not found");
374  continue;
375  }
376 
377  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
378  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
379  Double_t dr = TMath::Sqrt(dr2);
380  num=num+vp1->Pt()*TMath::Sqrt(dr);
381  den=den+vp1->Pt();
382  }
383  return num/den;
384 }
385 
386 //________________________________________________________________________
388 
389  if((fJetShapeSub==kDerivSub) && (jetContNb==0))
392  else
393  return AngularitySquared(jet, jetContNb);
394 
395 }
396 
397 
398 //________________________________________________________________________
400 
401  AliJetContainer *jetCont = GetJetContainer(jetContNb);
402  if (!jet->GetNumberOfTracks())
403  return 0;
404  Double_t den=0.;
405  Double_t num = 0.;
406  AliVParticle *vp1 = 0x0;
407  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
408  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
409 
410  if (!vp1){
411  Printf("AliVParticle associated to constituent not found");
412  continue;
413  }
414 
415  num=num+vp1->Pt()*vp1->Pt();
416  den=den+vp1->Pt();
417  }
418  return TMath::Sqrt(num)/den;
419 }
420 
421 //________________________________________________________________________
423  //calc subtracted jet mass
424  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
426  else return jet->GetShapeProperties()->GetSecondOrderSubtractedpTD();
427  else
428  return PTD(jet, jetContNb);
429 
430 }
431 
432 
433 
434 
436 
437  AliJetContainer *jetCont = GetJetContainer(jetContNb);
438  if (!jet->GetNumberOfTracks())
439  return 0;
440  Double_t den=0.;
441  Double_t num = 0.;
442  AliVParticle *vp1 = 0x0;
443 
444  std::vector<int> ordindex;
445  ordindex=jet->GetPtSortedTrackConstituentIndexes(jetCont->GetParticleContainer()->GetArray());
446  if(ordindex.size()<1) return -1;
447 
448  vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[0], jetCont->GetParticleContainer()->GetArray()));
449  if (!vp1){
450  Printf("AliVParticle associated to Leading constituent not found");
451  return -1;
452  }
453 
454 
455  num=vp1->Pt();
456 
457 
458 
459 return num;
460 }
461 
462 //________________________________________________________________________
464 
465  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
468  else
469  return HardTrack(jet, jetContNb);
470 
471 }
472 
473 
474 //________________________________________________________________________
476 
477  AliJetContainer *jetCont = GetJetContainer(jetContNb);
478  if (!jet->GetNumberOfTracks())
479  return 0;
480  Double_t den=0.;
481  Double_t num = 0.;
482  AliVParticle *vp1 = 0x0;
483  AliVParticle *vp2 = 0x0;
484  std::vector<int> ordindex;
485  ordindex=jet->GetPtSortedTrackConstituentIndexes(jetCont->GetParticleContainer()->GetArray());
486 
487 
488  if(ordindex.size()<2) return -1;
489 
490 
491  vp2 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[1], jetCont->GetParticleContainer()->GetArray()));
492  if (!vp2){
493  Printf("AliVParticle associated to Subleading constituent not found");
494  return -1;
495  }
496 
497 
498  den=vp2->Pt();
499  //Printf("vp1->Pt() =%f, vp2->Pt() =%f", vp1->Pt(), vp2->Pt());
500 
501 return den;
502 }
503 
504 //________________________________________________________________________
506  //calc subtracted jet mass
507 
508  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
511  else
512  return SecHardTrack(jet, jetContNb);
513 
514 }
515 
516 //----------------------------------------------------------------------
518  AliJetContainer *JetCont = GetJetContainer(JetContNb);
519  AliEmcalJet *SubJet=NULL;
520  Double_t DeltaR1=0;
521  Double_t DeltaR2=0;
522  AliVParticle *JetParticle=0x0;
523  Double_t SubJetiness_Numerator = 0;
524  Double_t SubJetiness_Denominator = 0;
525  Double_t Index=-2;
526  Bool_t Error=kFALSE;
527  if (!Jet->GetNumberOfTracks()) return -2;
528  if (Reclusterer->GetNumberOfJets() < N) return -2;
529  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
530  JetParticle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
531  for (Int_t j=1; j<=N; j++){
532  Index=SubJetOrdering(Jet,Reclusterer,j,0,kTRUE);
533  if(Index==-999){
534  Error = kTRUE;
535  i=Jet->GetNumberOfTracks();
536  break;
537  }
538  if(j==1){
539  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);
540  }
541  else{
542  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);
543  if (DeltaR2<DeltaR1) DeltaR1=DeltaR2;
544  }
545  }
546  SubJetiness_Numerator=SubJetiness_Numerator+(JetParticle->Pt()*DeltaR1);
547  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));
548  else SubJetiness_Denominator=SubJetiness_Denominator+(TMath::Power((Reclusterer->GetJet(SubJetOrdering(Jet,Reclusterer,N,0,kTRUE))->Pt()),A)*JetParticle->Pt()*TMath::Power(JetRadius,B));
549  }
550  if (SubJetiness_Denominator!=0 && !Error) return SubJetiness_Numerator/SubJetiness_Denominator;
551  else return -2;
552 }
553 //____________________________________________________________________________
554 
555 //----------------------------------------------------------------------
557  AliJetContainer *JetCont = GetJetContainer(JetContNb);
558  AliEmcalJet *SubJet=NULL;
559  Double_t DeltaR1=0;
560  Double_t DeltaR2=0;
561  AliVParticle *JetParticle=0x0;
562  Double_t SubJetiness_Numerator = 0;
563  Double_t SubJetiness_Denominator = 0;
564  Double_t Index=-2;
565  if (Reclusterer->GetNumberOfJets() < 1) return -2;
566  Index=SubJetOrdering(Jet,Reclusterer,1,0,kTRUE);
567  if(Index==-999) return -2;
568  SubJetiness_Numerator=(Reclusterer->GetJet(Index)->Pt());
569  SubJetiness_Denominator=Jet->Pt();
570  return SubJetiness_Numerator/SubJetiness_Denominator;
571 
572 
573 }
574 //__________________________________________________________________________________
575 
576 
577 
578 
579 
581 
582  AliJetContainer *jetCont = GetJetContainer(jetContNb);
583  if (!jet->GetNumberOfTracks())
584  return 0;
585  Double_t den=0.;
586  Double_t num = 0.;
587  AliVParticle *vp1 = 0x0;
588  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
589  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
590 
591  if (!vp1){
592  Printf("AliVParticle associated to constituent not found");
593  continue;
594  }
595 
596  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
597  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
598  Double_t dr = TMath::Sqrt(dr2);
599  if(dr<=fSubjetRadius) num=num+vp1->Pt();
600 
601  }
602  return num/jet->Pt();
603 }
604 
605 
606 
607 
608 //________________________________________________________________________
610  //calc subtracted jet mass
611 
612  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
615  else
616  return CoreFrac(jet, jetContNb);
617 
618 }
619 
620 
621 
622 
623 //----------------------------------------------------------------------
625  AliEmcalJet *SubJet=NULL;
626  Double_t SortingVariable;
627  Int_t ArraySize =N+1;
628  TArrayD *JetSorter = new TArrayD(ArraySize);
629  TArrayD *JetIndexSorter = new TArrayD(ArraySize);
630  for (Int_t i=0; i<ArraySize; i++){
631  JetSorter->SetAt(0,i);
632  }
633  for (Int_t i=0; i<ArraySize; i++){
634  JetIndexSorter->SetAt(0,i);
635  }
636  if(Reclusterer->GetNumberOfJets()<N) return -999;
637  for (Int_t i=0; i<Reclusterer->GetNumberOfJets(); i++){
638  SubJet=Reclusterer->GetJet(i);
639  if (Type==0) SortingVariable=SubJet->Pt();
640  else if (Type==1) SortingVariable=SubJet->E();
641  else if (Type==2) SortingVariable=SubJet->M();
642  for (Int_t j=0; j<N; j++){
643  if (SortingVariable>JetSorter->GetAt(j)){
644  for (Int_t k=N-1; k>=j; k--){
645  JetSorter->SetAt(JetSorter->GetAt(k),k+1);
646  JetIndexSorter->SetAt(JetIndexSorter->GetAt(k),k+1);
647  }
648  JetSorter->SetAt(SortingVariable,j);
649  JetIndexSorter->SetAt(i,j);
650  break;
651  }
652  }
653  }
654  if (!Index) return JetSorter->GetAt(N-1);
655  else return JetIndexSorter->GetAt(N-1);
656 }
657 
658 
659 
660 //returns -1 if the Nth hardest jet is requested where N>number of available jets
661 //type: 0=Pt 1=E 2=M
662 //Index TRUE=returns index FALSE=returns value of quantatiy in question
663 
664 
665 
666 
667 
668 
669 
670 
671 
672 
673 
674 
675 //________________________________________________________________________
677 
679  TClonesArray *tracksArray = partCont->GetArray();
680 
681  if(!partCont || !tracksArray) return -99999;
682  AliAODTrack *track = 0x0;
683  AliEmcalParticle *emcPart = 0x0;
684 
685 
686  TList *trackList = new TList();
687  Int_t triggers[100];
688  for (Int_t iTrigger=0; iTrigger<100; iTrigger++) triggers[iTrigger] = 0;
689  Int_t iTT = 0;
690 
691  for(Int_t iTrack=0; iTrack <= tracksArray->GetEntriesFast(); iTrack++){
692 
693 
694  if (fJetShapeSub == kConstSub){
695  emcPart = static_cast<AliEmcalParticle*>(tracksArray->At(iTrack));
696  if (!emcPart) continue;
697  if(TMath::Abs(emcPart->Eta())>0.9) continue;
698  if (emcPart->Pt()<0.15) continue;
699 
700  if ((emcPart->Pt() >= minpT) && (emcPart->Pt()< maxpT)) {
701  trackList->Add(emcPart);
702  triggers[iTT] = iTrack;
703  iTT++;
704  }
705  }
706  else{
707  track = static_cast<AliAODTrack*>(tracksArray->At(iTrack));
708  if (!track) continue;
709  if(TMath::Abs(track->Eta())>0.9) continue;
710  if (track->Pt()<0.15) continue;
711  if (!(track->TestFilterBit(768))) continue;
712 
713  if ((track->Pt() >= minpT) && (track->Pt()< maxpT)) {
714  trackList->Add(track);
715  triggers[iTT] = iTrack;
716  iTT++;
717 
718  }
719  }
720  }
721 
722  if (iTT == 0) return -99999;
723  Int_t nbRn = 0, index = 0 ;
724  TRandom3* random = new TRandom3(0);
725  nbRn = random->Integer(iTT);
726 
727  index = triggers[nbRn];
728  //Printf("iTT Total= %d, nbRn = %d, Index = %d",iTT, nbRn, index );
729  return index;
730 
731 }
732 
733 //----------------------------------------------------------------------
735 AliEmcalJetFinder *AliAnalysisTaskFakeJets::Recluster(AliEmcalJet *Jet, Int_t JetContNb, Double_t SubJetRadius, Double_t SubJetMinPt, Int_t Algorithm, const char* Name){
736 
737  AliJetContainer *JetCont = GetJetContainer(JetContNb);
738  AliEmcalJetFinder *Reclusterer = new AliEmcalJetFinder(Name); //JetFinder Object for reclustered jets
739  Reclusterer->SetRadius(SubJetRadius);
740  Reclusterer->SetJetMinPt(SubJetMinPt);
741  Reclusterer->SetJetAlgorithm(Algorithm); //0 for anti-kt 1 for kt
742  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
743  Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
744  if(Reclusterer->AliEmcalJetFinder::Filter(Jet, JetCont, dVtx)){;} //reclustering jet1 using the jetfinderobject Reclusterer
745  return Reclusterer;
746 }
747 
748 
749 
750 
751 //__________________________________________________________________________________
753 
754  if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
755  else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
756  if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
757  else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
758  double dphi = mphi-vphi;
759  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
760  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
761  return dphi;//dphi in [-Pi, Pi]
762 }
763 
764 
765 //________________________________________________________________________
767  //
768  // retrieve event objects
769  //
771  return kFALSE;
772 
773  return kTRUE;
774 }
775 
776 //_______________________________________________________________________
778 {
779  // Called once at the end of the analysis.
780 
781  // fTreeObservableTagging = dynamic_cast<TTree*>(GetOutputData(1));
782  // if (!fTreeObservableTagging){
783  // Printf("ERROR: fTreeObservableTagging not available");
784  // return;
785  // }
786 
787 }
788 
Float_t GetJetSecHardTrack(AliEmcalJet *jet, Int_t jetContNb)
void SetRadius(Double_t val)
Double_t GetFirstOrderSubtractedAngularity() const
double Double_t
Definition: External.C:58
void SetJetMinPt(Double_t val)
Float_t GetJetCoreFrac(AliEmcalJet *jet, Int_t jetContNb)
Definition: External.C:236
AliJetContainer * GetJetContainer(Int_t i=0) const
Double_t Eta() const
Definition: AliEmcalJet.h:108
Double_t Phi() const
Definition: AliEmcalJet.h:104
AliEmcalJetFinder * Recluster(AliEmcalJet *Jet, Int_t JetContNb, Double_t SubJetRadius, Double_t SubJetMinPt, Int_t Algorithm, const char *Name)
Int_t SelectTrigger(Float_t minpT, Float_t maxpT)
void Terminate(Option_t *option)
ClassImp(AliAnalysisTaskFakeJets) AliAnalysisTaskFakeJets
Double_t E() const
Definition: AliEmcalJet.h:106
Double_t NSubJettiness(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t A, Int_t B)
Float_t GetJetHardTrack(AliEmcalJet *jet, Int_t jetContNb)
Float_t CoreFrac(AliEmcalJet *jet, Int_t jetContNb)
Container for particles within the EMCAL framework.
Float_t Angularity(AliEmcalJet *jet, Int_t jetContNb)
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:126
Float_t HardTrack(AliEmcalJet *jet, Int_t jetContNb)
TString kData
Declare data MC or deltaAOD.
Double_t GetSubjetFraction(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, AliEmcalJetFinder *Reclusterer)
AliParticleContainer * GetParticleContainer(Int_t i=0) const
AliParticleContainer * GetParticleContainer() const
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 GetJetpTD(AliEmcalJet *jet, Int_t jetContNb)
float Float_t
Definition: External.C:68
Double_t GetSecondOrderSubtractedAngularity() const
Float_t AngularitySquared(AliEmcalJet *jet, Int_t jetContNb)
Double_t fCent
!event centrality
Float_t GetJetMass(AliEmcalJet *jet, Int_t jetContNb)
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)
AliEmcalJet * GetJet(Int_t index)
Double_t Pt() const
Definition: AliEmcalJet.h:96
Float_t GetJetAngularity(AliEmcalJet *jet, Int_t jetContNb)
AliEmcalList * fOutput
!output list
Float_t PTD(AliEmcalJet *jet, Int_t jetContNb)
Short_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:147
Double_t SubJetOrdering(AliEmcalJet *Jet, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t Type, Bool_t Index)
const Int_t nVar
void SetMakeGeneralHistograms(Bool_t g)
Float_t GetJetAngularitySquared(AliEmcalJet *jet, Int_t jetContNb)
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
Float_t GetJetNumberOfConstituents(AliEmcalJet *jet, Int_t jetContNb)
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
Float_t SecHardTrack(AliEmcalJet *jet, Int_t jetContNb)
AliEmcalJetShapeProperties * GetShapeProperties() const
Definition: AliEmcalJet.h:247
Double_t RelativePhi(Double_t mphi, Double_t vphi)
Double_t M() const
Definition: AliEmcalJet.h:107
Container for jet within the EMCAL jet framework.
std::vector< int > GetPtSortedTrackConstituentIndexes(TClonesArray *tracks) const