AliPhysics  b0c77bb (b0c77bb)
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  fCheckResolution(kFALSE),
75  fSubjetCutoff(0.1),
76  fMinPtConst(1),
77  fHardCutoff(0),
78  fDoTwoTrack(kFALSE),
79  fPhiCutValue(0.02),
80  fEtaCutValue(0.02),
81  fMagFieldPolarity(1),
82  fDerivSubtrOrder(0),
83  fh2ResponseUW(0x0),
84  fh2ResponseW(0x0),
85  fPhiJetCorr6(0x0),
86  fPhiJetCorr7(0x0),
87  fEtaJetCorr6(0x0),
88  fEtaJetCorr7(0x0),
89  fPtJetCorr(0x0),
90  fPtJet(0x0),
91  fhpTjetpT(0x0),
92  fhPt(0x0),
93  fhPhi(0x0),
94  fhTrackPhi(0x0),
95  fHLundIterative(0x0),
96  fHCheckResolutionSubjets(0x0),
97  fNbOfConstvspT(0x0),
98  fTreeObservableTagging(0)
99 
100 {
101  for(Int_t i=0;i<8;i++){
102  fShapesVar[i]=0;}
103  SetMakeGeneralHistograms(kTRUE);
104  DefineOutput(1, TList::Class());
105  DefineOutput(2, TTree::Class());
106 }
107 
108 //________________________________________________________________________
110  AliAnalysisTaskEmcalJet(name, kTRUE),
111  fContainer(0),
112  fMinFractionShared(0),
113  fJetShapeType(kData),
114  fJetShapeSub(kNoSub),
115  fJetSelection(kInclusive),
116  fPtThreshold(-9999.),
117  fRMatching(0.2),
118  fSelectedShapes(0),
119  fminpTTrig(20.),
120  fmaxpTTrig(50.),
121  fangWindowRecoil(0.6),
122  fSemigoodCorrect(0),
123  fHolePos(0),
124  fHoleWidth(0),
125  fCentSelectOn(kTRUE),
126  fCentMin(0),
127  fCentMax(10),
128  fOneConstSelectOn(kFALSE),
129  fTrackCheckPlots(kFALSE),
130  fCheckResolution(kFALSE),
131  fSubjetCutoff(0.1),
132  fMinPtConst(1),
133  fHardCutoff(0),
134  fDoTwoTrack(kFALSE),
135  fPhiCutValue(0.02),
136  fEtaCutValue(0.02),
137  fMagFieldPolarity(1),
138  fDerivSubtrOrder(0),
139  fh2ResponseUW(0x0),
140  fh2ResponseW(0x0),
141  fPhiJetCorr6(0x0),
142  fPhiJetCorr7(0x0),
143  fEtaJetCorr6(0x0),
144  fEtaJetCorr7(0x0),
145  fPtJetCorr(0x0),
146  fPtJet(0x0),
147  fhpTjetpT(0x0),
148  fhPt(0x0),
149  fhPhi(0x0),
150  fhTrackPhi(0x0),
151  fHLundIterative(0x0),
152  fHCheckResolutionSubjets(0x0),
153  fNbOfConstvspT(0x0),
154  fTreeObservableTagging(0)
155 
156 {
157  // Standard constructor.
158  for(Int_t i=0;i<8;i++){
159  fShapesVar[i]=0;}
161 
162  DefineOutput(1, TList::Class());
163  DefineOutput(2, TTree::Class());
164 }
165 
166 //________________________________________________________________________
168 {
169  // Destructor.
170 }
171 
172 //________________________________________________________________________
174 {
175  // Create user output.
176 
178 
179  Bool_t oldStatus = TH1::AddDirectoryStatus();
180  TH1::AddDirectory(kFALSE);
181 
182 
183 
184 
185  fh2ResponseUW= new TH2F("fh2ResponseUW", "fh2ResponseUW", 100, 0, 200, 100, 0, 200);
186  fOutput->Add(fh2ResponseUW);
187  fh2ResponseW= new TH2F("fh2ResponseW", "fh2ResponseW", 100, 0, 200, 100, 0, 200);
188  fOutput->Add(fh2ResponseW);
189  fPhiJetCorr6= new TH2F("fPhiJetCorr6", "fPhiJetCorr6", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
190  fOutput->Add(fPhiJetCorr6);
191  fEtaJetCorr6= new TH2F("fEtaJetCorr6", "fEtaJetCorr6", 50, -1.5, 1.5, 50, -1.5, 1.5);
192  fOutput->Add(fEtaJetCorr6);
193 
194  fPhiJetCorr7= new TH2F("fPhiJetCorr7", "fPhiJetCorr7", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
195  fOutput->Add(fPhiJetCorr7);
196  fEtaJetCorr7= new TH2F("fEtaJetCorr7", "fEtaJetCorr7", 50, -1.5, 1.5, 50, -1.5, 1.5);
197  fOutput->Add(fEtaJetCorr7);
198 
199  fPtJetCorr= new TH2F("fPtJetCorr", "fPtJetCorr", 100, 0, 200, 100, 0, 200);
200  fOutput->Add(fPtJetCorr);
201  fPtJet= new TH1F("fPtJet", "fPtJet", 100, 0, 200);
202  fOutput->Add(fPtJet);
203 
204  fhpTjetpT= new TH2F("fhpTjetpT", "fhpTjetpT", 200, 0, 200, 200, 0, 200);
205  fOutput->Add(fhpTjetpT);
206  fhPt= new TH1F("fhPt", "fhPt", 200, 0, 200);
207  fOutput->Add(fhPt);
208  fhPhi= new TH1F("fhPhi", "fhPhi", 100, -TMath::Pi(), TMath::Pi());
209  fOutput->Add(fhPhi);
210  fhTrackPhi= new TH1F("fhTrackPhi", "fhTrackPhi", 100, 0, 2*TMath::Pi());
211  fOutput->Add(fhTrackPhi);
212 
213 
214 
215  //log(1/theta),log(kt),jetpT,depth, algo, Eradiator//
216  const Int_t dimSpec = 6;
217  const Int_t nBinsSpec[6] = {50,100,100,20,100,2};
218  const Double_t lowBinSpec[6] = {0.,-10,0,0,0,0};
219  const Double_t hiBinSpec[6] = {5.,10.,100,20,100,2};
220  fHLundIterative = new THnSparseF("fHLundIterative",
221  "LundIterativePlot [log(1/theta),log(z*theta),pTjet,algo]",
222  dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
223  fOutput->Add(fHLundIterative);
224 
225 
227  const Int_t dimResol = 5;
228  const Int_t nBinsResol[5] = {10,10,80,80,80};
229  const Double_t lowBinResol[5] = {0,0,-1,-1,-1};
230  const Double_t hiBinResol[5] = {200,0.3,1,1,1};
231  fHCheckResolutionSubjets = new THnSparseF("fHCheckResolutionSubjets",
232  "Mom.Resolution of Subjets vs opening angle",
233  dimResol,nBinsResol,lowBinResol,hiBinResol);
235 
236 
237 
238 
239  fNbOfConstvspT=new TH2F("fNbOfConstvspT", "fNbOfConstvspT", 100, 0, 100, 200, 0, 200);
240  fOutput->Add(fNbOfConstvspT);
241 
242  // =========== Switch on Sumw2 for all histos ===========
243  for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
244  TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
245  if (h1){
246  h1->Sumw2();
247  continue;
248  }
249  THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
250  if(hn)hn->Sumw2();
251  }
252 
253 
254  TH1::AddDirectory(oldStatus);
255  const Int_t nVar = 8;
256  const char* nameoutput = GetOutputSlot(2)->GetContainer()->GetName();
257  fTreeObservableTagging = new TTree(nameoutput, nameoutput);
258 
259 
260  TString *fShapesVarNames = new TString [nVar];
261 
262  fShapesVarNames[0] = "partonCode";
263  fShapesVarNames[1] = "ptJet";
264  fShapesVarNames[2] = "ktAv";
265  fShapesVarNames[3] = "thetaAv";
266  fShapesVarNames[4] = "ptJetMatch";
267  fShapesVarNames[5] = "ktAvMatch";
268  fShapesVarNames[6] = "thetaAvMatch";
269  fShapesVarNames[7]="weightPythia";
270 
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 
383  // AliParticleContainer *partContAn = GetParticleContainer(0);
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  if(fCheckResolution) CheckSubjetResolution(jet1,jetCont,jet3,jetContTrue);
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  RecursiveParents(jet1,jetCont);
588 
589 
590  Float_t ptMatch=0.;
591  Int_t kMatched = 0;
592  Double_t ktAvMatch=0;;
593  Double_t thetaAvMatch=0;
594  Double_t aver1=0;
595  Double_t aver2=0;
596  if (fJetShapeType==kPythiaDef) {
597  kMatched =1;
598  if(fJetShapeSub==kConstSub) kMatched = 3;
599 
600  ptMatch=jet3->Pt();
601  RecursiveParentsMCAverage(jet3,kMatched, aver1, aver2);
602  ktAvMatch=aver1;
603  thetaAvMatch=aver2;
604 
605  }
606 
608  if(fJetShapeSub==kConstSub) kMatched = 3;
609  if(fJetShapeSub==kDerivSub) kMatched = 2;
610  ptMatch=jet3->Pt();
611 
612  RecursiveParentsMCAverage(jet3,kMatched, aver1, aver2);
613  ktAvMatch=aver1;
614  thetaAvMatch=aver2;
615 
616  }
617 
618 
619 
621  kMatched = 0;
622  ptMatch=0.;
623  ktAvMatch=0.;
624  thetaAvMatch=0.;
625 
626 
627  }
628 
629 
630 
631  fShapesVar[4] = ptMatch;
632  fShapesVar[5] = ktAvMatch;
633  fShapesVar[6] = thetaAvMatch;
634  fShapesVar[7] = kWeight;
635  fTreeObservableTagging->Fill();
636 
637 
638 
639 
640 
641 
642  }
643 
644  }
645 
646  return kTRUE;
647 }
648 
649 //________________________________________________________________________
651  //calc subtracted jet mass
652  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
654  else return jet->GetShapeProperties()->GetSecondOrderSubtracted();
655  else
656  return jet->M();
657 }
658 
659 //________________________________________________________________________
661 
662  AliJetContainer *jetCont = GetJetContainer(jetContNb);
663  if (!jet->GetNumberOfTracks())
664  return 0;
665  Double_t den=0.;
666  Double_t num = 0.;
667  AliVParticle *vp1 = 0x0;
668  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
669  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
670 
671  if (!vp1){
672  Printf("AliVParticle associated to constituent not found");
673  continue;
674  }
675 
676  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
677  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
678  Double_t dr = TMath::Sqrt(dr2);
679  num=num+vp1->Pt()*dr;
680  den=den+vp1->Pt();
681  }
682  return num/den;
683 }
684 
685 //________________________________________________________________________
687 
688  if((fJetShapeSub==kDerivSub) && (jetContNb==0))
691  else
692  return Angularity(jet, jetContNb);
693 
694 }
695 
696 //____________________________________________________________________________
697 
699 
700  AliTrackContainer *PartCont = NULL;
701  AliParticleContainer *PartContMC = NULL;
702 
703 
704  if (fJetShapeSub==kConstSub){
706  else PartCont = GetTrackContainer(1);
707  }
708  else{
710  else PartCont = GetTrackContainer(0);
711  }
712 
713  TClonesArray *TracksArray = NULL;
714  TClonesArray *TracksArrayMC = NULL;
715 
716  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) TracksArrayMC = PartContMC->GetArray();
717  else TracksArray = PartCont->GetArray();
718 
720  if(!PartContMC || !TracksArrayMC) return -2;
721  }
722  else {
723  if(!PartCont || !TracksArray) return -2;
724  }
725 
726 
727  AliAODTrack *Track = 0x0;
728  Float_t sumpt=0;
729  Int_t NTracks=0;
730  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) NTracks = TracksArrayMC->GetEntriesFast();
731  else NTracks = TracksArray->GetEntriesFast();
732 
733  for(Int_t i=0; i < NTracks; i++){
735  if((Track = static_cast<AliAODTrack*>(PartContMC->GetAcceptParticle(i)))){
736  if (!Track) continue;
737  if(TMath::Abs(Track->Eta())>0.9) continue;
738  Double_t dphi = RelativePhi(Track->Phi(),jet->Phi());
739  Double_t dr2 = (Track->Eta()-jet->Eta())*(Track->Eta()-jet->Eta()) + dphi*dphi;
740  Double_t dr = TMath::Sqrt(dr2);
741  if((dr>=0.8) && (dr<1)) sumpt=sumpt+Track->Pt();
742  }
743  }
744  else{
745  if((Track = static_cast<AliAODTrack*>(PartCont->GetAcceptTrack(i)))){
746  if (!Track) continue;
747  if(TMath::Abs(Track->Eta())>0.9) continue;
748  if (Track->Pt()<0.15) continue;
749  Double_t dphi = RelativePhi(Track->Phi(),jet->Phi());
750  Double_t dr2 = (Track->Eta()-jet->Eta())*(Track->Eta()-jet->Eta()) + dphi*dphi;
751  Double_t dr = TMath::Sqrt(dr2);
752  if((dr>=0.8) && (dr<1)) sumpt=sumpt+Track->Pt();
753 
754  }
755  }
756  }
757 
758 
759 
760  return sumpt;
761 
762 
763 
764 
765 }
766 
767 //________________________________________________________________________
769 
770  if((fJetShapeSub==kDerivSub) && (jetContNb==0)) return -2;
771  else
772  return Coronna(jet, jetContNb);
773 
774 }
775 
776 
777 
778 
779 
780 //________________________________________________________________________
782 
783  AliJetContainer *jetCont = GetJetContainer(jetContNb);
784  if (!jet->GetNumberOfTracks())
785  return 0;
786  Double_t den=0.;
787  Double_t num = 0.;
788  AliVParticle *vp1 = 0x0;
789  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
790  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
791 
792  if (!vp1){
793  Printf("AliVParticle associated to constituent not found");
794  continue;
795  }
796 
797  num=num+vp1->Pt()*vp1->Pt();
798  den=den+vp1->Pt();
799  }
800  return TMath::Sqrt(num)/den;
801 }
802 
803 //________________________________________________________________________
805  //calc subtracted jet mass
806  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
808  else return jet->GetShapeProperties()->GetSecondOrderSubtractedpTD();
809  else
810  return PTD(jet, jetContNb);
811 
812 }
813 
814 //_____________________________________________________________________________
816 
817  AliJetContainer *jetCont = GetJetContainer(jetContNb);
818  if (!jet->GetNumberOfTracks())
819  return 0;
820  Double_t mxx = 0.;
821  Double_t myy = 0.;
822  Double_t mxy = 0.;
823  int nc = 0;
824  Double_t sump2 = 0.;
825  Double_t pxjet=jet->Px();
826  Double_t pyjet=jet->Py();
827  Double_t pzjet=jet->Pz();
828 
829 
830  //2 general normalized vectors perpendicular to the jet
831  TVector3 ppJ1(pxjet, pyjet, pzjet);
832  TVector3 ppJ3(- pxjet* pzjet, - pyjet * pzjet, pxjet * pxjet + pyjet * pyjet);
833  ppJ3.SetMag(1.);
834  TVector3 ppJ2(-pyjet, pxjet, 0);
835  ppJ2.SetMag(1.);
836  AliVParticle *vp1 = 0x0;
837  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
838  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
839 
840  if (!vp1){
841  Printf("AliVParticle associated to constituent not found");
842  continue;
843  }
844 
845  TVector3 pp(vp1->Px(), vp1->Py(), vp1->Pz());
846 
847  //local frame
848  TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
849  TVector3 pPerp = pp - pLong;
850  //projection onto the two perpendicular vectors defined above
851 
852  Float_t ppjX = pPerp.Dot(ppJ2);
853  Float_t ppjY = pPerp.Dot(ppJ3);
854  Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
855  if(ppjT<=0) return 0;
856 
857  mxx += (ppjX * ppjX / ppjT);
858  myy += (ppjY * ppjY / ppjT);
859  mxy += (ppjX * ppjY / ppjT);
860  nc++;
861  sump2 += ppjT;}
862 
863  if(nc<2) return 0;
864  if(sump2==0) return 0;
865  // Sphericity Matrix
866  Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
867  TMatrixDSym m0(2,ele);
868 
869  // Find eigenvectors
870  TMatrixDSymEigen m(m0);
871  TVectorD eval(2);
872  TMatrixD evecm = m.GetEigenVectors();
873  eval = m.GetEigenValues();
874  // Largest eigenvector
875  int jev = 0;
876  // cout<<eval[0]<<" "<<eval[1]<<endl;
877  if (eval[0] < eval[1]) jev = 1;
878  TVectorD evec0(2);
879  // Principle axis
880  evec0 = TMatrixDColumn(evecm, jev);
881  Double_t compx=evec0[0];
882  Double_t compy=evec0[1];
883  TVector2 evec(compx, compy);
884  Double_t circ=0;
885  if(jev==1) circ=2*eval[0];
886  if(jev==0) circ=2*eval[1];
887 
888  return circ;
889 
890 
891 
892 }
893 
894 
895 
896 
897 //________________________________________________________________________
899  //calc subtracted jet mass
900 
901  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
904  else
905  return Circularity(jet, jetContNb);
906 
907 }
908 
909 //________________________________________________________________________
911 
912  AliJetContainer *jetCont = GetJetContainer(jetContNb);
913  if (!jet->GetNumberOfTracks())
914  return 0;
915  Double_t den=0.;
916  Double_t num = 0.;
917  AliVParticle *vp1 = 0x0;
918  AliVParticle *vp2 = 0x0;
919  std::vector<int> ordindex;
920  ordindex=jet->GetPtSortedTrackConstituentIndexes(jetCont->GetParticleContainer()->GetArray());
921  //Printf("Nbof const = %d", jet->GetNumberOfTracks());
922  //Printf("ordindex[0] = %d, ordindex[1] = %d", ordindex[0], ordindex[1]);
923 
924  if(ordindex.size()<2) return -1;
925 
926  vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[0], jetCont->GetParticleContainer()->GetArray()));
927  if (!vp1){
928  Printf("AliVParticle associated to Leading constituent not found");
929  return -1;
930  }
931 
932  vp2 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[1], jetCont->GetParticleContainer()->GetArray()));
933  if (!vp2){
934  Printf("AliVParticle associated to Subleading constituent not found");
935  return -1;
936  }
937 
938 
939  num=vp1->Pt();
940  den=vp2->Pt();
941  //Printf("vp1->Pt() =%f, vp2->Pt() =%f", vp1->Pt(), vp2->Pt());
942 
943 return num-den;
944 }
945 
946 //________________________________________________________________________
948  //calc subtracted jet mass
949 
950  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
953  else
954  return LeSub(jet, jetContNb);
955 
956 }
957 
958 //________________________________________________________________________
960  //calc subtracted jet mass
961 
962  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
965  else
966  return jet->GetNumberOfTracks();
967 
968  }
969 
970 
971 //______________________________________________________________________________
973 
974  AliJetContainer *jetCont = GetJetContainer(jetContNb);
975  if (!jet->GetNumberOfTracks())
976  return 0;
977  Double_t mxx = 0.;
978  Double_t myy = 0.;
979  Double_t mxy = 0.;
980  int nc = 0;
981  Double_t sump2 = 0.;
982 
983  AliVParticle *vp1 = 0x0;
984  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
985  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
986 
987  if (!vp1){
988  Printf("AliVParticle associated to constituent not found");
989  continue;
990  }
991 
992  Double_t ppt=vp1->Pt();
993  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
994 
995  Double_t deta = vp1->Eta()-jet->Eta();
996  mxx += ppt*ppt*deta*deta;
997  myy += ppt*ppt*dphi*dphi;
998  mxy -= ppt*ppt*deta*TMath::Abs(dphi);
999  nc++;
1000  sump2 += ppt*ppt;
1001 
1002  }
1003  if(nc<2) return 0;
1004  if(sump2==0) return 0;
1005  // Sphericity Matrix
1006  Double_t ele[4] = {mxx , mxy , mxy , myy };
1007  TMatrixDSym m0(2,ele);
1008 
1009  // Find eigenvectors
1010  TMatrixDSymEigen m(m0);
1011  TVectorD eval(2);
1012  TMatrixD evecm = m.GetEigenVectors();
1013  eval = m.GetEigenValues();
1014  // Largest eigenvector
1015  int jev = 0;
1016  // cout<<eval[0]<<" "<<eval[1]<<endl;
1017  if (eval[0] < eval[1]) jev = 1;
1018  TVectorD evec0(2);
1019  // Principle axis
1020  evec0 = TMatrixDColumn(evecm, jev);
1021  Double_t compx=evec0[0];
1022  Double_t compy=evec0[1];
1023  TVector2 evec(compx, compy);
1024  Double_t sig=0;
1025  if(jev==1) sig=TMath::Sqrt(TMath::Abs(eval[0])/sump2);
1026  if(jev==0) sig=TMath::Sqrt(TMath::Abs(eval[1])/sump2);
1027 
1028  return sig;
1029 
1030 }
1031 
1032 //________________________________________________________________________
1034  //calc subtracted jet mass
1035 
1036  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
1038  else return jet->GetShapeProperties()->GetSecondOrderSubtractedSigma2();
1039  else
1040  return Sigma2(jet, jetContNb);
1041 
1042 }
1043 
1044 //________________________________________________________________________
1046 
1047  AliTrackContainer *PartCont = NULL;
1048  AliParticleContainer *PartContMC = NULL;
1049 
1050 
1051  if (fJetShapeSub==kConstSub){
1053  else PartCont = GetTrackContainer(1);
1054  }
1055  else{
1057  else PartCont = GetTrackContainer(0);
1058  }
1059 
1060  TClonesArray *TracksArray = NULL;
1061  TClonesArray *TracksArrayMC = NULL;
1062 
1063  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) TracksArrayMC = PartContMC->GetArray();
1064  else TracksArray = PartCont->GetArray();
1065 
1067  if(!PartContMC || !TracksArrayMC) return -99999;
1068  }
1069  else {
1070  if(!PartCont || !TracksArray) return -99999;
1071  }
1072 
1073 
1074  AliAODTrack *Track = 0x0;
1075 
1076 
1077 
1078  TList trackList;
1079  Int_t triggers[100];
1080  for (Int_t iTrigger=0; iTrigger<100; iTrigger++) triggers[iTrigger] = 0;
1081  Int_t iTT = 0;
1082  Int_t NTracks=0;
1083  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) NTracks = TracksArrayMC->GetEntriesFast();
1084  else NTracks = TracksArray->GetEntriesFast();
1085 
1086  for(Int_t i=0; i < NTracks; i++){
1088  if((Track = static_cast<AliAODTrack*>(PartContMC->GetAcceptParticle(i)))){
1089  if (!Track) continue;
1090  if(TMath::Abs(Track->Eta())>0.9) continue;
1091  if (Track->Pt()<0.15) continue;
1092  if ((Track->Pt() >= minpT) && (Track->Pt()< maxpT)) {
1093  triggers[iTT] = i;
1094  iTT++;
1095  }
1096  }
1097  }
1098  else{
1099  if((Track = static_cast<AliAODTrack*>(PartCont->GetAcceptTrack(i)))){
1100  if (!Track) continue;
1101  if(TMath::Abs(Track->Eta())>0.9) continue;
1102  if (Track->Pt()<0.15) continue;
1103  if ((Track->Pt() >= minpT) && (Track->Pt()< maxpT)) {
1104  triggers[iTT] = i;
1105  iTT++;
1106  }
1107  }
1108  }
1109  }
1110 
1111 
1112  if (iTT == 0) return -99999;
1113  Int_t nbRn = 0, index = 0 ;
1114  TRandom3 random(0);
1115  nbRn = random.Integer(iTT);
1116  index = triggers[nbRn];
1117  //Printf("iTT Total= %d, nbRn = %d, Index = %d",iTT, nbRn, index );
1118  return index;
1119 
1120 }
1121 
1122 //__________________________________________________________________________________
1124 
1125  if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1126  else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1127  if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1128  else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1129  double dphi = mphi-vphi;
1130  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1131  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1132  return dphi;//dphi in [-Pi, Pi]
1133 }
1134 
1135 
1136 //_________________________________________________________________________
1138 
1139  std::vector<fastjet::PseudoJet> fInputVectors;
1140  fInputVectors.clear();
1141  fastjet::PseudoJet PseudoTracks;
1142 
1143  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
1144 
1145  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
1146  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
1147  if (!fTrk) continue;
1148  if(fDoTwoTrack==kTRUE && CheckClosePartner(i,fJet,fTrk,fTrackCont)) continue;
1149  PseudoTracks.reset(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
1150  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
1151  fInputVectors.push_back(PseudoTracks);
1152 
1153  }
1154  fastjet::JetAlgorithm jetalgo(fastjet::cambridge_algorithm);
1155 
1156 
1157 
1158  fastjet::JetDefinition fJetDef(jetalgo, 1., static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
1159 
1160  try {
1161  fastjet::ClusterSequence fClustSeqSA(fInputVectors, fJetDef);
1162  std::vector<fastjet::PseudoJet> fOutputJets;
1163  fOutputJets.clear();
1164  fOutputJets=fClustSeqSA.inclusive_jets(0);
1165 
1166  fastjet::PseudoJet jj;
1167  fastjet::PseudoJet j1;
1168  fastjet::PseudoJet j2;
1169  jj=fOutputJets[0];
1170  double ktaverage=0;
1171  double thetaverage=0;
1172  double nall=0;
1173  double flagSubjet=0;
1174  while(jj.has_parents(j1,j2)){
1175  nall=nall+1;
1176  if(j1.perp() < j2.perp()) swap(j1,j2);
1177  flagSubjet=0;
1178  double delta_R=j1.delta_R(j2);
1179  double z=j2.perp()/(j1.perp()+j2.perp());
1180  double y =log(1.0/delta_R);
1181  double lnpt_rel=log(j2.perp()*delta_R);
1182  double yh=j1.e()+j2.e();
1183  vector < fastjet::PseudoJet > constitj1 = sorted_by_pt(j1.constituents());
1184  if(constitj1[0].perp()>fMinPtConst) flagSubjet=1;
1185  if(z>fHardCutoff){
1186  ktaverage=ktaverage+lnpt_rel;
1187  thetaverage=thetaverage+delta_R;
1188  Double_t LundEntries[6] = {y,lnpt_rel,fOutputJets[0].perp(),nall,yh,flagSubjet};
1189  fHLundIterative->Fill(LundEntries);}
1190  jj=j1;}
1191 
1192 
1193  fShapesVar[2]=ktaverage/nall;
1194  fShapesVar[3]=thetaverage/nall;
1195 
1196 
1197  } catch (fastjet::Error) {
1198  AliError(" [w] FJ Exception caught.");
1199  //return -1;
1200  }
1201 
1202 
1203 
1204 
1205  return;
1206 
1207 
1208 }
1209 //_________________________________________________________________________
1211  AliJetContainer *jetCont = GetJetContainer(km);
1212  std::vector<fastjet::PseudoJet> fInputVectors;
1213  fInputVectors.clear();
1214  fastjet::PseudoJet PseudoTracks;
1215 
1216  AliParticleContainer *fTrackCont = jetCont->GetParticleContainer();
1217 
1218  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
1219  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
1220  if (!fTrk) continue;
1221 
1222  PseudoTracks.reset(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
1223  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
1224  fInputVectors.push_back(PseudoTracks);
1225 
1226  }
1227  fastjet::JetAlgorithm jetalgo(fastjet::cambridge_algorithm);
1228 
1229 
1230 
1231  fastjet::JetDefinition fJetDef(jetalgo, 1., static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
1232 
1233  try {
1234  fastjet::ClusterSequence fClustSeqSA(fInputVectors, fJetDef);
1235  std::vector<fastjet::PseudoJet> fOutputJets;
1236  fOutputJets.clear();
1237  fOutputJets=fClustSeqSA.inclusive_jets(0);
1238 
1239  fastjet::PseudoJet jj;
1240  fastjet::PseudoJet j1;
1241  fastjet::PseudoJet j2;
1242  jj=fOutputJets[0];
1243  double ktaverage=0;
1244  double thetaverage=0;
1245  double nall=0;
1246  double flagSubjet=0;
1247  while(jj.has_parents(j1,j2)){
1248  nall=nall+1;
1249  if(j1.perp() < j2.perp()) swap(j1,j2);
1250 
1251  double delta_R=j1.delta_R(j2);
1252 
1253  double lnpt_rel=log(j2.perp()*delta_R);
1254 
1255  ktaverage=ktaverage+lnpt_rel;
1256  thetaverage=thetaverage+delta_R;
1257 
1258  jj=j1;}
1259 
1260 
1261  average1=ktaverage/nall;
1262  average2=thetaverage/nall;
1263 
1264 
1265  } catch (fastjet::Error) {
1266  AliError(" [w] FJ Exception caught.");
1267  //return -1;
1268  }
1269 
1270 
1271 
1272 
1273  return;
1274 
1275 
1276  }
1277 
1278 //_________________________________________________________________________
1280 
1281  std::vector<fastjet::PseudoJet> fInputVectors;
1282  fInputVectors.clear();
1283  fastjet::PseudoJet PseudoTracks;
1284 
1285  std::vector<fastjet::PseudoJet> fInputVectorsM;
1286  fInputVectorsM.clear();
1287  fastjet::PseudoJet PseudoTracksM;
1288 
1289  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
1290  AliParticleContainer *fTrackContM = fJetContM->GetParticleContainer();
1291 
1292  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
1293  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
1294  if (!fTrk) continue;
1295  PseudoTracks.reset(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
1296  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
1297  fInputVectors.push_back(PseudoTracks);
1298 
1299  }
1300  fastjet::JetAlgorithm jetalgo(fastjet::cambridge_algorithm);
1301  fastjet::JetDefinition fJetDef(jetalgo, 1., static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
1302 
1303  if (fTrackContM) for (Int_t i=0; i<fJetM->GetNumberOfTracks(); i++) {
1304  AliVParticle *fTrk = fJetM->TrackAt(i, fTrackContM->GetArray());
1305  if (!fTrk) continue;
1306  PseudoTracksM.reset(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
1307  PseudoTracksM.set_user_index(fJetM->TrackAt(i)+100);
1308  fInputVectorsM.push_back(PseudoTracksM);
1309 
1310  }
1311  fastjet::JetAlgorithm jetalgoM(fastjet::cambridge_algorithm);
1312  fastjet::JetDefinition fJetDefM(jetalgoM, 1., static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
1313 
1314 
1315  try {
1316  fastjet::ClusterSequence fClustSeqSA(fInputVectors, fJetDef);
1317  std::vector<fastjet::PseudoJet> fOutputJets;
1318  fOutputJets.clear();
1319  fOutputJets=fClustSeqSA.inclusive_jets(0);
1320 
1321  fastjet::ClusterSequence fClustSeqSAM(fInputVectorsM, fJetDefM);
1322  std::vector<fastjet::PseudoJet> fOutputJetsM;
1323  fOutputJetsM.clear();
1324  fOutputJetsM=fClustSeqSAM.inclusive_jets(0);
1325 
1326 
1327 
1328  fastjet::PseudoJet jj,jjM;
1329  fastjet::PseudoJet j1,j1M;
1330  fastjet::PseudoJet j2,j2M;
1331  jj=fOutputJets[0];
1332  jjM=fOutputJetsM[0];
1333 
1334  double z1=0;
1335  double z2=0;
1336  double zcut=0.1;
1337  while((jj.has_parents(j1,j2)) && (z1<zcut)){
1338  if(j1.perp() < j2.perp()) swap(j1,j2);
1339 
1340  z1=j2.perp()/(j1.perp()+j2.perp());
1341  jj=j1;}
1342  if(z1<zcut) return;
1343 
1344 
1345 
1346  while((jjM.has_parents(j1M,j2M)) && (z2<zcut)){
1347  if(j1M.perp() < j2M.perp()) swap(j1M,j2M);
1348 
1349  z2=j2M.perp()/(j1M.perp()+j2M.perp());
1350  jjM=j1M;}
1351  if(z2<zcut) return;
1352 
1353 
1354 
1355  double delta_R1=j1.delta_R(j1M);
1356  double delta_R2=j2.delta_R(j2M);
1357  double delta_R=j1.delta_R(j2);
1358  double residz=(z1-z2)/z2;
1359  double resid1=(j1.perp()-j1M.perp())/j1M.perp();
1360  double resid2=(j2.perp()-j2M.perp())/j2M.perp();
1361 
1362  if((delta_R1<fSubjetCutoff) && (delta_R2<fSubjetCutoff)){
1363  Double_t ResolEntries[5] = {fOutputJets[0].perp(),delta_R,resid1,resid2,residz};
1364  fHCheckResolutionSubjets->Fill(ResolEntries);}
1365 
1366 
1367  } catch (fastjet::Error) {
1368  AliError(" [w] FJ Exception caught.");
1369  //return -1;
1370  }
1371 
1372 
1373 
1374 
1375  return;
1376 
1377 
1378 
1379 
1380 
1381 
1382 
1383 }
1384 
1385 
1387  //check if tracks are close//
1388  for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
1389  AliVParticle *fTrk2 = fJet->TrackAt(i, fTrackCont->GetArray());
1390  if (!fTrk2) continue;
1391  if(i==index) continue;
1392  Double_t phi1 = fTrk1->Phi();
1393  Double_t phi2 = fTrk2->Phi();
1394  Double_t chg1 = fTrk1->Charge();
1395  Double_t chg2 = fTrk2->Charge();
1396  Double_t ptv1 = fTrk1->Pt();
1397  Double_t ptv2 = fTrk2->Pt();
1398  Double_t deta=fTrk2->Eta()-fTrk1->Eta();
1399  const Float_t kLimit = fPhiCutValue* 3;
1400 
1401  if (TMath::Abs(fTrk1->Eta()-fTrk2->Eta()) < fEtaCutValue*2.5*3){
1402  Float_t initdpsinner = (phi2 - TMath::ASin(0.075*chg2*fMagFieldPolarity*0.8/ptv2) -
1403  (phi1 - TMath::ASin(0.075*chg1*fMagFieldPolarity*0.8/ptv1)));
1404 
1405  Float_t initdpsouter = (phi2 - TMath::ASin(0.075*chg2*fMagFieldPolarity*2.5/ptv2) -
1406  (phi1 - TMath::ASin(0.075*chg1*fMagFieldPolarity*2.5/ptv1)));
1407 
1408  initdpsinner = TVector2::Phi_mpi_pi(initdpsinner);
1409  initdpsouter = TVector2::Phi_mpi_pi(initdpsouter);
1410 
1411 
1412 
1413  if (TMath::Abs(initdpsinner) < kLimit ||
1414  TMath::Abs(initdpsouter) < kLimit || initdpsinner * initdpsouter < 0 ) {
1415  Double_t mindps = 1e5;
1416 
1417 
1418 
1419  for (Double_t rad = 0.8; rad < 2.51; rad += 0.01) {
1420  Double_t dps = (phi2 - TMath::ASin(0.075*chg2*fMagFieldPolarity*rad/ptv2) - (phi1 - TMath::ASin(0.075*chg1*fMagFieldPolarity*rad/ptv1)));
1421  dps = TVector2::Phi_mpi_pi(dps);
1422  if (TMath::Abs(dps) < TMath::Abs(mindps))
1423  mindps = dps;
1424  }
1425  if(TMath::Abs(mindps)<fPhiCutValue && TMath::Abs(deta)<fEtaCutValue) return kTRUE;
1426  } }
1427  return kFALSE;
1428  }}
1429 
1430 
1431 
1432 
1433 
1434 //________________________________________________________________________
1436  //
1437  // retrieve event objects
1438  //
1440  return kFALSE;
1441 
1442  return kTRUE;
1443 }
1444 
1445 //_______________________________________________________________________
1447 {
1448  // Called once at the end of the analysis.
1449 
1450  // fTreeObservableTagging = dynamic_cast<TTree*>(GetOutputData(1));
1451  // if (!fTreeObservableTagging){
1452  // Printf("ERROR: fTreeObservableTagging not available");
1453  // return;
1454  // }
1455 
1456 }
1457 
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:341
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
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
void RecursiveParentsMCAverage(AliEmcalJet *fJet, Int_t km, Double_t &aver1, Double_t &aver2)
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
void RecursiveParents(AliEmcalJet *fJet, AliJetContainer *fJetCont)
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:375
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