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