AliPhysics  c7b8e89 (c7b8e89)
AliAnalysisTaskEmcalQGTagging.cxx
Go to the documentation of this file.
1 //
2 // Jet QG tagging analysis task.
3 //
4 // Author: D. Caffarri, L. Cunqueiro
5 //just testing
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 
42 
43 
44 #include "AliAODEvent.h"
46 
47 using std::cout;
48 using std::endl;
49 
51 
52 //________________________________________________________________________
55  fContainer(0),
56  fMinFractionShared(0),
57  fJetShapeType(kData),
58  fJetShapeSub(kNoSub),
59  fJetSelection(kInclusive),
60  fPtThreshold(-9999.),
61  fRMatching(0.2),
62  fSelectedShapes(0),
63  fminpTTrig(20.),
64  fmaxpTTrig(50.),
65  fangWindowRecoil(0.6),
66  fSemigoodCorrect(0),
67  fHolePos(0),
68  fHoleWidth(0),
69  fCentSelectOn(kTRUE),
70  fCentMin(0),
71  fCentMax(10),
72  fOneConstSelectOn(kFALSE),
73  fTrackCheckPlots(kFALSE),
74  fSubjetCutoff(0.1),
75  fMinPtConst(1),
76  fHardCutoff(0),
77  fDoTwoTrack(kFALSE),
78  fPhiCutValue(0.02),
79  fEtaCutValue(0.02),
80  fMagFieldPolarity(1),
81  fDerivSubtrOrder(0),
82  fh2ResponseUW(0x0),
83  fh2ResponseW(0x0),
84  fPhiJetCorr6(0x0),
85  fPhiJetCorr7(0x0),
86  fEtaJetCorr6(0x0),
87  fEtaJetCorr7(0x0),
88  fPtJetCorr(0x0),
89  fPtJet(0x0),
90  fhpTjetpT(0x0),
91  fhPt(0x0),
92  fhPhi(0x0),
93  fhTrackPhi(0x0),
94  fHLundIterative(0x0),
95  fHCheckResolutionSubjets(0x0),
96  fNbOfConstvspT(0x0),
97  fTreeObservableTagging(0)
98 
99 {
100  for(Int_t i=0;i<17;i++){
101  fShapesVar[i]=0;}
102  SetMakeGeneralHistograms(kTRUE);
103  DefineOutput(1, TList::Class());
104  DefineOutput(2, TTree::Class());
105 }
106 
107 //________________________________________________________________________
109  AliAnalysisTaskEmcalJet(name, kTRUE),
110  fContainer(0),
111  fMinFractionShared(0),
112  fJetShapeType(kData),
113  fJetShapeSub(kNoSub),
114  fJetSelection(kInclusive),
115  fPtThreshold(-9999.),
116  fRMatching(0.2),
117  fSelectedShapes(0),
118  fminpTTrig(20.),
119  fmaxpTTrig(50.),
120  fangWindowRecoil(0.6),
121  fSemigoodCorrect(0),
122  fHolePos(0),
123  fHoleWidth(0),
124  fCentSelectOn(kTRUE),
125  fCentMin(0),
126  fCentMax(10),
127  fOneConstSelectOn(kFALSE),
128  fTrackCheckPlots(kFALSE),
129  fSubjetCutoff(0.1),
130  fMinPtConst(1),
131  fHardCutoff(0),
132  fDoTwoTrack(kFALSE),
133  fPhiCutValue(0.02),
134  fEtaCutValue(0.02),
135  fMagFieldPolarity(1),
136  fDerivSubtrOrder(0),
137  fh2ResponseUW(0x0),
138  fh2ResponseW(0x0),
139  fPhiJetCorr6(0x0),
140  fPhiJetCorr7(0x0),
141  fEtaJetCorr6(0x0),
142  fEtaJetCorr7(0x0),
143  fPtJetCorr(0x0),
144  fPtJet(0x0),
145  fhpTjetpT(0x0),
146  fhPt(0x0),
147  fhPhi(0x0),
148  fhTrackPhi(0x0),
149  fHLundIterative(0x0),
150  fHCheckResolutionSubjets(0x0),
151  fNbOfConstvspT(0x0),
152  fTreeObservableTagging(0)
153 
154 {
155  // Standard constructor.
156  for(Int_t i=0;i<17;i++){
157  fShapesVar[i]=0;}
159 
160  DefineOutput(1, TList::Class());
161  DefineOutput(2, TTree::Class());
162 }
163 
164 //________________________________________________________________________
166 {
167  // Destructor.
168 }
169 
170 //________________________________________________________________________
172 {
173  // Create user output.
174 
176 
177  Bool_t oldStatus = TH1::AddDirectoryStatus();
178  TH1::AddDirectory(kFALSE);
179 
180 
181 
182 
183  fh2ResponseUW= new TH2F("fh2ResponseUW", "fh2ResponseUW", 100, 0, 200, 100, 0, 200);
184  fOutput->Add(fh2ResponseUW);
185  fh2ResponseW= new TH2F("fh2ResponseW", "fh2ResponseW", 100, 0, 200, 100, 0, 200);
186  fOutput->Add(fh2ResponseW);
187  fPhiJetCorr6= new TH2F("fPhiJetCorr6", "fPhiJetCorr6", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
188  fOutput->Add(fPhiJetCorr6);
189  fEtaJetCorr6= new TH2F("fEtaJetCorr6", "fEtaJetCorr6", 50, -1.5, 1.5, 50, -1.5, 1.5);
190  fOutput->Add(fEtaJetCorr6);
191 
192  fPhiJetCorr7= new TH2F("fPhiJetCorr7", "fPhiJetCorr7", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
193  fOutput->Add(fPhiJetCorr7);
194  fEtaJetCorr7= new TH2F("fEtaJetCorr7", "fEtaJetCorr7", 50, -1.5, 1.5, 50, -1.5, 1.5);
195  fOutput->Add(fEtaJetCorr7);
196 
197  fPtJetCorr= new TH2F("fPtJetCorr", "fPtJetCorr", 100, 0, 200, 100, 0, 200);
198  fOutput->Add(fPtJetCorr);
199  fPtJet= new TH1F("fPtJet", "fPtJet", 100, 0, 200);
200  fOutput->Add(fPtJet);
201 
202  fhpTjetpT= new TH2F("fhpTjetpT", "fhpTjetpT", 200, 0, 200, 200, 0, 200);
203  fOutput->Add(fhpTjetpT);
204  fhPt= new TH1F("fhPt", "fhPt", 200, 0, 200);
205  fOutput->Add(fhPt);
206  fhPhi= new TH1F("fhPhi", "fhPhi", 100, -TMath::Pi(), TMath::Pi());
207  fOutput->Add(fhPhi);
208  fhTrackPhi= new TH1F("fhTrackPhi", "fhTrackPhi", 100, 0, 2*TMath::Pi());
209  fOutput->Add(fhTrackPhi);
210 
211 
212 
213  //log(1/theta),log(kt),jetpT,depth, algo, Eradiator//
214  const Int_t dimSpec = 7;
215  const Int_t nBinsSpec[7] = {50,100,100,20,3,100,2};
216  const Double_t lowBinSpec[7] = {0.,-10,0,0,0,0,0};
217  const Double_t hiBinSpec[7] = {5.,10.,100,20,3,100,2};
218  fHLundIterative = new THnSparseF("fHLundIterative",
219  "LundIterativePlot [log(1/theta),log(z*theta),pTjet,algo]",
220  dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
221  fOutput->Add(fHLundIterative);
222 
223 
225  const Int_t dimResol = 5;
226  const Int_t nBinsResol[5] = {10,10,80,80,80};
227  const Double_t lowBinResol[5] = {0,0,-1,-1,-1};
228  const Double_t hiBinResol[5] = {200,0.3,1,1,1};
229  fHCheckResolutionSubjets = new THnSparseF("fHCheckResolutionSubjets",
230  "Mom.Resolution of Subjets vs opening angle",
231  dimResol,nBinsResol,lowBinResol,hiBinResol);
233 
234 
235 
236 
237  fNbOfConstvspT=new TH2F("fNbOfConstvspT", "fNbOfConstvspT", 100, 0, 100, 200, 0, 200);
238  fOutput->Add(fNbOfConstvspT);
239 
240  // =========== Switch on Sumw2 for all histos ===========
241  for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
242  TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
243  if (h1){
244  h1->Sumw2();
245  continue;
246  }
247  THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
248  if(hn)hn->Sumw2();
249  }
250 
251 
252  TH1::AddDirectory(oldStatus);
253  const Int_t nVar = 14;
254  const char* nameoutput = GetOutputSlot(2)->GetContainer()->GetName();
255  fTreeObservableTagging = new TTree(nameoutput, nameoutput);
256 
257 
258  TString *fShapesVarNames = new TString [nVar];
259 
260  fShapesVarNames[0] = "partonCode";
261  fShapesVarNames[1] = "ptJet";
262  fShapesVarNames[2] = "ptDJet";
263  fShapesVarNames[3] = "phiJet";
264  // fShapesVarNames[4] = "nbOfConst";
265  fShapesVarNames[4] = "angularity";
266  //fShapesVarNames[5] = "circularity";
267  fShapesVarNames[5] = "lesub";
268  //fShapesVarNames[7] = "coronna";
269 
270  fShapesVarNames[6] = "ptJetMatch";
271  fShapesVarNames[7] = "ptDJetMatch";
272  fShapesVarNames[8] = "phiJetMatch";
273  // fShapesVarNames[12] = "nbOfConstMatch";
274  fShapesVarNames[9] = "angularityMatch";
275  //fShapesVarNames[12] = "circularityMatch";
276  fShapesVarNames[10] = "lesubMatch";
277  //fShapesVarNames[14] = "coronnaMatch";
278  fShapesVarNames[11]="weightPythia";
279  fShapesVarNames[12]="nsd";
280  fShapesVarNames[13]="nall";
281  //fShapesVarNames[14]="ntrksEvt";
282  //fShapesVarNames[16]="rhoVal";
283  //fShapesVarNames[17]="rhoMassVal";
284  //fShapesVarNames[12]="ptUnsub";
285 
286  for(Int_t ivar=0; ivar < nVar; ivar++){
287  cout<<"looping over variables"<<endl;
288  fTreeObservableTagging->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/F", fShapesVarNames[ivar].Data()));}
289 
290  PostData(1,fOutput);
291  PostData(2,fTreeObservableTagging);
292 
293  delete [] fShapesVarNames;
294 }
295 
296 //________________________________________________________________________
298 {
299  // Run analysis code here, if needed. It will be executed before FillHistograms().
300 
301  return kTRUE;
302 }
303 
304 //________________________________________________________________________
306 {
307  // Fill histograms.
308  //cout<<"base container"<<endl;
309  AliEmcalJet* jet1 = NULL;
310  AliJetContainer *jetCont = GetJetContainer(0);
311 
312  Float_t kWeight=1;
313  if (fCentSelectOn)
314  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
315 
316  AliAODTrack *triggerHadron = 0x0;
317 
318 
319  AliTrackContainer *PartCont =NULL;
320  AliParticleContainer *PartContMC=NULL;
321 
322  if (fJetShapeSub==kConstSub){
324  else PartCont = GetTrackContainer(1);
325  }
326  else{
328  else PartCont = GetTrackContainer(0);
329  }
330 
331  TClonesArray *TrackArray = NULL;
332  TClonesArray *TrackArrayMC = NULL;
333  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) TrackArrayMC = PartContMC->GetArray();
334  else TrackArray = PartCont->GetArray();
335 
336  Int_t NTracks=0;
337  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) NTracks = TrackArrayMC->GetEntriesFast();
338  else NTracks = TrackArray->GetEntriesFast();
339 
340 
341 
342 
343  if (fJetSelection == kRecoil) {
344  //Printf("Recoil jets!!!, fminpTTrig = %f, fmaxpTTrig = %f", fminpTTrig, fmaxpTTrig);
345  Int_t triggerHadronLabel = SelectTrigger(fminpTTrig, fmaxpTTrig);
346 
347 
348  if (triggerHadronLabel==-99999) {
349  //Printf ("Trigger Hadron not found, return");
350  return 0;}
351 
352  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) triggerHadron = static_cast<AliAODTrack*>(TrackArrayMC->At(triggerHadronLabel));
353  else triggerHadron = static_cast<AliAODTrack*>(TrackArray->At(triggerHadronLabel));
354 
356  if (!triggerHadron) {
357  //Printf("No Trigger hadron with the found label!!");
358  return 0;
359  }
360 
361  if(fSemigoodCorrect){
362  Double_t disthole=RelativePhi(triggerHadron->Phi(),fHolePos);
363  if(TMath::Abs(disthole)+fHoleWidth>TMath::Pi()-fangWindowRecoil){
364  return 0;}
365  }
366 
367 
368  fhPt->Fill(triggerHadron->Pt());
369 
370  }
371 
372 
373  if(fTrackCheckPlots){
374  //here check tracks//
375  AliAODTrack *Track = 0x0;
376  for(Int_t i=0; i < NTracks; i++){
378  if((Track = static_cast<AliAODTrack*>(PartContMC->GetAcceptParticle(i)))){
379  if (!Track) continue;
380  if(TMath::Abs(Track->Eta())>0.9) continue;
381  if (Track->Pt()<0.15) continue;
382  fhTrackPhi->Fill(Track->Phi());
383  }
384  }
385  else{
386  if((Track = static_cast<AliAODTrack*>(PartCont->GetAcceptTrack(i)))){
387  if (!Track) continue;
388  if(TMath::Abs(Track->Eta())>0.9) continue;
389  if (Track->Pt()<0.15) continue;
390  fhTrackPhi->Fill(Track->Phi());
391  }
392  }
393  }
394  }
395 
396 
397  // AliParticleContainer *partContAn = GetParticleContainer(0);
398  //TClonesArray *trackArrayAn = partContAn->GetArray();
399 
400 
401  Float_t rhoVal=0, rhoMassVal = 0.;
402  if(jetCont) {
403 
404  jetCont->ResetCurrentID();
406  //rho
407  AliRhoParameter* rhoParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject("RhoSparseR020"));
408  if (!rhoParam) {
409  Printf("%s: Could not retrieve rho %s (some histograms will be filled with zero)!", GetName(), jetCont->GetRhoName().Data());
410  } else rhoVal = rhoParam->GetVal();
411  //rhom
412  AliRhoParameter* rhomParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject("RhoMassSparseR020"));
413  if (!rhomParam) {
414  Printf("%s: Could not retrieve rho_m %s (some histograms will be filled with zero)!", GetName(), jetCont->GetRhoMassName().Data());
415  } else rhoMassVal = rhomParam->GetVal();
416  }
417 
418  while((jet1 = jetCont->GetNextAcceptJet())) {
419  if (!jet1) continue;
420  AliEmcalJet* jet2 = 0x0;
421  AliEmcalJet* jet3 = 0x0;
422  fPtJet->Fill(jet1->Pt());
423  AliEmcalJet *jetUS = NULL;
424  Int_t ifound=0, jfound=0;
425  Int_t ilab=-1, jlab=-1;
426 
428  Double_t disthole=RelativePhi(jet1->Phi(),fHolePos);
429  if(TMath::Abs(disthole)<fHoleWidth){
430  continue;}
431  }
432 
433  Float_t dphiRecoil = 0.;
434  if (fJetSelection == kRecoil){
435  dphiRecoil = RelativePhi(triggerHadron->Phi(), jet1->Phi());
436  if (TMath::Abs(dphiRecoil) < (TMath::Pi() - fangWindowRecoil)) {
437  // Printf("Recoil jets back to back not found! continuing");
438  continue;
439  }
440 
441  fhpTjetpT->Fill(triggerHadron->Pt(), jet1->Pt());
442  //Printf(" ************ FILLING HISTOS****** shapeSub = %d, triggerHadron = %f, jet1 = %f", fJetShapeSub, triggerHadron->Pt(), jet1->Pt());
443  fhPhi->Fill(RelativePhi(triggerHadron->Phi(), jet1->Phi()));
444 
445  }
446 
447 
448  fShapesVar[0] = 0.;
450  AliJetContainer *jetContTrue = GetJetContainer(1);
451  AliJetContainer *jetContUS = GetJetContainer(2);
452 
453  if(fJetShapeSub==kConstSub){
454  for(Int_t i = 0; i<jetContUS->GetNJets(); i++) {
455  jetUS = jetContUS->GetJet(i);
456  if(jetUS->GetLabel()==jet1->GetLabel()) {
457  ifound++;
458  if(ifound==1) ilab = i;
459  }
460  }
461  if(ilab==-1) continue;
462  jetUS=jetContUS->GetJet(ilab);
463  jet2=jetUS->ClosestJet();
464  }
465 
466  if(!(fJetShapeSub==kConstSub)) jet2 = jet1->ClosestJet();
467  if (!jet2) {
468  Printf("jet2 does not exist, returning");
469  continue;
470  }
471 
472  AliJetContainer *jetContPart=GetJetContainer(3);
473  jet3=jet2->ClosestJet();
474 
475  if(!jet3){
476  Printf("jet3 does not exist, returning");
477  continue;
478  }
479  cout<<"jet 3 exists"<<jet3->Pt()<<endl;
480 
481 
482  fh2ResponseUW->Fill(jet1->Pt(),jet2->Pt());
483 
484  Double_t fraction=0;
485  if(!(fJetShapeSub==kConstSub)) fraction = jetCont->GetFractionSharedPt(jet1);
486  if(fJetShapeSub==kConstSub) fraction = jetContUS->GetFractionSharedPt(jetUS);
487  //if (fraction > 0.1) cout<<"***** hey a jet matched with fraction"<<fraction<<" "<<jet1->Pt()<<" "<<jet2->Pt()<<" "<<fCent<<endl;
488 
489  if(fraction<fMinFractionShared) continue;
490  //InputEvent()->Print();
491 
492  }
493 
494 
495 
496  if (fJetShapeType == kPythiaDef){
497 
498  AliJetContainer *jetContTrue = GetJetContainer(1);
499  AliJetContainer *jetContUS = GetJetContainer(2);
500  AliJetContainer *jetContPart = GetJetContainer(3);
501 
502  if(fJetShapeSub==kConstSub){
503 
504  for(Int_t i = 0; i<jetContUS->GetNJets(); i++) {
505  jetUS = jetContUS->GetJet(i);
506  if(jetUS->GetLabel()==jet1->GetLabel()) {
507  ifound++;
508  if(ifound==1) ilab = i;
509  }
510  }
511  if(ilab==-1) continue;
512  jetUS=jetContUS->GetJet(ilab);
513  jet2=jetUS->ClosestJet();
514 
515  if (!jet2) {
516  Printf("jet2 does not exist, returning");
517  continue;
518  }
519 
520  for(Int_t j=0; j<jetContPart->GetNJets(); j++) {
521 
522  jet3 = jetContPart->GetJet(j);
523  if(!jet3) continue;
524  if(jet3->GetLabel()==jet2->GetLabel()) {
525  jfound++;
526  if(jfound==1) jlab = j;
527  }
528  }
529  if(jlab==-1) continue;
530  jet3=jetContPart->GetJet(jlab);
531  if(!jet3){
532  Printf("jet3 does not exist, returning");
533  continue;
534  }
535  }
536  if(!(fJetShapeSub==kConstSub)) jet3 = jet1->ClosestJet();
537  if (!jet3) {
538  Printf("jet3 does not exist, returning");
539  continue;
540  }
541 
542 
543  fh2ResponseUW->Fill(jet1->Pt(),jet3->Pt());
544  CheckSubjetResolution(jet1,jetCont,jet3,jetContTrue);
545 
546  }
547 
548 
549  if (fJetShapeType == kGenOnTheFly){
550  const AliEmcalPythiaInfo *partonsInfo = 0x0;
551  partonsInfo = GetPythiaInfo();
552  Double_t jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi6());
553  Double_t detap1=(jet1->Eta())-(partonsInfo->GetPartonEta6());
554  kWeight=partonsInfo->GetPythiaEventWeight();
555  fh2ResponseW->Fill(jet1->Pt(),jet1->Pt(),kWeight);
556 
557  Float_t dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
558  fEtaJetCorr6->Fill(jet1->Eta(), partonsInfo->GetPartonEta6());
559  fPhiJetCorr6->Fill(jet1->Phi(), partonsInfo->GetPartonPhi6());
560  if(dRp1 < fRMatching) {
561  fShapesVar[0] = partonsInfo->GetPartonFlag6();
562  fPtJetCorr ->Fill(partonsInfo->GetPartonPt6(), jet1->Pt());
563  }
564  else {
565  jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi7());
566  detap1=(jet1->Eta())-(partonsInfo->GetPartonEta7());
567  dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
568  fEtaJetCorr7->Fill(jet1->Eta(), partonsInfo->GetPartonEta7());
569  fPhiJetCorr7->Fill(jet1->Phi(), partonsInfo->GetPartonPhi7());
570  if(dRp1 < fRMatching) {
571  fShapesVar[0] = partonsInfo->GetPartonFlag7();
572  fPtJetCorr->Fill(partonsInfo->GetPartonPt7(), jet1->Pt());
573  }
574  else fShapesVar[0]=0;
575  }
576  }
577 
578 
579 
580 
581  Double_t ptSubtracted = 0;
582  if (fJetShapeSub==kConstSub) ptSubtracted= jet1->Pt();
583 
584  else if (fJetShapeSub==kDerivSub) {
585  ptSubtracted=jet1->Pt()-GetRhoVal(0)*jet1->Area();
586  }
587 
588  else if (fJetShapeSub==kNoSub) {
589  if ((fJetShapeType==kData) || (fJetShapeType==kDetEmbPartPythia)) ptSubtracted=jet1->Pt()-GetRhoVal(0)*jet1->Area();
590  else if ((fJetShapeType==kPythiaDef) || (fJetShapeType==kMCTrue) || (fJetShapeType==kGenOnTheFly)) ptSubtracted= jet1->Pt();
591  }
592 
593  //Printf("ptSubtracted=%f,fPtThreshold =%f ", ptSubtracted, fPtThreshold);
594  if (ptSubtracted < fPtThreshold) continue;
595 
596  if (fOneConstSelectOn == kTRUE) fNbOfConstvspT->Fill(GetJetNumberOfConstituents(jet1,0), ptSubtracted);
597  if ((fCentSelectOn == kFALSE) && (jet1->GetNumberOfTracks() <= 1)) continue;
598 
599 
600  fShapesVar[1] = ptSubtracted;
601  fShapesVar[2] = GetJetpTD(jet1,0);
602  fShapesVar[3] =jet1->Phi();
603  if(fJetShapeType==kData && fJetSelection == kRecoil) fShapesVar[3]=RelativePhi(triggerHadron->Phi(), jet1->Phi());
604  //GetJetMass(jet1,0);
605  fShapesVar[4] = GetJetAngularity(jet1,0);
606  //fShapesVar[5] = GetJetCircularity(jet1,0);
607  fShapesVar[5] = GetJetLeSub(jet1,0);
608  //fShapesVar[6] = GetJetCoronna(jet1,0);
609  RecursiveParents(jet1,jetCont,0);
610  RecursiveParents(jet1,jetCont,1);
611  RecursiveParents(jet1,jetCont,2);
612 
613  Float_t ptMatch=0., ptDMatch=0., massMatch=0.,angulMatch=0.,lesubMatch=0;
614  Int_t kMatched = 0;
615 
616  if (fJetShapeType==kPythiaDef) {
617  kMatched =1;
618  if(fJetShapeSub==kConstSub) kMatched = 3;
619 
620  ptMatch=jet3->Pt();
621  ptDMatch=GetJetpTD(jet3, kMatched);
622  massMatch=jet3->Phi();
623  // GetJetMass(jet3,kMatched);
624  //constMatch=1.*GetJetNumberOfConstituents(jet2,kMatched);
625  angulMatch=GetJetAngularity(jet3, kMatched);
626  //circMatch=GetJetCircularity(jet3, kMatched);
627  lesubMatch=GetJetLeSub(jet3, kMatched);
628  //coronnaMatch=GetJetCoronna(jet3,kMatched);
629  //sigma2Match = GetSigma2(jet2, kMatched);
630  }
631 
633  if(fJetShapeSub==kConstSub) kMatched = 3;
634  if(fJetShapeSub==kDerivSub) kMatched = 2;
635  ptMatch=jet3->Pt();
636  ptDMatch=GetJetpTD(jet3, kMatched);
637  massMatch=jet2->Pt();
638  //GetJetMass(jet3,kMatched);
639  // constMatch=1.*GetJetNumberOfConstituents(jet3,kMatched);
640  angulMatch=GetJetAngularity(jet3, kMatched);
641  // circMatch=GetJetCircularity(jet3, kMatched);
642  lesubMatch=GetJetLeSub(jet3, kMatched);
643  //coronnaMatch = GetJetCoronna(jet3, kMatched);
644 
645  }
646 
647 
648 
650  kMatched = 0;
651  ptMatch=0.;
652  ptDMatch=0.;
653  massMatch=0.;
654  //constMatch=0.;
655  angulMatch=0.;
656  // circMatch=0.;
657  lesubMatch=0.;
658  //coronnaMatch =0.;
659 
660  }
661 
662 
663 
664  fShapesVar[6] = ptMatch;
665  fShapesVar[7] = ptDMatch;
666  fShapesVar[8] = massMatch;
667  fShapesVar[9] = angulMatch;
668  //fShapesVar[12] = circMatch;
669  fShapesVar[10] = lesubMatch;
670  // fShapesVar[14] = coronnaMatch;
671  fShapesVar[11] = kWeight;
672  //fShapesVar[16] = ntracksEvt;
673  // fShapesVar[16] = rhoVal;
674  //fShapesVar[17] = rhoMassVal;
675  //fShapesVar[16] = jet1->Pt();
676 
677 
678  fTreeObservableTagging->Fill();
679 
680 
681 
682 
683 
684 
685  }
686 
687  }
688 
689  return kTRUE;
690 }
691 
692 //________________________________________________________________________
694  //calc subtracted jet mass
695  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
697  else return jet->GetShapeProperties()->GetSecondOrderSubtracted();
698  else
699  return jet->M();
700 }
701 
702 //________________________________________________________________________
704 
705  AliJetContainer *jetCont = GetJetContainer(jetContNb);
706  if (!jet->GetNumberOfTracks())
707  return 0;
708  Double_t den=0.;
709  Double_t num = 0.;
710  AliVParticle *vp1 = 0x0;
711  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
712  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
713 
714  if (!vp1){
715  Printf("AliVParticle associated to constituent not found");
716  continue;
717  }
718 
719  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
720  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
721  Double_t dr = TMath::Sqrt(dr2);
722  num=num+vp1->Pt()*dr;
723  den=den+vp1->Pt();
724  }
725  return num/den;
726 }
727 
728 //________________________________________________________________________
730 
731  if((fJetShapeSub==kDerivSub) && (jetContNb==0))
734  else
735  return Angularity(jet, jetContNb);
736 
737 }
738 
739 //____________________________________________________________________________
740 
742 
743  AliTrackContainer *PartCont = NULL;
744  AliParticleContainer *PartContMC = NULL;
745 
746 
747  if (fJetShapeSub==kConstSub){
749  else PartCont = GetTrackContainer(1);
750  }
751  else{
753  else PartCont = GetTrackContainer(0);
754  }
755 
756  TClonesArray *TracksArray = NULL;
757  TClonesArray *TracksArrayMC = NULL;
758 
759  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) TracksArrayMC = PartContMC->GetArray();
760  else TracksArray = PartCont->GetArray();
761 
763  if(!PartContMC || !TracksArrayMC) return -2;
764  }
765  else {
766  if(!PartCont || !TracksArray) return -2;
767  }
768 
769 
770  AliAODTrack *Track = 0x0;
771  Float_t sumpt=0;
772  Int_t NTracks=0;
773  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) NTracks = TracksArrayMC->GetEntriesFast();
774  else NTracks = TracksArray->GetEntriesFast();
775 
776  for(Int_t i=0; i < NTracks; i++){
778  if((Track = static_cast<AliAODTrack*>(PartContMC->GetAcceptParticle(i)))){
779  if (!Track) continue;
780  if(TMath::Abs(Track->Eta())>0.9) continue;
781  Double_t dphi = RelativePhi(Track->Phi(),jet->Phi());
782  Double_t dr2 = (Track->Eta()-jet->Eta())*(Track->Eta()-jet->Eta()) + dphi*dphi;
783  Double_t dr = TMath::Sqrt(dr2);
784  if((dr>=0.8) && (dr<1)) sumpt=sumpt+Track->Pt();
785  }
786  }
787  else{
788  if((Track = static_cast<AliAODTrack*>(PartCont->GetAcceptTrack(i)))){
789  if (!Track) continue;
790  if(TMath::Abs(Track->Eta())>0.9) continue;
791  if (Track->Pt()<0.15) continue;
792  Double_t dphi = RelativePhi(Track->Phi(),jet->Phi());
793  Double_t dr2 = (Track->Eta()-jet->Eta())*(Track->Eta()-jet->Eta()) + dphi*dphi;
794  Double_t dr = TMath::Sqrt(dr2);
795  if((dr>=0.8) && (dr<1)) sumpt=sumpt+Track->Pt();
796 
797  }
798  }
799  }
800 
801 
802 
803  return sumpt;
804 
805 
806 
807 
808 }
809 
810 //________________________________________________________________________
812 
813  if((fJetShapeSub==kDerivSub) && (jetContNb==0)) return -2;
814  else
815  return Coronna(jet, jetContNb);
816 
817 }
818 
819 
820 
821 
822 
823 //________________________________________________________________________
825 
826  AliJetContainer *jetCont = GetJetContainer(jetContNb);
827  if (!jet->GetNumberOfTracks())
828  return 0;
829  Double_t den=0.;
830  Double_t num = 0.;
831  AliVParticle *vp1 = 0x0;
832  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
833  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
834 
835  if (!vp1){
836  Printf("AliVParticle associated to constituent not found");
837  continue;
838  }
839 
840  num=num+vp1->Pt()*vp1->Pt();
841  den=den+vp1->Pt();
842  }
843  return TMath::Sqrt(num)/den;
844 }
845 
846 //________________________________________________________________________
848  //calc subtracted jet mass
849  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
851  else return jet->GetShapeProperties()->GetSecondOrderSubtractedpTD();
852  else
853  return PTD(jet, jetContNb);
854 
855 }
856 
857 //_____________________________________________________________________________
859 
860  AliJetContainer *jetCont = GetJetContainer(jetContNb);
861  if (!jet->GetNumberOfTracks())
862  return 0;
863  Double_t mxx = 0.;
864  Double_t myy = 0.;
865  Double_t mxy = 0.;
866  int nc = 0;
867  Double_t sump2 = 0.;
868  Double_t pxjet=jet->Px();
869  Double_t pyjet=jet->Py();
870  Double_t pzjet=jet->Pz();
871 
872 
873  //2 general normalized vectors perpendicular to the jet
874  TVector3 ppJ1(pxjet, pyjet, pzjet);
875  TVector3 ppJ3(- pxjet* pzjet, - pyjet * pzjet, pxjet * pxjet + pyjet * pyjet);
876  ppJ3.SetMag(1.);
877  TVector3 ppJ2(-pyjet, pxjet, 0);
878  ppJ2.SetMag(1.);
879  AliVParticle *vp1 = 0x0;
880  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
881  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
882 
883  if (!vp1){
884  Printf("AliVParticle associated to constituent not found");
885  continue;
886  }
887 
888  TVector3 pp(vp1->Px(), vp1->Py(), vp1->Pz());
889 
890  //local frame
891  TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
892  TVector3 pPerp = pp - pLong;
893  //projection onto the two perpendicular vectors defined above
894 
895  Float_t ppjX = pPerp.Dot(ppJ2);
896  Float_t ppjY = pPerp.Dot(ppJ3);
897  Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
898  if(ppjT<=0) return 0;
899 
900  mxx += (ppjX * ppjX / ppjT);
901  myy += (ppjY * ppjY / ppjT);
902  mxy += (ppjX * ppjY / ppjT);
903  nc++;
904  sump2 += ppjT;}
905 
906  if(nc<2) return 0;
907  if(sump2==0) return 0;
908  // Sphericity Matrix
909  Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
910  TMatrixDSym m0(2,ele);
911 
912  // Find eigenvectors
913  TMatrixDSymEigen m(m0);
914  TVectorD eval(2);
915  TMatrixD evecm = m.GetEigenVectors();
916  eval = m.GetEigenValues();
917  // Largest eigenvector
918  int jev = 0;
919  // cout<<eval[0]<<" "<<eval[1]<<endl;
920  if (eval[0] < eval[1]) jev = 1;
921  TVectorD evec0(2);
922  // Principle axis
923  evec0 = TMatrixDColumn(evecm, jev);
924  Double_t compx=evec0[0];
925  Double_t compy=evec0[1];
926  TVector2 evec(compx, compy);
927  Double_t circ=0;
928  if(jev==1) circ=2*eval[0];
929  if(jev==0) circ=2*eval[1];
930 
931  return circ;
932 
933 
934 
935 }
936 
937 
938 
939 
940 //________________________________________________________________________
942  //calc subtracted jet mass
943 
944  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
947  else
948  return Circularity(jet, jetContNb);
949 
950 }
951 
952 //________________________________________________________________________
954 
955  AliJetContainer *jetCont = GetJetContainer(jetContNb);
956  if (!jet->GetNumberOfTracks())
957  return 0;
958  Double_t den=0.;
959  Double_t num = 0.;
960  AliVParticle *vp1 = 0x0;
961  AliVParticle *vp2 = 0x0;
962  std::vector<int> ordindex;
963  ordindex=jet->GetPtSortedTrackConstituentIndexes(jetCont->GetParticleContainer()->GetArray());
964  //Printf("Nbof const = %d", jet->GetNumberOfTracks());
965  //Printf("ordindex[0] = %d, ordindex[1] = %d", ordindex[0], ordindex[1]);
966 
967  if(ordindex.size()<2) return -1;
968 
969  vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[0], jetCont->GetParticleContainer()->GetArray()));
970  if (!vp1){
971  Printf("AliVParticle associated to Leading constituent not found");
972  return -1;
973  }
974 
975  vp2 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[1], jetCont->GetParticleContainer()->GetArray()));
976  if (!vp2){
977  Printf("AliVParticle associated to Subleading constituent not found");
978  return -1;
979  }
980 
981 
982  num=vp1->Pt();
983  den=vp2->Pt();
984  //Printf("vp1->Pt() =%f, vp2->Pt() =%f", vp1->Pt(), vp2->Pt());
985 
986 return num-den;
987 }
988 
989 //________________________________________________________________________
991  //calc subtracted jet mass
992 
993  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
996  else
997  return LeSub(jet, jetContNb);
998 
999 }
1000 
1001 //________________________________________________________________________
1003  //calc subtracted jet mass
1004 
1005  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
1008  else
1009  return jet->GetNumberOfTracks();
1010 
1011  }
1012 
1013 
1014 //______________________________________________________________________________
1016 
1017  AliJetContainer *jetCont = GetJetContainer(jetContNb);
1018  if (!jet->GetNumberOfTracks())
1019  return 0;
1020  Double_t mxx = 0.;
1021  Double_t myy = 0.;
1022  Double_t mxy = 0.;
1023  int nc = 0;
1024  Double_t sump2 = 0.;
1025 
1026  AliVParticle *vp1 = 0x0;
1027  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
1028  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
1029 
1030  if (!vp1){
1031  Printf("AliVParticle associated to constituent not found");
1032  continue;
1033  }
1034 
1035  Double_t ppt=vp1->Pt();
1036  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
1037 
1038  Double_t deta = vp1->Eta()-jet->Eta();
1039  mxx += ppt*ppt*deta*deta;
1040  myy += ppt*ppt*dphi*dphi;
1041  mxy -= ppt*ppt*deta*TMath::Abs(dphi);
1042  nc++;
1043  sump2 += ppt*ppt;
1044 
1045  }
1046  if(nc<2) return 0;
1047  if(sump2==0) return 0;
1048  // Sphericity Matrix
1049  Double_t ele[4] = {mxx , mxy , mxy , myy };
1050  TMatrixDSym m0(2,ele);
1051 
1052  // Find eigenvectors
1053  TMatrixDSymEigen m(m0);
1054  TVectorD eval(2);
1055  TMatrixD evecm = m.GetEigenVectors();
1056  eval = m.GetEigenValues();
1057  // Largest eigenvector
1058  int jev = 0;
1059  // cout<<eval[0]<<" "<<eval[1]<<endl;
1060  if (eval[0] < eval[1]) jev = 1;
1061  TVectorD evec0(2);
1062  // Principle axis
1063  evec0 = TMatrixDColumn(evecm, jev);
1064  Double_t compx=evec0[0];
1065  Double_t compy=evec0[1];
1066  TVector2 evec(compx, compy);
1067  Double_t sig=0;
1068  if(jev==1) sig=TMath::Sqrt(TMath::Abs(eval[0])/sump2);
1069  if(jev==0) sig=TMath::Sqrt(TMath::Abs(eval[1])/sump2);
1070 
1071  return sig;
1072 
1073 }
1074 
1075 //________________________________________________________________________
1077  //calc subtracted jet mass
1078 
1079  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
1081  else return jet->GetShapeProperties()->GetSecondOrderSubtractedSigma2();
1082  else
1083  return Sigma2(jet, jetContNb);
1084 
1085 }
1086 
1087 //________________________________________________________________________
1089 
1090  AliTrackContainer *PartCont = NULL;
1091  AliParticleContainer *PartContMC = NULL;
1092 
1093 
1094  if (fJetShapeSub==kConstSub){
1096  else PartCont = GetTrackContainer(1);
1097  }
1098  else{
1100  else PartCont = GetTrackContainer(0);
1101  }
1102 
1103  TClonesArray *TracksArray = NULL;
1104  TClonesArray *TracksArrayMC = NULL;
1105 
1106  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) TracksArrayMC = PartContMC->GetArray();
1107  else TracksArray = PartCont->GetArray();
1108 
1110  if(!PartContMC || !TracksArrayMC) return -99999;
1111  }
1112  else {
1113  if(!PartCont || !TracksArray) return -99999;
1114  }
1115 
1116 
1117  AliAODTrack *Track = 0x0;
1118 
1119 
1120 
1121  TList trackList;
1122  Int_t triggers[100];
1123  for (Int_t iTrigger=0; iTrigger<100; iTrigger++) triggers[iTrigger] = 0;
1124  Int_t iTT = 0;
1125  Int_t NTracks=0;
1126  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) NTracks = TracksArrayMC->GetEntriesFast();
1127  else NTracks = TracksArray->GetEntriesFast();
1128 
1129  for(Int_t i=0; i < NTracks; i++){
1131  if((Track = static_cast<AliAODTrack*>(PartContMC->GetAcceptParticle(i)))){
1132  if (!Track) continue;
1133  if(TMath::Abs(Track->Eta())>0.9) continue;
1134  if (Track->Pt()<0.15) continue;
1135  if ((Track->Pt() >= minpT) && (Track->Pt()< maxpT)) {
1136  triggers[iTT] = i;
1137  iTT++;
1138  }
1139  }
1140  }
1141  else{
1142  if((Track = static_cast<AliAODTrack*>(PartCont->GetAcceptTrack(i)))){
1143  if (!Track) continue;
1144  if(TMath::Abs(Track->Eta())>0.9) continue;
1145  if (Track->Pt()<0.15) continue;
1146  if ((Track->Pt() >= minpT) && (Track->Pt()< maxpT)) {
1147  triggers[iTT] = i;
1148  iTT++;
1149  }
1150  }
1151  }
1152  }
1153 
1154 
1155  if (iTT == 0) return -99999;
1156  Int_t nbRn = 0, index = 0 ;
1157  TRandom3 random(0);
1158  nbRn = random.Integer(iTT);
1159  index = triggers[nbRn];
1160  //Printf("iTT Total= %d, nbRn = %d, Index = %d",iTT, nbRn, index );
1161  return index;
1162 
1163 }
1164 
1165 //__________________________________________________________________________________
1167 
1168  if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1169  else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1170  if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1171  else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1172  double dphi = mphi-vphi;
1173  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1174  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1175  return dphi;//dphi in [-Pi, Pi]
1176 }
1177 
1178 
1179 //_________________________________________________________________________
1181 
1182  std::vector<fastjet::PseudoJet> fInputVectors;
1183  fInputVectors.clear();
1184  fastjet::PseudoJet PseudoTracks;
1185  double xflagalgo=0;
1186  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
1187 
1188  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
1189  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
1190  if (!fTrk) continue;
1191  if(fDoTwoTrack==kTRUE && CheckClosePartner(i,fJet,fTrk,fTrackCont)) continue;
1192  PseudoTracks.reset(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
1193  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
1194  fInputVectors.push_back(PseudoTracks);
1195 
1196  }
1197  fastjet::JetAlgorithm jetalgo(fastjet::antikt_algorithm);
1198 
1199  if(ReclusterAlgo==0){ xflagalgo=0.5;
1200  jetalgo=fastjet::kt_algorithm ;}
1201 
1202  if(ReclusterAlgo==1){ xflagalgo=1.5;
1203  jetalgo=fastjet::cambridge_algorithm;}
1204  if(ReclusterAlgo==2){ xflagalgo=2.5;
1205  jetalgo=fastjet::antikt_algorithm;}
1206 
1207  fastjet::JetDefinition fJetDef(jetalgo, 1., static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
1208 
1209  try {
1210  fastjet::ClusterSequence fClustSeqSA(fInputVectors, fJetDef);
1211  std::vector<fastjet::PseudoJet> fOutputJets;
1212  fOutputJets.clear();
1213  fOutputJets=fClustSeqSA.inclusive_jets(0);
1214 
1215  fastjet::PseudoJet jj;
1216  fastjet::PseudoJet j1;
1217  fastjet::PseudoJet j2;
1218  jj=fOutputJets[0];
1219  double ndepth=0;
1220  double nall=0;
1221  double flagSubjet=0;
1222  while(jj.has_parents(j1,j2)){
1223  nall=nall+1;
1224  if(j1.perp() < j2.perp()) swap(j1,j2);
1225  flagSubjet=0;
1226  double delta_R=j1.delta_R(j2);
1227  double z=j2.perp()/(j1.perp()+j2.perp());
1228  double y =log(1.0/delta_R);
1229  double lnpt_rel=log(j2.perp()*delta_R);
1230  double yh=j1.e()+j2.e();
1231  vector < fastjet::PseudoJet > constitj1 = sorted_by_pt(j1.constituents());
1232  if(constitj1[0].perp()>fMinPtConst) flagSubjet=1;
1233  if(z>fHardCutoff){
1234  ndepth=ndepth+1;
1235  Double_t LundEntries[7] = {y,lnpt_rel,fOutputJets[0].perp(),nall,xflagalgo,yh,flagSubjet};
1236  fHLundIterative->Fill(LundEntries);}
1237  jj=j1;}
1238 
1239  if(ReclusterAlgo==1){
1240  fShapesVar[12]=nall;
1241  fShapesVar[13]=ndepth;
1242  }
1243 
1244  } catch (fastjet::Error) {
1245  AliError(" [w] FJ Exception caught.");
1246  //return -1;
1247  }
1248 
1249 
1250 
1251 
1252  return;
1253 
1254 
1255 }
1256 
1257 //_________________________________________________________________________
1259 
1260  std::vector<fastjet::PseudoJet> fInputVectors;
1261  fInputVectors.clear();
1262  fastjet::PseudoJet PseudoTracks;
1263 
1264  std::vector<fastjet::PseudoJet> fInputVectorsM;
1265  fInputVectorsM.clear();
1266  fastjet::PseudoJet PseudoTracksM;
1267 
1268  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
1269  AliParticleContainer *fTrackContM = fJetContM->GetParticleContainer();
1270 
1271  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
1272  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
1273  if (!fTrk) continue;
1274  PseudoTracks.reset(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
1275  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
1276  fInputVectors.push_back(PseudoTracks);
1277 
1278  }
1279  fastjet::JetAlgorithm jetalgo(fastjet::cambridge_algorithm);
1280  fastjet::JetDefinition fJetDef(jetalgo, 1., static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
1281 
1282  if (fTrackContM) for (Int_t i=0; i<fJetM->GetNumberOfTracks(); i++) {
1283  AliVParticle *fTrk = fJetM->TrackAt(i, fTrackContM->GetArray());
1284  if (!fTrk) continue;
1285  PseudoTracksM.reset(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
1286  PseudoTracksM.set_user_index(fJetM->TrackAt(i)+100);
1287  fInputVectorsM.push_back(PseudoTracksM);
1288 
1289  }
1290  fastjet::JetAlgorithm jetalgoM(fastjet::cambridge_algorithm);
1291  fastjet::JetDefinition fJetDefM(jetalgoM, 1., static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
1292 
1293 
1294  try {
1295  fastjet::ClusterSequence fClustSeqSA(fInputVectors, fJetDef);
1296  std::vector<fastjet::PseudoJet> fOutputJets;
1297  fOutputJets.clear();
1298  fOutputJets=fClustSeqSA.inclusive_jets(0);
1299 
1300  fastjet::ClusterSequence fClustSeqSAM(fInputVectorsM, fJetDefM);
1301  std::vector<fastjet::PseudoJet> fOutputJetsM;
1302  fOutputJetsM.clear();
1303  fOutputJetsM=fClustSeqSAM.inclusive_jets(0);
1304 
1305 
1306 
1307  fastjet::PseudoJet jj,jjM;
1308  fastjet::PseudoJet j1,j1M;
1309  fastjet::PseudoJet j2,j2M;
1310  jj=fOutputJets[0];
1311  jjM=fOutputJetsM[0];
1312 
1313  double z1=0;
1314  double z2=0;
1315  double zcut=0.1;
1316  while((jj.has_parents(j1,j2)) && (z1<zcut)){
1317  if(j1.perp() < j2.perp()) swap(j1,j2);
1318 
1319  z1=j2.perp()/(j1.perp()+j2.perp());
1320  jj=j1;}
1321  if(z1<zcut) return;
1322 
1323 
1324 
1325  while((jjM.has_parents(j1M,j2M)) && (z2<zcut)){
1326  if(j1M.perp() < j2M.perp()) swap(j1M,j2M);
1327 
1328  z2=j2M.perp()/(j1M.perp()+j2M.perp());
1329  jjM=j1M;}
1330  if(z2<zcut) return;
1331 
1332 
1333 
1334  double delta_R1=j1.delta_R(j1M);
1335  double delta_R2=j2.delta_R(j2M);
1336  double delta_R=j1.delta_R(j2);
1337  double residz=(z1-z2)/z2;
1338  double resid1=(j1.perp()-j1M.perp())/j1M.perp();
1339  double resid2=(j2.perp()-j2M.perp())/j2M.perp();
1340 
1341  if((delta_R1<fSubjetCutoff) && (delta_R2<fSubjetCutoff)){
1342  Double_t ResolEntries[5] = {fOutputJets[0].perp(),delta_R,resid1,resid2,residz};
1343  fHCheckResolutionSubjets->Fill(ResolEntries);}
1344 
1345 
1346  } catch (fastjet::Error) {
1347  AliError(" [w] FJ Exception caught.");
1348  //return -1;
1349  }
1350 
1351 
1352 
1353 
1354  return;
1355 
1356 
1357 
1358 
1359 
1360 
1361 
1362 }
1363 
1364 
1366  //check if tracks are close//
1367  for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
1368  AliVParticle *fTrk2 = fJet->TrackAt(i, fTrackCont->GetArray());
1369  if (!fTrk2) continue;
1370  if(i==index) continue;
1371  Double_t phi1 = fTrk1->Phi();
1372  Double_t phi2 = fTrk2->Phi();
1373  Double_t chg1 = fTrk1->Charge();
1374  Double_t chg2 = fTrk2->Charge();
1375  Double_t ptv1 = fTrk1->Pt();
1376  Double_t ptv2 = fTrk2->Pt();
1377  Double_t deta=fTrk2->Eta()-fTrk1->Eta();
1378  const Float_t kLimit = fPhiCutValue* 3;
1379 
1380  if (TMath::Abs(fTrk1->Eta()-fTrk2->Eta()) < fEtaCutValue*2.5*3){
1381  Float_t initdpsinner = (phi2 - TMath::ASin(0.075*chg2*fMagFieldPolarity*0.8/ptv2) -
1382  (phi1 - TMath::ASin(0.075*chg1*fMagFieldPolarity*0.8/ptv1)));
1383 
1384  Float_t initdpsouter = (phi2 - TMath::ASin(0.075*chg2*fMagFieldPolarity*2.5/ptv2) -
1385  (phi1 - TMath::ASin(0.075*chg1*fMagFieldPolarity*2.5/ptv1)));
1386 
1387  initdpsinner = TVector2::Phi_mpi_pi(initdpsinner);
1388  initdpsouter = TVector2::Phi_mpi_pi(initdpsouter);
1389 
1390 
1391 
1392  if (TMath::Abs(initdpsinner) < kLimit ||
1393  TMath::Abs(initdpsouter) < kLimit || initdpsinner * initdpsouter < 0 ) {
1394  Double_t mindps = 1e5;
1395 
1396 
1397 
1398  for (Double_t rad = 0.8; rad < 2.51; rad += 0.01) {
1399  Double_t dps = (phi2 - TMath::ASin(0.075*chg2*fMagFieldPolarity*rad/ptv2) - (phi1 - TMath::ASin(0.075*chg1*fMagFieldPolarity*rad/ptv1)));
1400  dps = TVector2::Phi_mpi_pi(dps);
1401  if (TMath::Abs(dps) < TMath::Abs(mindps))
1402  mindps = dps;
1403  }
1404  if(TMath::Abs(mindps)<fPhiCutValue && TMath::Abs(deta)<fEtaCutValue) return kTRUE;
1405  } }
1406  return kFALSE;
1407  }}
1408 
1409 
1410 
1411 
1412 
1413 //________________________________________________________________________
1415  //
1416  // retrieve event objects
1417  //
1419  return kFALSE;
1420 
1421  return kTRUE;
1422 }
1423 
1424 //_______________________________________________________________________
1426 {
1427  // Called once at the end of the analysis.
1428 
1429  // fTreeObservableTagging = dynamic_cast<TTree*>(GetOutputData(1));
1430  // if (!fTreeObservableTagging){
1431  // Printf("ERROR: fTreeObservableTagging not available");
1432  // return;
1433  // }
1434 
1435 }
1436 
Int_t SelectTrigger(Float_t minpT, Float_t maxpT)
Float_t Coronna(AliEmcalJet *jet, Int_t jetContNb)
Double_t GetFirstOrderSubtractedAngularity() const
Double_t Area() const
Definition: AliEmcalJet.h:130
const TString & GetRhoName() const
double Double_t
Definition: External.C:58
Double_t GetSecondOrderSubtractedSigma2() const
Definition: External.C:236
Bool_t RetrieveEventObjects()
Retrieve common objects from event.
Float_t GetPartonEta6() const
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:327
AliJetContainer * GetJetContainer(Int_t i=0) const
Bool_t FillHistograms()
Function filling histograms.
Float_t GetPartonEta7() const
Double_t GetSecondOrderSubtractedConstituent() const
Double_t Eta() const
Definition: AliEmcalJet.h:121
Double_t Py() const
Definition: AliEmcalJet.h:107
Double_t Phi() const
Definition: AliEmcalJet.h:117
Float_t GetJetMass(AliEmcalJet *jet, Int_t jetContNb)
Float_t GetPythiaEventWeight() const
Container with name, TClonesArray and cuts for particles.
Int_t GetLabel() const
Definition: AliEmcalJet.h:124
Float_t GetPartonPhi7() const
Float_t GetJetAngularity(AliEmcalJet *jet, Int_t jetContNb)
Int_t GetPartonFlag7() const
Container for particles within the EMCAL framework.
Float_t GetPartonPhi6() const
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
Float_t GetJetCoronna(AliEmcalJet *jet, Int_t jetContNb)
TString kData
Declare data MC or deltaAOD.
Double_t Px() const
Definition: AliEmcalJet.h:106
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
const TString & GetRhoMassName() const
Int_t GetPartonFlag6() const
AliParticleContainer * GetParticleContainer() const
Double_t GetFirstOrderSubtractedConstituent() const
int Int_t
Definition: External.C:63
Float_t GetJetCircularity(AliEmcalJet *jet, Int_t jetContNb)
unsigned int UInt_t
Definition: External.C:33
Float_t PTD(AliEmcalJet *jet, Int_t jetContNb)
float Float_t
Definition: External.C:68
Int_t GetNJets() const
Double_t GetSecondOrderSubtractedAngularity() const
Float_t GetSigma2(AliEmcalJet *jet, Int_t jetContNb)
virtual AliVParticle * GetAcceptParticle(Int_t i=-1) const
Double_t fCent
!event centrality
void RecursiveParents(AliEmcalJet *fJet, AliJetContainer *fJetCont, Int_t ReclusterAlgo)
AliEmcalJet * GetNextAcceptJet()
Float_t Angularity(AliEmcalJet *jet, Int_t jetContNb)
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 LeSub(AliEmcalJet *jet, Int_t jetContNb)
Double_t Pt() const
Definition: AliEmcalJet.h:109
Double_t GetRhoVal(Int_t i=0) const
Double_t GetFirstOrderSubtractedCircularity() const
void Track()
AliEmcalList * fOutput
!output list
Float_t Sigma2(AliEmcalJet *jet, Int_t jetContNb)
Double_t RelativePhi(Double_t mphi, Double_t vphi)
Double_t GetSecondOrderSubtractedCircularity() const
Float_t GetJetNumberOfConstituents(AliEmcalJet *jet, Int_t jetContNb)
AliTrackContainer * GetTrackContainer(Int_t i=0) const
Bool_t CheckClosePartner(Int_t index, AliEmcalJet *fJet, AliVParticle *fTrack, AliParticleContainer *fTrackCont)
Store some informaion about a Pythia eventThis class is used to store some information about a Pythia...
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:51
const AliEmcalPythiaInfo * GetPythiaInfo() const
Float_t GetJetpTD(AliEmcalJet *jet, Int_t jetContNb)
Double_t Pz() const
Definition: AliEmcalJet.h:108
Double_t GetFractionSharedPt(const AliEmcalJet *jet, AliParticleContainer *cont2=0x0) const
Declaration of class AliEmcalPythiaInfo.
const char Option_t
Definition: External.C:48
void UserCreateOutputObjects()
Main initialization function on the worker.
Double_t GetFirstOrderSubtractedSigma2() const
Float_t Circularity(AliEmcalJet *jet, Int_t jetContNb)
bool Bool_t
Definition: External.C:53
AliEmcalJetShapeProperties * GetShapeProperties() const
Definition: AliEmcalJet.h:361
Float_t GetJetLeSub(AliEmcalJet *jet, Int_t jetContNb)
Float_t GetPartonPt7() const
void swap(PWGJE::EMCALJetTasks::AliAnalysisTaskEmcalJetHPerformance &first, PWGJE::EMCALJetTasks::AliAnalysisTaskEmcalJetHPerformance &second)
Double_t M() const
Definition: AliEmcalJet.h:120
Container for jet within the EMCAL jet framework.
Bool_t Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.
Definition: External.C:196
std::vector< int > GetPtSortedTrackConstituentIndexes(TClonesArray *tracks) const
Float_t GetPartonPt6() const
void CheckSubjetResolution(AliEmcalJet *fJet, AliJetContainer *fJetCont, AliEmcalJet *fJetM, AliJetContainer *fJetContM)
AliEmcalJet * GetJet(Int_t i) const