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