AliPhysics  7f2a7c4 (7f2a7c4)
 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  fhGroomedPtvJetPt(0x0),
112  fhDroppedBranches(0x0),
113  fhPtTriggerHadron(0x0),
114  fh2PtTriggerHadronJet(0x0),
115  fhPhiTriggerHadronJet(0x0),
116  fhPhiTriggerHadronEventPlane(0x0),
117  fhPhiTriggerHadronEventPlaneTPC(0x0)
118 
119 {
120  for(Int_t i=0;i<nBranch;i++){
121  fJetInfoVar[i]=0;
122  }
123  SetMakeGeneralHistograms(kTRUE);
124  DefineOutput(1, TList::Class());
125  DefineOutput(2, TTree::Class());
126 }
127 
128 //________________________________________________________________________
130  AliAnalysisTaskEmcalJet(name, kTRUE),
131  fContainer(0),
132  fMinFractionShared(0),
133  fJetShapeType(kData),
134  fJetShapeSub(kNoSub),
135  fJetSelection(kInclusive),
136  fPtThreshold(-9999.),
137  fRMatching(0.2),
138  fPtMinTriggerHadron(20.),
139  fPtMaxTriggerHadron(50.),
140  fRecoilAngularWindow(0.6),
141  fSemigoodCorrect(0),
142  fHolePos(0),
143  fHoleWidth(0),
144  fCentSelectOn(kTRUE),
145  fCentMin(0),
146  fCentMax(10),
147  fSubJetAlgorithm(0),
148  fSubJetRadius(0.1),
149  fSubJetMinPt(1),
150  fJetRadius(0.4),
151  fRMatched(0.2),
152  fSharedFractionPtMin(0.5),
153  fDerivSubtrOrder(0),
154  fFullTree(kFALSE),
155  fBeta_SD(0),
156  fZCut(0.1),
157  fRho(1e-6),
158  fRhoParam(0),
159  fNsubMeasure(kFALSE),
160  fDoSoftDrop(kFALSE),
161  fhJetPt(0x0),
162  fhJetPhi(0x0),
163  fhJetEta(0x0),
164  fhJetMass(0x0),
165  fhJetRadius(0x0),
166  fhJetCounter(0x0),
167  fhNumberOfJetTracks(0x0),
168  fTreeJetInfo(0),
169  fhJetArea(0x0),
170  fhTrackPt(0x0),
171  fhGroomedPtvJetPt(0x0),
172  fhDroppedBranches(0x0),
173  fhPtTriggerHadron(0x0),
174  fh2PtTriggerHadronJet(0x0),
175  fhPhiTriggerHadronJet(0x0),
176  fhPhiTriggerHadronEventPlane(0x0),
177  fhPhiTriggerHadronEventPlaneTPC(0x0)
178 
179 {
180  // Standard constructor.
181  for(Int_t i=0;i<nBranch;i++){
182  fJetInfoVar[i]=0;
183  }
185  DefineOutput(1, TList::Class());
186  DefineOutput(2, TTree::Class());
187 }
188 
189 //________________________________________________________________________
191 {
192  // Destructor.
193 }
194 
195 //________________________________________________________________________
197 {
198  // Create user output.
199 
201 
202  Bool_t oldStatus = TH1::AddDirectoryStatus();
203  TH1::AddDirectory(kFALSE);
204  TH1::AddDirectory(oldStatus);
205  const char* nameoutput = GetOutputSlot(2)->GetContainer()->GetName();
206  fTreeJetInfo = new TTree(nameoutput, nameoutput);
207 
208  if (!fFullTree){
209  const Int_t nVarMin = 18;
210  TString *fJetInfoVarNames = new TString [nVarMin];
211 
212  fJetInfoVarNames[0] = "Pt";
213  fJetInfoVarNames[1] = "Pt_Truth";
214  fJetInfoVarNames[2] = "SymParam";
215  fJetInfoVarNames[3] = "SymParam_Truth";
216  fJetInfoVarNames[4] = "Tau1";
217  fJetInfoVarNames[5] = "Tau1_Truth";
218  fJetInfoVarNames[6] = "Tau2";
219  fJetInfoVarNames[7] = "Tau2_Truth";
220  fJetInfoVarNames[8] = "PTD";
221  fJetInfoVarNames[9] = "PTD_Truth";
222  fJetInfoVarNames[10] = "Angularity";
223  fJetInfoVarNames[11] = "Angularity_Truth";
224  fJetInfoVarNames[12] = "DelR";
225  fJetInfoVarNames[13] = "DelR_Truth";
226  fJetInfoVarNames[14] = "N_Groomed_Branches";
227  fJetInfoVarNames[15] = "N_Groomed_Branches_Truth";
228  fJetInfoVarNames[16] = "Groomed_Jet_Pt";
229  fJetInfoVarNames[17] = "Groomed_Jet_Pt_Truth";
230 
231 
232 
233  for(Int_t ivar=0; ivar < nVarMin; ivar++){
234  cout<<"looping over variables"<<endl;
235  fTreeJetInfo->Branch(fJetInfoVarNames[ivar].Data(), &fJetInfoVar[ivar], Form("%s/D", fJetInfoVarNames[ivar].Data()));
236  }
237  }
238 
240 
241  fhJetPt= new TH1F("fhJetPt", "Jet Pt",1500,-0.5,149.5 );
242  fOutput->Add(fhJetPt);
243  fhJetPhi= new TH1F("fhJetPhi", "Jet Phi",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
244  fOutput->Add(fhJetPhi);
245  fhJetEta= new TH1F("fhJetEta", "Jet Eta", Eta_Bins_1, Eta_Lower, Eta_Upper);
246  fOutput->Add(fhJetEta);
247  fhJetMass= new TH1F("fhJetMass", "Jet Mass", 4000,-0.5, 39.5);
248  fOutput->Add(fhJetMass);
249  fhJetRadius= new TH1F("fhJetRadius", "Jet Radius", 100, -0.05,0.995);
250  fOutput->Add(fhJetRadius);
251  fhNumberOfJetTracks= new TH1F("fhNumberOfJetTracks", "Number of Tracks within a Jet", 300, -0.5,299.5);
253  fhJetCounter= new TH1F("fhJetCounter", "Jet Counter", 150, -0.5, 149.5);
254  fOutput->Add(fhJetCounter);
255  fhJetArea= new TH1F("fhJetArea", "Jet Area", 400,-0.5, 1.95);
256  fOutput->Add(fhJetArea);
257  fhTrackPt= new TH1F("fhTrackPt", "Track Pt",600,-0.5,59.5);
258  fOutput->Add(fhTrackPt);
259  fhGroomedPtvJetPt= new TH2F("fhGroomedPtvJetPt","Groomed Jet p_{T} v Original Jet p_{T}",150,0,150,150,0,150);
261  fhDroppedBranches= new TH1F("fhDroppedBranches","Number of Softdropped branches",50,0,50);
263 
264  if(fJetSelection == kRecoil){
265  fhPtTriggerHadron= new TH1F("fhPtTriggerHadron", "fhPtTriggerHadron",1500,-0.5,149.5);
267  fh2PtTriggerHadronJet= new TH2F("fh2PtTriggerHadronJet", "fh2PtTriggerHadronJet",1500,-0.5,149.5,1500,-0.5,149.5);
269  fhPhiTriggerHadronJet= new TH1F("fhPhiTriggerHadronJet", "fhPhiTriggerHadronJet",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
271  fhPhiTriggerHadronEventPlane= new TH1F("fhPhiTriggerHadronEventPlane", "fhPhiTriggerHadronEventPlane",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
273  fhPhiTriggerHadronEventPlaneTPC= new TH1F("fhPhiTriggerHadronEventPlaneTPC", "fhPhiTriggerHadronEventPlaneTPC",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
275 
276 
277  }
278  }
279  PostData(1,fOutput);
280  PostData(2,fTreeJetInfo);
281  // delete [] fShapesVarNames;
282 }
283 
284 //________________________________________________________________________
286 {
287  // Run analysis code here, if needed. It will be executed before FillHistograms().
288 
289  return kTRUE;
290 }
291 
292 //________________________________________________________________________
294 {
295  AliMultSelection *MultSelection = 0x0;
296  AliAODEvent *fEvent = dynamic_cast<AliAODEvent*>(InputEvent());
297  if (!fEvent) {
298  Error("UserExec","AOD not available");
299  return 0;
300  }
301  MultSelection = (AliMultSelection * ) fEvent->FindListObject("MultSelection");
302  if(!MultSelection) {
303  AliWarning("AliMultSelection object not found!");
304  }
305  else{
306  Percentile_1 = MultSelection->GetMultiplicityPercentile("V0M");
307  }
308  if (fCentSelectOn){
309  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
310  }
311 
313  AliAODTrack *TriggerHadron = 0x0;
314  if (fJetSelection == kRecoil) {
315  //you have to set a flag and the limits of the pT interval for your trigger
317  if (TriggerHadronLabel==-99999) return 0; //Trigger Hadron Not Found
318  AliTrackContainer *PartCont =NULL;
319  if (fJetShapeSub==kConstSub) PartCont = GetTrackContainer(1);
320  else PartCont = GetTrackContainer(0);
321  TClonesArray *TrackArray = PartCont->GetArray();
322  TriggerHadron = static_cast<AliAODTrack*>(TrackArray->At(TriggerHadronLabel));
323  if (!TriggerHadron) return 0;//No trigger hadron with label found
324  if(fSemigoodCorrect){
325  Double_t HoleDistance=RelativePhi(TriggerHadron->Phi(),fHolePos);
326  if(TMath::Abs(HoleDistance)+fHoleWidth+fJetRadius>TMath::Pi()-fRecoilAngularWindow) return 0;
327  }
328  fhPtTriggerHadron->Fill(TriggerHadron->Pt()); //Needed for per trigger Normalisation
329  fhPhiTriggerHadronEventPlane->Fill(TMath::Abs(RelativePhiEventPlane(fEPV0,TriggerHadron->Phi()))); //fEPV0 is the event plane from AliAnalysisTaskEmcal
330  fhPhiTriggerHadronEventPlaneTPC->Fill(TMath::Abs(RelativePhiEventPlane(((AliVAODHeader*)(InputEvent()->GetHeader()))->GetEventplane(),TriggerHadron->Phi()))); //TPC event plane
331  }
332 
333 
336  AliEmcalJet *Jet1 = NULL; //Original Jet in the event
337  AliJetContainer *JetCont= GetJetContainer(0); //Jet Container for event
338  Int_t JetCounter=0; //Counts number of jets in event
339  Double_t JetPhi=0;
340  Bool_t EventCounter=kFALSE;
341  Double_t JetPt_ForThreshold=0;
342  AliEmcalJet *Jet2= NULL;
343  if(JetCont) {
344  JetCont->ResetCurrentID();
345  while((Jet1=JetCont->GetNextAcceptJet())) {
346  if(!Jet1) continue;
347  if (fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) JetPt_ForThreshold = Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
348  else JetPt_ForThreshold = Jet1->Pt();
349  if(JetPt_ForThreshold<fPtThreshold) {
350  continue;
351  }
352  else {
353  Float_t RecoilDeltaPhi = 0.;
354  if (fJetSelection == kRecoil){
355  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
356  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
357  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
358  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
359  }
360  if (!EventCounter){
361  EventCounter=kTRUE;
362  }
363  JetCounter++;
364  fhJetPt->Fill(Jet1->Pt());
365  fhJetArea->Fill(Jet1->Area());
366  JetPhi=Jet1->Phi();
367  if(JetPhi < -1*TMath::Pi()) JetPhi += (2*TMath::Pi());
368  else if (JetPhi > TMath::Pi()) JetPhi -= (2*TMath::Pi());
369  fhJetPhi->Fill(JetPhi);
370  fhJetEta->Fill(Jet1->Eta());
371  fhJetMass->Fill(Jet1->M());
372  fhJetRadius->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
373  fhNumberOfJetTracks->Fill(Jet1->GetNumberOfTracks());
374  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fJetInfoVar[0]= Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
375  else fJetInfoVar[0]=Jet1->Pt();
376  fJetInfoVar[1]=0;
377  if(fDoSoftDrop) {
378  SoftDrop(Jet1,JetCont,fZCut,fBeta_SD,kFALSE);
379  SoftDrop(Jet1,JetCont,fZCut,fBeta_SD,kTRUE);
380  }
381  else{
382  fJetInfoVar[2]=0;
383  fJetInfoVar[3]=0;
384  fJetInfoVar[4]=0;
385  fJetInfoVar[5]=0;
386  fJetInfoVar[6]=0;
387  fJetInfoVar[7]=0;
388  fJetInfoVar[12]=0;
389  fJetInfoVar[13]=0;
390  fJetInfoVar[14]=0;
391  fJetInfoVar[15]=0;
392  fJetInfoVar[16]=0;
393  fJetInfoVar[17]=0;
394  }
395  fJetInfoVar[8]=PTD(Jet1,0);
396  fJetInfoVar[9]=0;
397  fJetInfoVar[10]=Angularity(Jet1,0);
398  fJetInfoVar[11]=0;
399  fTreeJetInfo->Fill();
400  }
401  }
402  fhJetCounter->Fill(JetCounter); //Number of Jets in Each Event
403  }
404  }
405 
407  if(fJetShapeType == kDetEmbPart){
408  AliEmcalJet *JetHybridS = NULL; //Subtracted hybrid Jet
409  AliEmcalJet *JetHybridUS = NULL; //Unsubtracted Hybrid Jet
410  AliEmcalJet *JetPythDet = NULL; //Detector Level Pythia Jet
411  AliEmcalJet *JetPythTrue = NULL; //Particle Level Pyhtia Jet
412  AliJetContainer *JetContHybridS= GetJetContainer(0); //Jet Container for Subtracted Hybrid Jets
413  AliJetContainer *JetContHybridUS= GetJetContainer(1); //Jet Container for Unsubtracted Hybrid Jets
414  AliJetContainer *JetContPythDet= GetJetContainer(2); //Jet Container for Detector Level Pyhtia Jets
415  AliJetContainer *JetContPythTrue= GetJetContainer(3); //Jet Container for Particle Level Pythia Jets
416 
417  Bool_t JetsMatched = kFALSE;
418  Double_t JetPtThreshold;
419  JetContHybridS->ResetCurrentID();
420  if(fJetShapeSub==kConstSub){
421  while((JetHybridS = JetContHybridS->GetNextAcceptJet())){
422  if (fJetShapeSub==kConstSub) JetPtThreshold=JetHybridS->Pt();
423  else JetPtThreshold=JetHybridS->Pt()-(GetRhoVal(0)*JetHybridS->Area());
424  if ( (!JetHybridS) || (JetPtThreshold<fPtThreshold)) continue;
425  Int_t JetNumber=-1;
426  for(Int_t i = 0; i<JetContHybridUS->GetNJets(); i++) {
427  JetHybridUS = JetContHybridUS->GetJet(i);
428  if(JetHybridUS->GetLabel()==JetHybridS->GetLabel()) {
429  JetNumber=i;
430  }
431  }
432  if(JetNumber==-1) continue;
433  JetHybridUS=JetContHybridUS->GetJet(JetNumber);
434  if (JetContHybridUS->AliJetContainer::GetFractionSharedPt(JetHybridUS)<fSharedFractionPtMin) continue;
435  JetPythDet=JetHybridUS->ClosestJet();
436 
437  if(!(fJetShapeSub==kConstSub)) JetHybridUS = JetHybridS->ClosestJet();
438  if (!JetHybridUS) {
439  Printf("Unsubtracted embedded jet does not exist, returning");
440  continue;
441  }
442  if (!JetPythDet) continue;
443  JetPythTrue=JetPythDet->ClosestJet();
444  if(!JetPythTrue) continue;
445  JetsMatched=kTRUE;
446 
447  fJetInfoVar[0]=JetHybridS->Pt();
448  fJetInfoVar[1]=JetPythTrue->Pt();
449  if(fDoSoftDrop) {
450  SoftDrop(JetHybridS,JetContHybridS,fZCut,fBeta_SD,kFALSE);
451  SoftDrop(JetPythTrue,JetContPythTrue,fZCut,fBeta_SD,kTRUE);
452  }
453  else{
454  fJetInfoVar[2]=0;
455  fJetInfoVar[3]=0;
456  fJetInfoVar[4]=0;
457  fJetInfoVar[5]=0;
458  fJetInfoVar[6]=0;
459  fJetInfoVar[7]=0;
460  fJetInfoVar[12]=0;
461  fJetInfoVar[13]=0;
462  fJetInfoVar[14]=0;
463  fJetInfoVar[15]=0;
464  fJetInfoVar[16]=0;
465  fJetInfoVar[17]=0;
466  }
467  fJetInfoVar[8]=PTD(JetHybridS,0);
468  fJetInfoVar[9]=PTD(JetPythTrue,0);
469  fJetInfoVar[10]=Angularity(JetHybridS,0);
470  fJetInfoVar[11]=Angularity(JetPythTrue,0);
471  fTreeJetInfo->Fill();
472  }
473  }
474  }
475  return kTRUE;
476 }
477 
478 //________________________________________________________________________
480 
481  if(Phi < -1*TMath::Pi()) Phi += (2*TMath::Pi());
482  else if (Phi > TMath::Pi()) Phi -= (2*TMath::Pi());
483  Double_t DeltaPhi=Phi-EventPlane;
484  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
485  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
486  return DeltaPhi;
487 }
488 //________________________________________________________________________
490 
491  if(Phi1 < -1*TMath::Pi()) Phi1 += (2*TMath::Pi());
492  else if (Phi1 > TMath::Pi()) Phi1 -= (2*TMath::Pi());
493  if(Phi2 < -1*TMath::Pi()) Phi2 += (2*TMath::Pi());
494  else if (Phi2 > TMath::Pi()) Phi2 -= (2*TMath::Pi());
495  Double_t DeltaPhi=Phi2-Phi1;
496  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
497  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
498  return DeltaPhi;
499 }
500 
501 
502 
503 
504 //--------------------------------------------------------------------------
506 
507  AliJetContainer *JetCont = GetJetContainer(JetContNb);
508  Double_t Angularity_Numerator=0;
509  Double_t Angularity_Denominator=0;
510  AliVParticle *Particle=0x0;
511  Double_t DeltaPhi=0;
512 
513 
514  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks
515  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
516  if(!Particle) continue;
517  DeltaPhi=RelativePhi(Jet->Phi(),Particle->Phi());
518  Angularity_Numerator=Angularity_Numerator+(Particle->Pt()*TMath::Sqrt(((Particle->Eta()-Jet->Eta())*(Particle->Eta()-Jet->Eta()))+(DeltaPhi*DeltaPhi)));
519  Angularity_Denominator= Angularity_Denominator+Particle->Pt();
520  }
521  if(Angularity_Denominator!=0) return Angularity_Numerator/Angularity_Denominator;
522  else return -1;
523 
524 }
525 
526 
527 
528 //--------------------------------------------------------------------------
530 
531  AliJetContainer *JetCont = GetJetContainer(JetContNb);
532  Double_t PTD_Numerator=0; //Reset these values
533  Double_t PTD_Denominator=0;
534  AliVParticle *Particle=0x0;
535  Double_t DeltaPhi=0;
536  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks
537  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
538  if(!Particle) continue;
539  PTD_Numerator=PTD_Numerator+(Particle->Pt()*Particle->Pt());
540  PTD_Denominator=PTD_Denominator+Particle->Pt();
541  }
542  if(PTD_Denominator!=0) return TMath::Sqrt(PTD_Numerator)/PTD_Denominator;
543  else return -1;
544 
545 }
546 
547 //_________________________________________________________________________
548  void AliAnalysisTaskRecoilJetYield::SoftDrop(AliEmcalJet *fJet,AliJetContainer *fJetCont, double zcut, double beta, Bool_t fTruthJet){
549  std::vector<fastjet::PseudoJet> fInputVectors;
550  Double_t JetInvMass=0, PseudJetInvMass=0, TrackMom = 0, TrackEnergy = 0;
551  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
552  Double_t JetEta=fJet->Eta(),JetPhi=fJet->Phi();
553  Double_t FJTrackEta[9999],FJTrackPhi[9999],FJTrackPt[9999],EmcalJetTrackEta[9999],EmcalJetTrackPhi[9999],EmcalJetTrackPt[9999];
554  UShort_t FJNTracks=0,EmcalJetNTracks=0;
555  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
556  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
557  JetInvMass += fTrk->M();
558  if (!fTrk) continue;
559  fastjet::PseudoJet PseudoTracks(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
560  TrackMom += TMath::Sqrt(TMath::Power(fTrk->Px(),2)+TMath::Power(fTrk->Py(),2)+TMath::Power(fTrk->Pz(),2));
561  TrackEnergy += fTrk->E();
562  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
563  PseudJetInvMass += PseudoTracks.m();
564  fInputVectors.push_back(PseudoTracks);
565  EmcalJetTrackEta[i]=fTrk->Eta();
566  EmcalJetTrackPhi[i]=fTrk->Phi();
567  EmcalJetTrackPt[i]=fTrk->Pt();
568  EmcalJetNTracks++;
569  }
570  fastjet::JetDefinition *fJetDef;
571  fastjet::ClusterSequence *fClustSeqSA;
572 
573 
574 
575  fJetDef = new fastjet::JetDefinition(fastjet::antikt_algorithm, fJetRadius*2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
576 
577  try {
578  fClustSeqSA = new fastjet::ClusterSequence(fInputVectors, *fJetDef);
579  } catch (fastjet::Error) {
580  AliError(" [w] FJ Exception caught.");
581  //return -1;
582  }
583 
584  std::vector<fastjet::PseudoJet> fOutputJets;
585  fOutputJets.clear();
586  fOutputJets=fClustSeqSA->inclusive_jets(0);
587 
588 
589  std::vector<fastjet::PseudoJet> jet_constituents = fOutputJets[0].constituents();
590  Float_t NSubjettinessResult1, NSubjettinessResult2, NSubBeta = 1, R0=0.4;
591  std::vector<fastjet::PseudoJet> Subjet_Axes;
592  fastjet::PseudoJet SubJet1_Axis,SubJet2_Axis;
593  if(jet_constituents.size()>=2){
594  fastjet::contrib::Nsubjettiness nSub1(1,fastjet::contrib::KT_Axes(),fastjet::contrib::NormalizedMeasure(NSubBeta,R0));
595  fastjet::contrib::Nsubjettiness nSub2(2,fastjet::contrib::KT_Axes(),fastjet::contrib::NormalizedMeasure(NSubBeta,R0));
596 
597  NSubjettinessResult1 = nSub1.result(fOutputJets[0]);
598  NSubjettinessResult2 = nSub2.result(fOutputJets[0]);
599  Subjet_Axes = nSub2.currentAxes();
600  SubJet1_Axis = Subjet_Axes[0];
601  SubJet2_Axis = Subjet_Axes[1];
602 
603  Double_t SubJet1_Eta=SubJet1_Axis.pseudorapidity();
604  Double_t SubJet2_Eta=SubJet2_Axis.pseudorapidity();
605  Double_t SubJet1_Phi=SubJet1_Axis.phi();
606  if(SubJet1_Phi < -1*TMath::Pi()) SubJet1_Phi += (2*TMath::Pi());
607  else if (SubJet1_Phi > TMath::Pi()) SubJet1_Phi -= (2*TMath::Pi());
608  Double_t SubJet2_Phi=SubJet2_Axis.phi();
609  if(SubJet1_Phi < -1*TMath::Pi()) SubJet1_Phi += (2*TMath::Pi());
610  else if (SubJet1_Phi > TMath::Pi()) SubJet1_Phi -= (2*TMath::Pi());
611 
612  Double_t DeltaPhiSubJets,DeltaEtaSubJets;
613  DeltaPhiSubJets=SubJet1_Phi-SubJet2_Phi;
614  DeltaEtaSubJets=SubJet1_Eta-SubJet2_Eta;
615  if(DeltaPhiSubJets < -1*TMath::Pi()) DeltaPhiSubJets += (2*TMath::Pi());
616  else if (DeltaPhiSubJets > TMath::Pi()) DeltaPhiSubJets -= (2*TMath::Pi());
617  Double_t DelR=0;
618  DelR = TMath::Power(TMath::Power(DeltaPhiSubJets,2)+TMath::Power(DeltaEtaSubJets,2),0.5);
619 
620  if(!fTruthJet) fJetInfoVar[12]=DelR;
621  else fJetInfoVar[13]=DelR;
622  }
623 
624  fastjet::contrib::SoftDrop softdrop(beta, zcut);
625  softdrop.set_verbose_structure(kTRUE);
626  fastjet::PseudoJet finaljet = softdrop(fOutputJets[0]);
627  if(!fTruthJet){
628  fJetInfoVar[4]=NSubjettinessResult1;
629  fJetInfoVar[6]=NSubjettinessResult2;
630  }
631  else {
632  fJetInfoVar[5]=NSubjettinessResult1;
633  fJetInfoVar[7]=NSubjettinessResult2;
634  }
635  AliEmcalJet* jet = new AliEmcalJet(finaljet.perp(), finaljet.eta(), finaljet.phi(), finaljet.m());
636  std::vector<fastjet::PseudoJet> fSDTracks=finaljet.constituents();
637  Double_t FastjetTrackDelR,EmcalTrackDelR;
638  for(Int_t i=0;i<fJet->GetNumberOfConstituents();i++){
639  if(i<=finaljet.constituents().size()){
640  FastjetTrackDelR = TMath::Sqrt(TMath::Power(fSDTracks[i].eta()-JetEta,2)+TMath::Power(fSDTracks[i].phi()-JetPhi,2));
641  FJTrackEta[i]=fSDTracks[i].eta();
642  FJTrackPhi[i]=fSDTracks[i].phi();
643  FJTrackPt[i]=fSDTracks[i].perp();
644  FJNTracks++;
645  }
646  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
647  EmcalTrackDelR = TMath::Sqrt(TMath::Power(fTrk->Eta()-JetEta,2)+TMath::Power(fTrk->Phi()-JetPhi,2));
648  }
649  Int_t NDroppedTracks = fJet->GetNumberOfTracks()-finaljet.constituents().size();
650  Int_t nConstituents(fClustSeqSA->constituents(finaljet).size());
651  jet->SetNumberOfTracks(nConstituents);
652  Double_t SymParam, Mu, DeltaR;
653  SymParam=(finaljet.structure_of<fastjet::contrib::SoftDrop>().symmetry());
654  Mu=(finaljet.structure_of<fastjet::contrib::SoftDrop>().mu());
655  DeltaR=(finaljet.structure_of<fastjet::contrib::SoftDrop>().delta_R());
656  fhGroomedPtvJetPt->Fill(finaljet.perp(),fJet->Pt());
657  fhDroppedBranches->Fill(finaljet.structure_of<fastjet::contrib::SoftDrop>().dropped_count());
658  if(!fTruthJet) fJetInfoVar[2]=SymParam;
659  else fJetInfoVar[3]=SymParam;
660  if(!fTruthJet) fJetInfoVar[14]=finaljet.structure_of<fastjet::contrib::SoftDrop>().dropped_count();
661  else fJetInfoVar[15]=finaljet.structure_of<fastjet::contrib::SoftDrop>().dropped_count();
662  if(!fTruthJet) fJetInfoVar[16]=finaljet.perp();
663  else fJetInfoVar[17]=finaljet.perp();
664 
665  return;
666 
667 
668 }
669 
670 //________________________________________________________________________
672  AliTrackContainer *PartCont = NULL;
673  if (fJetShapeSub==kConstSub) PartCont = GetTrackContainer(1);
674  else PartCont = GetTrackContainer(0);
675  TClonesArray *TracksArray = PartCont->GetArray();
676  if(!PartCont || !TracksArray) return -99999;
677  AliAODTrack *Track = 0x0;
678  Int_t Trigger_Index[10000];
679  for (Int_t i=0; i<100; i++) Trigger_Index[i] = 0;
680  Int_t Trigger_Counter = 0;
681  for(Int_t i=0; i < TracksArray->GetEntriesFast(); i++){
682  if((Track = static_cast<AliAODTrack*>(PartCont->GetAcceptTrack(i)))){
683  if (!Track) continue;
684  fhTrackPt->Fill(Track->Pt());
685  if(TMath::Abs(Track->Eta())>0.9) continue;
686  if (Track->Pt()<0.15) continue;
687  if ((Track->Pt() >= PtMin) && (Track->Pt()< PtMax)) {
688  Trigger_Index[Trigger_Counter] = i;
689  Trigger_Counter++;
690  }
691  }
692  }
693  if (Trigger_Counter == 0) return -99999;
694  Int_t RandomNumber = 0, Index = 0 ;
695  TRandom3* Random = new TRandom3(0);
696  RandomNumber = Random->Integer(Trigger_Counter);
697  Index = Trigger_Index[RandomNumber];
698  return Index;
699 }
700 
701 //________________________________________________________________________
703  //
704  // retrieve event objects
705  //
707  return kFALSE;
708 
709  return kTRUE;
710 }
711 
712 
713 //_______________________________________________________________________
715 {
716  // Called once at the end of the analysis.
717 
718 
719 }
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
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:219
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
Int_t GetLabel() const
Definition: AliEmcalJet.h:117
Double_t Eta_Lower
Float_t Percentile_1
Double_t Eta_Upper
void SoftDrop(AliEmcalJet *fJet, AliJetContainer *fJetCont, double zcut, double beta, Bool_t fTruthJet)
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)
Int_t GetNJets() const
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
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.
AliEmcalJet * GetJet(Int_t i) const