AliPhysics  6cf2591 (6cf2591)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskRecoilJetYield.cxx
Go to the documentation of this file.
1 //
2 // My analysis task
3 //
4 // Author: Harry Andrews
5 
6 #include <TClonesArray.h>
7 #include <TH1F.h>
8 #include <TH2F.h>
9 #include <TH3F.h>
10 #include <TCanvas.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 <AliAnalysisDataSlot.h>
21 #include <AliAnalysisDataContainer.h>
22 #include "TMatrixD.h"
23 #include "TMatrixDSym.h"
24 #include "TMatrixDSymEigen.h"
25 #include "TVector3.h"
26 #include "TVector2.h"
27 #include "AliVCluster.h"
28 #include "AliVTrack.h"
29 #include "AliEmcalJet.h"
30 #include "AliRhoParameter.h"
31 #include "AliLog.h"
32 #include "AliEmcalParticle.h"
33 #include "AliMCEvent.h"
34 #include "AliGenPythiaEventHeader.h"
35 #include "AliAODMCHeader.h"
36 #include "AliMCEvent.h"
37 #include "AliAnalysisManager.h"
38 #include "AliJetContainer.h"
39 #include "AliParticleContainer.h"
40 //#include "AliPythiaInfo.h"
41 #include "TRandom3.h"
42 #include "AliPicoTrack.h"
43 #include "AliEmcalJetFinder.h"
44 #include "AliAODEvent.h"
46 #include "AliMultSelection.h"
47 
48 //Globals
49 
50 
51 
55 Double_t Phi_Upper=2*(TMath::Pi());
56 Double_t Phi_Lower=(-1*(TMath::Pi()));
59 
60 
61 
62 
63 using std::cout;
64 using std::endl;
65 
67 
68 //________________________________________________________________________
71  fContainer(0),
72  fMinFractionShared(0),
73  fJetShapeType(kData),
74  fJetShapeSub(kNoSub),
75  fJetSelection(kInclusive),
76  fPtThreshold(-9999.),
77  fRMatching(0.2),
78  fPtMinTriggerHadron(20.),
79  fPtMaxTriggerHadron(50.),
80  fRecoilAngularWindow(0.6),
81  fSemigoodCorrect(0),
82  fHolePos(0),
83  fHoleWidth(0),
84  fCentSelectOn(kTRUE),
85  fCentMin(0),
86  fCentMax(10),
87  fSubJetAlgorithm(0),
88  fSubJetRadius(0.1),
89  fSubJetMinPt(1),
90  fJetRadius(0.4),
91  fRMatched(0.2),
92  fSharedFractionPtMin(0.5),
93  fDerivSubtrOrder(0),
94  fFullTree(kFALSE),
95  fBeta_SD(0),
96  fZCut(0.1),
97  fRho(1e-6),
98  fRhoParam(0),
99  fNsubMeasure(kFALSE),
100  fDoSoftDrop(kFALSE),
101  fhJetPt(0x0),
102  fhJetPhi(0x0),
103  fhJetEta(0x0),
104  fhJetMass(0x0),
105  fhJetRadius(0x0),
106  fhJetCounter(0x0),
107  fhNumberOfJetTracks(0x0),
108  fTreeJetInfo(0),
109  fhJetArea(0x0),
110  fhTrackPt(0x0),
111  fhPtTriggerHadron(0x0),
112  fh2PtTriggerHadronJet(0x0),
113  fhPhiTriggerHadronJet(0x0),
114  fhPhiTriggerHadronEventPlane(0x0),
115  fhPhiTriggerHadronEventPlaneTPC(0x0)
116 
117 {
118  for(Int_t i=0;i<nBranch;i++){
119  fJetInfoVar[i]=0;
120  }
121  SetMakeGeneralHistograms(kTRUE);
122  DefineOutput(1, TList::Class());
123  DefineOutput(2, TTree::Class());
124 }
125 
126 //________________________________________________________________________
128  AliAnalysisTaskEmcalJet(name, kTRUE),
129  fContainer(0),
130  fMinFractionShared(0),
131  fJetShapeType(kData),
132  fJetShapeSub(kNoSub),
133  fJetSelection(kInclusive),
134  fPtThreshold(-9999.),
135  fRMatching(0.2),
136  fPtMinTriggerHadron(20.),
137  fPtMaxTriggerHadron(50.),
138  fRecoilAngularWindow(0.6),
139  fSemigoodCorrect(0),
140  fHolePos(0),
141  fHoleWidth(0),
142  fCentSelectOn(kTRUE),
143  fCentMin(0),
144  fCentMax(10),
145  fSubJetAlgorithm(0),
146  fSubJetRadius(0.1),
147  fSubJetMinPt(1),
148  fJetRadius(0.4),
149  fRMatched(0.2),
150  fSharedFractionPtMin(0.5),
151  fDerivSubtrOrder(0),
152  fFullTree(kFALSE),
153  fBeta_SD(0),
154  fZCut(0.1),
155  fRho(1e-6),
156  fRhoParam(0),
157  fNsubMeasure(kFALSE),
158  fDoSoftDrop(kFALSE),
159  fhJetPt(0x0),
160  fhJetPhi(0x0),
161  fhJetEta(0x0),
162  fhJetMass(0x0),
163  fhJetRadius(0x0),
164  fhJetCounter(0x0),
165  fhNumberOfJetTracks(0x0),
166  fTreeJetInfo(0),
167  fhJetArea(0x0),
168  fhTrackPt(0x0),
169  fhPtTriggerHadron(0x0),
170  fh2PtTriggerHadronJet(0x0),
171  fhPhiTriggerHadronJet(0x0),
172  fhPhiTriggerHadronEventPlane(0x0),
173  fhPhiTriggerHadronEventPlaneTPC(0x0)
174 
175 {
176  // Standard constructor.
177  for(Int_t i=0;i<nBranch;i++){
178  fJetInfoVar[i]=0;
179  }
181  DefineOutput(1, TList::Class());
182  DefineOutput(2, TTree::Class());
183 }
184 
185 //________________________________________________________________________
187 {
188  // Destructor.
189 }
190 
191 //________________________________________________________________________
193 {
194  // Create user output.
195 
197 
198  Bool_t oldStatus = TH1::AddDirectoryStatus();
199  TH1::AddDirectory(kFALSE);
200  TH1::AddDirectory(oldStatus);
201  const char* nameoutput = GetOutputSlot(2)->GetContainer()->GetName();
202  fTreeJetInfo = new TTree(nameoutput, nameoutput);
203 
204  if (!fFullTree){
205  const Int_t nVarMin = 13;
206  TString *fJetInfoVarNames = new TString [nVarMin];
207 
208  fJetInfoVarNames[0] = "Pt";
209  fJetInfoVarNames[1] = "Phi";
210  fJetInfoVarNames[2] = "Eta";
211  fJetInfoVarNames[3] = "SymParam";
212  fJetInfoVarNames[4] = "Tau1";
213  fJetInfoVarNames[5] = "Tau2";
214  fJetInfoVarNames[6] = "Mass";
215  fJetInfoVarNames[7] = "NTracks";
216  fJetInfoVarNames[8] = "JetPtNoCut";
217  fJetInfoVarNames[9] = "PTD";
218  fJetInfoVarNames[10] = "Angularity";
219  fJetInfoVarNames[11] = "Centrality";
220  fJetInfoVarNames[12] = "DelR";
221 
222 
223  for(Int_t ivar=0; ivar < nVarMin; ivar++){
224  cout<<"looping over variables"<<endl;
225  fTreeJetInfo->Branch(fJetInfoVarNames[ivar].Data(), &fJetInfoVar[ivar], Form("%s/D", fJetInfoVarNames[ivar].Data()));
226  }
227  }
228 
230 
231  fhJetPt= new TH1F("fhJetPt", "Jet Pt",1500,-0.5,149.5 );
232  fOutput->Add(fhJetPt);
233  fhJetPhi= new TH1F("fhJetPhi", "Jet Phi",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
234  fOutput->Add(fhJetPhi);
235  fhJetEta= new TH1F("fhJetEta", "Jet Eta", Eta_Bins_1, Eta_Lower, Eta_Upper);
236  fOutput->Add(fhJetEta);
237  fhJetMass= new TH1F("fhJetMass", "Jet Mass", 4000,-0.5, 39.5);
238  fOutput->Add(fhJetMass);
239  fhJetRadius= new TH1F("fhJetRadius", "Jet Radius", 100, -0.05,0.995);
240  fOutput->Add(fhJetRadius);
241  fhNumberOfJetTracks= new TH1F("fhNumberOfJetTracks", "Number of Tracks within a Jet", 300, -0.5,299.5);
243  fhJetCounter= new TH1F("fhJetCounter", "Jet Counter", 150, -0.5, 149.5);
244  fOutput->Add(fhJetCounter);
245  fhJetArea= new TH1F("fhJetArea", "Jet Area", 400,-0.5, 1.95);
246  fOutput->Add(fhJetArea);
247  fhTrackPt= new TH1F("fhTrackPt", "Track Pt",600,-0.5,59.5);
248  fOutput->Add(fhTrackPt);
249 
250  if(fJetSelection == kRecoil){
251  fhPtTriggerHadron= new TH1F("fhPtTriggerHadron", "fhPtTriggerHadron",1500,-0.5,149.5);
253  fh2PtTriggerHadronJet= new TH2F("fh2PtTriggerHadronJet", "fh2PtTriggerHadronJet",1500,-0.5,149.5,1500,-0.5,149.5);
255  fhPhiTriggerHadronJet= new TH1F("fhPhiTriggerHadronJet", "fhPhiTriggerHadronJet",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
257  fhPhiTriggerHadronEventPlane= new TH1F("fhPhiTriggerHadronEventPlane", "fhPhiTriggerHadronEventPlane",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
259  fhPhiTriggerHadronEventPlaneTPC= new TH1F("fhPhiTriggerHadronEventPlaneTPC", "fhPhiTriggerHadronEventPlaneTPC",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
261 
262 
263  }
264  }
265  PostData(1,fOutput);
266  PostData(2,fTreeJetInfo);
267  // delete [] fShapesVarNames;
268 }
269 
270 //________________________________________________________________________
272 {
273  // Run analysis code here, if needed. It will be executed before FillHistograms().
274 
275  return kTRUE;
276 }
277 
278 //________________________________________________________________________
280 {
281  AliMultSelection *MultSelection = 0x0;
282  AliAODEvent *fEvent = dynamic_cast<AliAODEvent*>(InputEvent());
283  if (!fEvent) {
284  Error("UserExec","AOD not available");
285  return 0;
286  }
287  MultSelection = (AliMultSelection * ) fEvent->FindListObject("MultSelection");
288  if(!MultSelection) {
289  AliWarning("AliMultSelection object not found!");
290  }
291  else{
292  Percentile_1 = MultSelection->GetMultiplicityPercentile("V0M");
293  }
294  if (fCentSelectOn){
295  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
296  }
297 
298  AliAODTrack *TriggerHadron = 0x0;
299  if (fJetSelection == kRecoil) {
300  //you have to set a flag and the limits of the pT interval for your trigger
302  if (TriggerHadronLabel==-99999) return 0; //Trigger Hadron Not Found
303  AliTrackContainer *PartCont =NULL;
304  if (fJetShapeSub==kConstSub) PartCont = GetTrackContainer("tracksConstSub");
305  else PartCont = GetTrackContainer("tracks");
306  TClonesArray *TrackArray = PartCont->GetArray();
307  TriggerHadron = static_cast<AliAODTrack*>(TrackArray->At(TriggerHadronLabel));
308  if (!TriggerHadron) return 0;//No trigger hadron with label found
309  if(fSemigoodCorrect){
310  Double_t HoleDistance=RelativePhi(TriggerHadron->Phi(),fHolePos);
311  if(TMath::Abs(HoleDistance)+fHoleWidth+fJetRadius>TMath::Pi()-fRecoilAngularWindow) return 0;
312  }
313  fhPtTriggerHadron->Fill(TriggerHadron->Pt()); //Needed for per trigger Normalisation
314  fhPhiTriggerHadronEventPlane->Fill(TMath::Abs(RelativePhiEventPlane(fEPV0,TriggerHadron->Phi()))); //fEPV0 is the event plane from AliAnalysisTaskEmcal
315  fhPhiTriggerHadronEventPlaneTPC->Fill(TMath::Abs(RelativePhiEventPlane(((AliVAODHeader*)(InputEvent()->GetHeader()))->GetEventplane(),TriggerHadron->Phi()))); //TPC event plane
316  }
318  AliEmcalJet *Jet1 = NULL; //Original Jet in the event
319  AliJetContainer *JetCont= GetJetContainer(0); //Jet Container for event
320  Int_t JetCounter=0; //Counts number of jets in event
321  Double_t JetPhi=0;
322  Bool_t EventCounter=kFALSE;
323  Double_t JetPt_ForThreshold=0;
324  AliEmcalJet *Jet2= NULL;
325  if(JetCont) {
326  JetCont->ResetCurrentID();
327  while((Jet1=JetCont->GetNextAcceptJet())) {
328  if(!Jet1) continue;
329  fJetInfoVar[8]=Jet1->Pt();
330  if (fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) JetPt_ForThreshold = Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
331  else JetPt_ForThreshold = Jet1->Pt();
332  if(JetPt_ForThreshold<fPtThreshold) {
333  continue;
334  }
335  else {
336  Float_t RecoilDeltaPhi = 0.;
337  if (fJetSelection == kRecoil){
338  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
339  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
340  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
341  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
342  }
343  if (!EventCounter){
344  EventCounter=kTRUE;
345  }
346  JetCounter++;
347  fhJetPt->Fill(Jet1->Pt());
348  fhJetArea->Fill(Jet1->Area());
349  JetPhi=Jet1->Phi();
350  if(JetPhi < -1*TMath::Pi()) JetPhi += (2*TMath::Pi());
351  else if (JetPhi > TMath::Pi()) JetPhi -= (2*TMath::Pi());
352  fhJetPhi->Fill(JetPhi);
353  fhJetEta->Fill(Jet1->Eta());
354  fhJetMass->Fill(Jet1->M());
355  fhJetRadius->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
356  fhNumberOfJetTracks->Fill(Jet1->GetNumberOfTracks());
357  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fJetInfoVar[0]= Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
358  else fJetInfoVar[0]=Jet1->Pt();
359  fJetInfoVar[1]=JetPhi;
360  fJetInfoVar[2]=Jet1->Eta();
361  if(fDoSoftDrop) SoftDrop(Jet1,JetCont,fZCut,fBeta_SD);
362  else{
363  fJetInfoVar[3]=0;
364  fJetInfoVar[4]=0;
365  fJetInfoVar[5]=0;
366  fJetInfoVar[12]=0;
367  }
368  fJetInfoVar[6]=Jet1->M();
370  fJetInfoVar[9]=PTD(Jet1,0);
371  fJetInfoVar[10]=Angularity(Jet1,0);
372  fJetInfoVar[11]=fCent;
373  fTreeJetInfo->Fill();
374  }
375  }
376  fhJetCounter->Fill(JetCounter); //Number of Jets in Each Event
377  }
378  }
379 
380  return kTRUE;
381 }
382 
383 //________________________________________________________________________
385 
386  if(Phi < -1*TMath::Pi()) Phi += (2*TMath::Pi());
387  else if (Phi > TMath::Pi()) Phi -= (2*TMath::Pi());
388  Double_t DeltaPhi=Phi-EventPlane;
389  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
390  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
391  return DeltaPhi;
392 }
393 //________________________________________________________________________
395 
396  if(Phi1 < -1*TMath::Pi()) Phi1 += (2*TMath::Pi());
397  else if (Phi1 > TMath::Pi()) Phi1 -= (2*TMath::Pi());
398  if(Phi2 < -1*TMath::Pi()) Phi2 += (2*TMath::Pi());
399  else if (Phi2 > TMath::Pi()) Phi2 -= (2*TMath::Pi());
400  Double_t DeltaPhi=Phi2-Phi1;
401  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
402  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
403  return DeltaPhi;
404 }
405 
406 
407 
408 
409 //--------------------------------------------------------------------------
411 
412  AliJetContainer *JetCont = GetJetContainer(JetContNb);
413  Double_t Angularity_Numerator=0;
414  Double_t Angularity_Denominator=0;
415  AliVParticle *Particle=0x0;
416  Double_t DeltaPhi=0;
417 
418 
419  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks
420  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
421  if(!Particle) continue;
422  DeltaPhi=RelativePhi(Jet->Phi(),Particle->Phi());
423  Angularity_Numerator=Angularity_Numerator+(Particle->Pt()*TMath::Sqrt(((Particle->Eta()-Jet->Eta())*(Particle->Eta()-Jet->Eta()))+(DeltaPhi*DeltaPhi)));
424  Angularity_Denominator= Angularity_Denominator+Particle->Pt();
425  }
426  if(Angularity_Denominator!=0) return Angularity_Numerator/Angularity_Denominator;
427  else return -1;
428 
429 }
430 
431 
432 
433 //--------------------------------------------------------------------------
435 
436  AliJetContainer *JetCont = GetJetContainer(JetContNb);
437  Double_t PTD_Numerator=0; //Reset these values
438  Double_t PTD_Denominator=0;
439  AliVParticle *Particle=0x0;
440  Double_t DeltaPhi=0;
441  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks
442  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
443  if(!Particle) continue;
444  PTD_Numerator=PTD_Numerator+(Particle->Pt()*Particle->Pt());
445  PTD_Denominator=PTD_Denominator+Particle->Pt();
446  }
447  if(PTD_Denominator!=0) return TMath::Sqrt(PTD_Numerator)/PTD_Denominator;
448  else return -1;
449 
450 }
451 
452 //_________________________________________________________________________
453 void AliAnalysisTaskRecoilJetYield::SoftDrop(AliEmcalJet *fJet,AliJetContainer *fJetCont, double zcut, double beta){
454  std::vector<fastjet::PseudoJet> fInputVectors;
455  Double_t JetInvMass=0, PseudJetInvMass=0, TrackMom = 0, TrackEnergy = 0;
456  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
457  Double_t JetEta=fJet->Eta(),JetPhi=fJet->Phi();
458  Double_t FJTrackEta[9999],FJTrackPhi[9999],FJTrackPt[9999],EmcalJetTrackEta[9999],EmcalJetTrackPhi[9999],EmcalJetTrackPt[9999];
459  UShort_t FJNTracks=0,EmcalJetNTracks=0;
460  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
461  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
462  JetInvMass += fTrk->M();
463  if (!fTrk) continue;
464  fastjet::PseudoJet PseudoTracks(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
465  TrackMom += TMath::Sqrt(TMath::Power(fTrk->Px(),2)+TMath::Power(fTrk->Py(),2)+TMath::Power(fTrk->Pz(),2));
466  TrackEnergy += fTrk->E();
467  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
468  PseudJetInvMass += PseudoTracks.m();
469  fInputVectors.push_back(PseudoTracks);
470  EmcalJetTrackEta[i]=fTrk->Eta();
471  EmcalJetTrackPhi[i]=fTrk->Phi();
472  EmcalJetTrackPt[i]=fTrk->Pt();
473  EmcalJetNTracks++;
474  }
475  fastjet::JetDefinition *fJetDef;
476  fastjet::ClusterSequence *fClustSeqSA;
477 
478 
479 
480  fJetDef = new fastjet::JetDefinition(fastjet::antikt_algorithm, fJetRadius*2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
481 
482  try {
483  fClustSeqSA = new fastjet::ClusterSequence(fInputVectors, *fJetDef);
484  } catch (fastjet::Error) {
485  AliError(" [w] FJ Exception caught.");
486  //return -1;
487  }
488 
489  std::vector<fastjet::PseudoJet> fOutputJets;
490  fOutputJets.clear();
491  fOutputJets=fClustSeqSA->inclusive_jets(0);
492 
493 
494  std::vector<fastjet::PseudoJet> jet_constituents = fOutputJets[0].constituents();
495  Float_t NSubjettinessResult1, NSubjettinessResult2, NSubBeta = 1, R0=0.4;
496  std::vector<fastjet::PseudoJet> Subjet_Axes;
497  fastjet::PseudoJet SubJet1_Axis,SubJet2_Axis;
498  if(jet_constituents.size()>=2){
499  fastjet::contrib::Nsubjettiness nSub1(1,fastjet::contrib::KT_Axes(),fastjet::contrib::NormalizedMeasure(NSubBeta,R0));
500  fastjet::contrib::Nsubjettiness nSub2(2,fastjet::contrib::KT_Axes(),fastjet::contrib::NormalizedMeasure(NSubBeta,R0));
501 
502  NSubjettinessResult1 = nSub1.result(fOutputJets[0]);
503  NSubjettinessResult2 = nSub2.result(fOutputJets[0]);
504  Subjet_Axes = nSub2.currentAxes();
505  SubJet1_Axis = Subjet_Axes[0];
506  SubJet2_Axis = Subjet_Axes[1];
507 
508  Double_t SubJet1_Eta=SubJet1_Axis.pseudorapidity();
509  Double_t SubJet2_Eta=SubJet2_Axis.pseudorapidity();
510  Double_t SubJet1_Phi=SubJet1_Axis.phi();
511  if(SubJet1_Phi < -1*TMath::Pi()) SubJet1_Phi += (2*TMath::Pi());
512  else if (SubJet1_Phi > TMath::Pi()) SubJet1_Phi -= (2*TMath::Pi());
513  Double_t SubJet2_Phi=SubJet2_Axis.phi();
514  if(SubJet1_Phi < -1*TMath::Pi()) SubJet1_Phi += (2*TMath::Pi());
515  else if (SubJet1_Phi > TMath::Pi()) SubJet1_Phi -= (2*TMath::Pi());
516 
517  Double_t DeltaPhiSubJets,DeltaEtaSubJets;
518  DeltaPhiSubJets=SubJet1_Phi-SubJet2_Phi;
519  DeltaEtaSubJets=SubJet1_Eta-SubJet2_Eta;
520  if(DeltaPhiSubJets < -1*TMath::Pi()) DeltaPhiSubJets += (2*TMath::Pi());
521  else if (DeltaPhiSubJets > TMath::Pi()) DeltaPhiSubJets -= (2*TMath::Pi());
522  Double_t DelR=0;
523  DelR = TMath::Power(TMath::Power(DeltaPhiSubJets,2)+TMath::Power(DeltaEtaSubJets,2),0.5);
524 
525  fJetInfoVar[12]=DelR;
526  }
527 
528  fastjet::contrib::SoftDrop softdrop(beta, zcut);
529  fastjet::PseudoJet finaljet = softdrop(fOutputJets[0]);
530  fJetInfoVar[4]=NSubjettinessResult1;
531  fJetInfoVar[5]=NSubjettinessResult2;
532  AliEmcalJet* jet = new AliEmcalJet(finaljet.perp(), finaljet.eta(), finaljet.phi(), finaljet.m());
533  std::vector<fastjet::PseudoJet> fSDTracks=finaljet.constituents();
534  Double_t FastjetTrackDelR,EmcalTrackDelR;
535  for(Int_t i=0;i<fJet->GetNumberOfConstituents();i++){
536  if(i<=finaljet.constituents().size()){
537  FastjetTrackDelR = TMath::Sqrt(TMath::Power(fSDTracks[i].eta()-JetEta,2)+TMath::Power(fSDTracks[i].phi()-JetPhi,2));
538  FJTrackEta[i]=fSDTracks[i].eta();
539  FJTrackPhi[i]=fSDTracks[i].phi();
540  FJTrackPt[i]=fSDTracks[i].perp();
541  FJNTracks++;
542  }
543  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
544  EmcalTrackDelR = TMath::Sqrt(TMath::Power(fTrk->Eta()-JetEta,2)+TMath::Power(fTrk->Phi()-JetPhi,2));
545  }
546  Int_t NDroppedTracks = fJet->GetNumberOfTracks()-finaljet.constituents().size();
547  Int_t nConstituents(fClustSeqSA->constituents(finaljet).size());
548  jet->SetNumberOfTracks(nConstituents);
549  Double_t SymParam, Mu, DeltaR;
550  SymParam=(finaljet.structure_of<fastjet::contrib::SoftDrop>().symmetry());
551  Mu=(finaljet.structure_of<fastjet::contrib::SoftDrop>().mu());
552  DeltaR=(finaljet.structure_of<fastjet::contrib::SoftDrop>().delta_R());
553  fJetInfoVar[3]=SymParam;
554  return;
555 
556 
557 }
558 
559 //________________________________________________________________________
561 
562  AliTrackContainer *PartCont = NULL;
563  if (fJetShapeSub==kConstSub) PartCont = GetTrackContainer("tracksConstSub");
564  else PartCont = GetTrackContainer("tracks");
565  TClonesArray *TracksArray = PartCont->GetArray();
566  if(!PartCont || !TracksArray) return -99999;
567  AliAODTrack *Track = 0x0;
568  Int_t Trigger_Index[100];
569  for (Int_t i=0; i<100; i++) Trigger_Index[i] = 0;
570  Int_t Trigger_Counter = 0;
571  for(Int_t i=0; i < TracksArray->GetEntriesFast(); i++){
572  if((Track = static_cast<AliAODTrack*>(PartCont->GetAcceptTrack(i)))){
573  if (!Track) continue;
574  fhTrackPt->Fill(Track->Pt());
575  if(TMath::Abs(Track->Eta())>0.9) continue;
576  if (Track->Pt()<0.15) continue;
577  if ((Track->Pt() >= PtMin) && (Track->Pt()< PtMax)) {
578  Trigger_Index[Trigger_Counter] = i;
579  Trigger_Counter++;
580  }
581  }
582  }
583  if (Trigger_Counter == 0) return -99999;
584  Int_t RandomNumber = 0, Index = 0 ;
585  TRandom3* Random = new TRandom3(0);
586  RandomNumber = Random->Integer(Trigger_Counter);
587  Index = Trigger_Index[RandomNumber];
588  return Index;
589 }
590 
591 //________________________________________________________________________
593  //
594  // retrieve event objects
595  //
597  return kFALSE;
598 
599  return kTRUE;
600 }
601 
602 
603 //_______________________________________________________________________
605 {
606  // Called once at the end of the analysis.
607 
608 
609 }
Double_t Area() const
Definition: AliEmcalJet.h:123
double Double_t
Definition: External.C:58
Double_t PTD(AliEmcalJet *Jet, Int_t JetContNb)
static Double_t DeltaPhi(Double_t phia, Double_t phib, Double_t rMin=-TMath::Pi()/2, Double_t rMax=3 *TMath::Pi()/2)
Definition: External.C:236
AliJetContainer * GetJetContainer(Int_t i=0) const
Double_t Eta() const
Definition: AliEmcalJet.h:114
Double_t Phi_Upper
Double_t Phi() const
Definition: AliEmcalJet.h:110
ClassImp(AliAnalysisTaskRecoilJetYield) AliAnalysisTaskRecoilJetYield
Container with name, TClonesArray and cuts for particles.
Double_t fEPV0
!event plane V0
Double_t Eta_Lower
Float_t Percentile_1
Double_t Eta_Upper
Double_t Phi_Lower
UShort_t GetNumberOfConstituents() const
Definition: AliEmcalJet.h:133
Container for particles within the EMCAL framework.
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:132
TString kData
Declare data MC or deltaAOD.
AliParticleContainer * GetParticleContainer() const
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
Int_t SelectTriggerHadron(Float_t PtMin, Float_t PtMax)
Double_t fCent
!event centrality
const Int_t nBranch
Double_t Angularity(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)
Double_t RelativePhi(Double_t Phi1, Double_t Phi2)
Double_t Pt() const
Definition: AliEmcalJet.h:102
Double_t GetRhoVal(Int_t i=0) const
AliEmcalList * fOutput
!output list
Short_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:153
AliTrackContainer * GetTrackContainer(Int_t i=0) const
Double_t DeltaR(const AliVParticle *part1, const AliVParticle *part2)
void SetMakeGeneralHistograms(Bool_t g)
virtual AliVTrack * GetAcceptTrack(Int_t i=-1) const
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
void SetNumberOfTracks(Int_t n)
Definition: AliEmcalJet.h:192
unsigned short UShort_t
Definition: External.C:28
void SoftDrop(AliEmcalJet *fJet, AliJetContainer *fJetCont, double zcut, double beta)
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
Double_t RelativePhiEventPlane(Double_t EventPlane, Double_t Phi)
Double_t M() const
Definition: AliEmcalJet.h:113
Container for jet within the EMCAL jet framework.