AliPhysics  d497547 (d497547)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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  fDerivSubtrOrder(0),
74  fh2ResponseUW(0x0),
75  fh2ResponseW(0x0),
76  fPhiJetCorr6(0x0),
77  fPhiJetCorr7(0x0),
78  fEtaJetCorr6(0x0),
79  fEtaJetCorr7(0x0),
80  fPtJetCorr(0x0),
81  fPtJet(0x0),
82  fhpTjetpT(0x0),
83  fhPt(0x0),
84  fhPhi(0x0),
85  fNbOfConstvspT(0x0),
86  fTreeObservableTagging(0)
87 
88 {
89  for(Int_t i=0;i<17;i++){
90  fShapesVar[i]=0;}
91  SetMakeGeneralHistograms(kTRUE);
92  DefineOutput(1, TList::Class());
93  DefineOutput(2, TTree::Class());
94 }
95 
96 //________________________________________________________________________
98  AliAnalysisTaskEmcalJet(name, kTRUE),
99  fContainer(0),
100  fMinFractionShared(0),
101  fJetShapeType(kData),
102  fJetShapeSub(kNoSub),
103  fJetSelection(kInclusive),
104  fPtThreshold(-9999.),
105  fRMatching(0.2),
106  fSelectedShapes(0),
107  fminpTTrig(20.),
108  fmaxpTTrig(50.),
109  fangWindowRecoil(0.6),
110  fSemigoodCorrect(0),
111  fHolePos(0),
112  fHoleWidth(0),
113  fCentSelectOn(kTRUE),
114  fCentMin(0),
115  fCentMax(10),
116  fOneConstSelectOn(kFALSE),
117  fDerivSubtrOrder(0),
118  fh2ResponseUW(0x0),
119  fh2ResponseW(0x0),
120  fPhiJetCorr6(0x0),
121  fPhiJetCorr7(0x0),
122  fEtaJetCorr6(0x0),
123  fEtaJetCorr7(0x0),
124  fPtJetCorr(0x0),
125  fPtJet(0x0),
126  fhpTjetpT(0x0),
127  fhPt(0x0),
128  fhPhi(0x0),
129  fNbOfConstvspT(0x0),
130  fTreeObservableTagging(0)
131 
132 {
133  // Standard constructor.
134  for(Int_t i=0;i<17;i++){
135  fShapesVar[i]=0;}
137 
138  DefineOutput(1, TList::Class());
139  DefineOutput(2, TTree::Class());
140 }
141 
142 //________________________________________________________________________
144 {
145  // Destructor.
146 }
147 
148 //________________________________________________________________________
150 {
151  // Create user output.
152 
154 
155  Bool_t oldStatus = TH1::AddDirectoryStatus();
156  TH1::AddDirectory(kFALSE);
157 
158 
159 
160 
161  fh2ResponseUW= new TH2F("fh2ResponseUW", "fh2ResponseUW", 100, 0, 200, 100, 0, 200);
162  fOutput->Add(fh2ResponseUW);
163  fh2ResponseW= new TH2F("fh2ResponseW", "fh2ResponseW", 100, 0, 200, 100, 0, 200);
164  fOutput->Add(fh2ResponseW);
165  fPhiJetCorr6= new TH2F("fPhiJetCorr6", "fPhiJetCorr6", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
166  fOutput->Add(fPhiJetCorr6);
167  fEtaJetCorr6= new TH2F("fEtaJetCorr6", "fEtaJetCorr6", 50, -1.5, 1.5, 50, -1.5, 1.5);
168  fOutput->Add(fEtaJetCorr6);
169 
170  fPhiJetCorr7= new TH2F("fPhiJetCorr7", "fPhiJetCorr7", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
171  fOutput->Add(fPhiJetCorr7);
172  fEtaJetCorr7= new TH2F("fEtaJetCorr7", "fEtaJetCorr7", 50, -1.5, 1.5, 50, -1.5, 1.5);
173  fOutput->Add(fEtaJetCorr7);
174 
175  fPtJetCorr= new TH2F("fPtJetCorr", "fPtJetCorr", 100, 0, 200, 100, 0, 200);
176  fOutput->Add(fPtJetCorr);
177  fPtJet= new TH1F("fPtJet", "fPtJet", 100, 0, 200);
178  fOutput->Add(fPtJet);
179 
180  fhpTjetpT= new TH2F("fhpTjetpT", "fhpTjetpT", 200, 0, 200, 200, 0, 200);
181  fOutput->Add(fhpTjetpT);
182  fhPt= new TH1F("fhPt", "fhPt", 200, 0, 200);
183  fOutput->Add(fhPt);
184  fhPhi= new TH1F("fhPhi", "fhPhi", 100, -TMath::Pi(), TMath::Pi());
185  fOutput->Add(fhPhi);
186 
187  fNbOfConstvspT=new TH2F("fNbOfConstvspT", "fNbOfConstvspT", 100, 0, 100, 200, 0, 200);
188  fOutput->Add(fNbOfConstvspT);
189 
190  // =========== Switch on Sumw2 for all histos ===========
191  for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
192  TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
193  if (h1){
194  h1->Sumw2();
195  continue;
196  }
197  THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
198  if(hn)hn->Sumw2();
199  }
200 
201 
202  TH1::AddDirectory(oldStatus);
203  const Int_t nVar = 19;
204  const char* nameoutput = GetOutputSlot(2)->GetContainer()->GetName();
205  fTreeObservableTagging = new TTree(nameoutput, nameoutput);
206 
207 
208  TString *fShapesVarNames = new TString [nVar];
209 
210  fShapesVarNames[0] = "partonCode";
211  fShapesVarNames[1] = "ptJet";
212  fShapesVarNames[2] = "ptDJet";
213  fShapesVarNames[3] = "mJet";
214  // fShapesVarNames[4] = "nbOfConst";
215  fShapesVarNames[4] = "angularity";
216  fShapesVarNames[5] = "circularity";
217  fShapesVarNames[6] = "lesub";
218  fShapesVarNames[7] = "coronna";
219 
220  fShapesVarNames[8] = "ptJetMatch";
221  fShapesVarNames[9] = "ptDJetMatch";
222  fShapesVarNames[10] = "mJetMatch";
223  // fShapesVarNames[12] = "nbOfConstMatch";
224  fShapesVarNames[11] = "angularityMatch";
225  fShapesVarNames[12] = "circularityMatch";
226  fShapesVarNames[13] = "lesubMatch";
227  fShapesVarNames[14] = "coronnaMatch";
228  fShapesVarNames[15]="weightPythia";
229  //fShapesVarNames[14]="ntrksEvt";
230  fShapesVarNames[16]="rhoVal";
231  fShapesVarNames[17]="rhoMassVal";
232  fShapesVarNames[18]="ptUnsub";
233 
234  for(Int_t ivar=0; ivar < nVar; ivar++){
235  cout<<"looping over variables"<<endl;
236  fTreeObservableTagging->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/F", fShapesVarNames[ivar].Data()));}
237 
238  PostData(1,fOutput);
239  PostData(2,fTreeObservableTagging);
240 
241  delete [] fShapesVarNames;
242 }
243 
244 //________________________________________________________________________
246 {
247  // Run analysis code here, if needed. It will be executed before FillHistograms().
248 
249  return kTRUE;
250 }
251 
252 //________________________________________________________________________
254 {
255  // Fill histograms.
256  //cout<<"base container"<<endl;
257  AliEmcalJet* jet1 = NULL;
258  AliJetContainer *jetCont = GetJetContainer(0);
259 
260  Float_t kWeight=1;
261  if (fCentSelectOn)
262  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
263 
264  AliAODTrack *triggerHadron = 0x0;
265 
266  if (fJetSelection == kRecoil) {
267  //Printf("Recoil jets!!!, fminpTTrig = %f, fmaxpTTrig = %f", fminpTTrig, fmaxpTTrig);
268  Int_t triggerHadronLabel = SelectTrigger(fminpTTrig, fmaxpTTrig);
269 
270 
271  if (triggerHadronLabel==-99999) {
272  //Printf ("Trigger Hadron not found, return");
273  return 0;}
274 
275  AliTrackContainer *PartCont =NULL;
276  AliParticleContainer *PartContMC=NULL;
277 
278  if (fJetShapeSub==kConstSub){
280  else PartCont = GetTrackContainer(1);
281  }
282  else{
284  else PartCont = GetTrackContainer(0);
285  }
286  TClonesArray *TrackArray = NULL;
287  TClonesArray *TrackArrayMC = NULL;
288  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) TrackArrayMC = PartContMC->GetArray();
289  else TrackArray = PartCont->GetArray();
290  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) triggerHadron = static_cast<AliAODTrack*>(TrackArrayMC->At(triggerHadronLabel));
291  else triggerHadron = static_cast<AliAODTrack*>(TrackArray->At(triggerHadronLabel));
292 
293 
294 
295 
296 
297  if (!triggerHadron) {
298  //Printf("No Trigger hadron with the found label!!");
299  return 0;
300  }
301 
302  if(fSemigoodCorrect){
303  Double_t disthole=RelativePhi(triggerHadron->Phi(),fHolePos);
304  if(TMath::Abs(disthole)+fHoleWidth>TMath::Pi()-fangWindowRecoil){
305  return 0;}
306  }
307 
308  fhPt->Fill(triggerHadron->Pt());
309 
310  }
311 
312 
314  TClonesArray *trackArrayAn = partContAn->GetArray();
315  Int_t ntracksEvt = trackArrayAn->GetEntriesFast();
316 
317  Float_t rhoVal=0, rhoMassVal = 0.;
318  if(jetCont) {
319 
320  jetCont->ResetCurrentID();
322  //rho
323  AliRhoParameter* rhoParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject("RhoSparseR020"));
324  if (!rhoParam) {
325  Printf("%s: Could not retrieve rho %s (some histograms will be filled with zero)!", GetName(), jetCont->GetRhoName().Data());
326  } else rhoVal = rhoParam->GetVal();
327  //rhom
328  AliRhoParameter* rhomParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject("RhoMassSparseR020"));
329  if (!rhomParam) {
330  Printf("%s: Could not retrieve rho_m %s (some histograms will be filled with zero)!", GetName(), jetCont->GetRhoMassName().Data());
331  } else rhoMassVal = rhomParam->GetVal();
332  }
333 
334  while((jet1 = jetCont->GetNextAcceptJet())) {
335  if (!jet1) continue;
336  AliEmcalJet* jet2 = 0x0;
337  AliEmcalJet* jet3 = 0x0;
338  fPtJet->Fill(jet1->Pt());
339  AliEmcalJet *jetUS = NULL;
340  Int_t ifound=0, jfound=0;
341  Int_t ilab=-1, jlab=-1;
342 
344  Double_t disthole=RelativePhi(jet1->Phi(),fHolePos);
345  if(TMath::Abs(disthole)<fHoleWidth){
346  continue;}
347  }
348 
349  Float_t dphiRecoil = 0.;
350  if (fJetSelection == kRecoil){
351  dphiRecoil = RelativePhi(triggerHadron->Phi(), jet1->Phi());
352  if (TMath::Abs(dphiRecoil) < (TMath::Pi() - fangWindowRecoil)) {
353  // Printf("Recoil jets back to back not found! continuing");
354  continue;
355  }
356 
357  fhpTjetpT->Fill(triggerHadron->Pt(), jet1->Pt());
358  //Printf(" ************ FILLING HISTOS****** shapeSub = %d, triggerHadron = %f, jet1 = %f", fJetShapeSub, triggerHadron->Pt(), jet1->Pt());
359  fhPhi->Fill(RelativePhi(triggerHadron->Phi(), jet1->Phi()));
360 
361  }
362 
363 
364  fShapesVar[0] = 0.;
366  AliJetContainer *jetContTrue = GetJetContainer(1);
367  AliJetContainer *jetContUS = GetJetContainer(2);
368 
369  if(fJetShapeSub==kConstSub){
370  for(Int_t i = 0; i<jetContUS->GetNJets(); i++) {
371  jetUS = jetContUS->GetJet(i);
372  if(jetUS->GetLabel()==jet1->GetLabel()) {
373  ifound++;
374  if(ifound==1) ilab = i;
375  }
376  }
377  if(ilab==-1) continue;
378  jetUS=jetContUS->GetJet(ilab);
379  jet2=jetUS->ClosestJet();
380  }
381 
382  if(!(fJetShapeSub==kConstSub)) jet2 = jet1->ClosestJet();
383  if (!jet2) {
384  Printf("jet2 does not exist, returning");
385  continue;
386  }
387 
388  AliJetContainer *jetContPart=GetJetContainer(3);
389  jet3=jet2->ClosestJet();
390 
391  if(!jet3){
392  Printf("jet3 does not exist, returning");
393  continue;
394  }
395  cout<<"jet 3 exists"<<jet3->Pt()<<endl;
396 
397 
398  fh2ResponseUW->Fill(jet1->Pt(),jet2->Pt());
399 
400  Double_t fraction=0;
401  if(!(fJetShapeSub==kConstSub)) fraction = jetCont->GetFractionSharedPt(jet1);
402  if(fJetShapeSub==kConstSub) fraction = jetContUS->GetFractionSharedPt(jetUS);
403  //if (fraction > 0.1) cout<<"***** hey a jet matched with fraction"<<fraction<<" "<<jet1->Pt()<<" "<<jet2->Pt()<<" "<<fCent<<endl;
404 
405  if(fraction<fMinFractionShared) continue;
406  //InputEvent()->Print();
407 
408  }
409 
410 
411 
412  if (fJetShapeType == kPythiaDef){
413 
414  AliJetContainer *jetContTrue = GetJetContainer(1);
415  AliJetContainer *jetContUS = GetJetContainer(2);
416  AliJetContainer *jetContPart = GetJetContainer(3);
417 
418  if(fJetShapeSub==kConstSub){
419 
420  for(Int_t i = 0; i<jetContUS->GetNJets(); i++) {
421  jetUS = jetContUS->GetJet(i);
422  if(jetUS->GetLabel()==jet1->GetLabel()) {
423  ifound++;
424  if(ifound==1) ilab = i;
425  }
426  }
427  if(ilab==-1) continue;
428  jetUS=jetContUS->GetJet(ilab);
429  jet2=jetUS->ClosestJet();
430 
431  if (!jet2) {
432  Printf("jet2 does not exist, returning");
433  continue;
434  }
435 
436  for(Int_t j=0; j<jetContPart->GetNJets(); j++) {
437 
438  jet3 = jetContPart->GetJet(j);
439  if(!jet3) continue;
440  if(jet3->GetLabel()==jet2->GetLabel()) {
441  jfound++;
442  if(jfound==1) jlab = j;
443  }
444  }
445  if(jlab==-1) continue;
446  jet3=jetContPart->GetJet(jlab);
447  if(!jet3){
448  Printf("jet3 does not exist, returning");
449  continue;
450  }
451  }
452  if(!(fJetShapeSub==kConstSub)) jet3 = jet1->ClosestJet();
453  if (!jet3) {
454  Printf("jet3 does not exist, returning");
455  continue;
456  }
457 
458 
459  fh2ResponseUW->Fill(jet1->Pt(),jet3->Pt());
460 
461 
462  }
463 
464 
465  if (fJetShapeType == kGenOnTheFly){
466  const AliEmcalPythiaInfo *partonsInfo = 0x0;
467  partonsInfo = GetPythiaInfo();
468  Double_t jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi6());
469  Double_t detap1=(jet1->Eta())-(partonsInfo->GetPartonEta6());
470  kWeight=partonsInfo->GetPythiaEventWeight();
471  fh2ResponseW->Fill(jet1->Pt(),jet1->Pt(),kWeight);
472 
473  Float_t dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
474  fEtaJetCorr6->Fill(jet1->Eta(), partonsInfo->GetPartonEta6());
475  fPhiJetCorr6->Fill(jet1->Phi(), partonsInfo->GetPartonPhi6());
476  if(dRp1 < fRMatching) {
477  fShapesVar[0] = partonsInfo->GetPartonFlag6();
478  fPtJetCorr ->Fill(partonsInfo->GetPartonPt6(), jet1->Pt());
479  }
480  else {
481  jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi7());
482  detap1=(jet1->Eta())-(partonsInfo->GetPartonEta7());
483  dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
484  fEtaJetCorr7->Fill(jet1->Eta(), partonsInfo->GetPartonEta7());
485  fPhiJetCorr7->Fill(jet1->Phi(), partonsInfo->GetPartonPhi7());
486  if(dRp1 < fRMatching) {
487  fShapesVar[0] = partonsInfo->GetPartonFlag7();
488  fPtJetCorr->Fill(partonsInfo->GetPartonPt7(), jet1->Pt());
489  }
490  else fShapesVar[0]=0;
491  }
492  }
493 
494 
495 
496 
497  Double_t ptSubtracted = 0;
498  if (fJetShapeSub==kConstSub) ptSubtracted= jet1->Pt();
499 
500  else if (fJetShapeSub==kDerivSub) {
501  ptSubtracted=jet1->Pt()-GetRhoVal(0)*jet1->Area();
502  }
503 
504  else if (fJetShapeSub==kNoSub) {
505  if ((fJetShapeType==kData) || (fJetShapeType==kDetEmbPartPythia)) ptSubtracted=jet1->Pt()-GetRhoVal(0)*jet1->Area();
506  else if ((fJetShapeType==kPythiaDef) || (fJetShapeType==kMCTrue) || (fJetShapeType==kGenOnTheFly)) ptSubtracted= jet1->Pt();
507  }
508 
509  //Printf("ptSubtracted=%f,fPtThreshold =%f ", ptSubtracted, fPtThreshold);
510  if (ptSubtracted < fPtThreshold) continue;
511 
512  if (fOneConstSelectOn == kTRUE) fNbOfConstvspT->Fill(GetJetNumberOfConstituents(jet1,0), ptSubtracted);
513  if ((fCentSelectOn == kFALSE) && (jet1->GetNumberOfTracks() <= 1)) continue;
514 
515 
516  fShapesVar[1] = ptSubtracted;
517  fShapesVar[2] = GetJetpTD(jet1,0);
518  fShapesVar[3] = GetJetMass(jet1,0);
519  fShapesVar[4] = GetJetAngularity(jet1,0);
520  fShapesVar[5] = GetJetCircularity(jet1,0);
521  fShapesVar[6] = GetJetLeSub(jet1,0);
522  fShapesVar[6] = GetJetCoronna(jet1,0);
523 
524  Float_t ptMatch=0., ptDMatch=0., massMatch=0., constMatch=0.,angulMatch=0.,circMatch=0., lesubMatch=0., sigma2Match=0., coronnaMatch=0;
525  Int_t kMatched = 0;
526 
527  if (fJetShapeType==kPythiaDef) {
528  kMatched =1;
529  if(fJetShapeSub==kConstSub) kMatched = 3;
530 
531  ptMatch=jet3->Pt();
532  ptDMatch=GetJetpTD(jet3, kMatched);
533  massMatch=GetJetMass(jet3,kMatched);
534  //constMatch=1.*GetJetNumberOfConstituents(jet2,kMatched);
535  angulMatch=GetJetAngularity(jet3, kMatched);
536  circMatch=GetJetCircularity(jet3, kMatched);
537  lesubMatch=GetJetLeSub(jet3, kMatched);
538  coronnaMatch=GetJetCoronna(jet3,kMatched);
539  //sigma2Match = GetSigma2(jet2, kMatched);
540  }
541 
543  if(fJetShapeSub==kConstSub) kMatched = 3;
544  if(fJetShapeSub==kDerivSub) kMatched = 2;
545  ptMatch=jet3->Pt();
546  ptDMatch=GetJetpTD(jet3, kMatched);
547  massMatch=GetJetMass(jet3,kMatched);
548  // constMatch=1.*GetJetNumberOfConstituents(jet3,kMatched);
549  angulMatch=GetJetAngularity(jet3, kMatched);
550  circMatch=GetJetCircularity(jet3, kMatched);
551  lesubMatch=GetJetLeSub(jet3, kMatched);
552  coronnaMatch = GetJetCoronna(jet3, kMatched);
553 
554  }
555 
556 
557 
559  kMatched = 0;
560  ptMatch=0.;
561  ptDMatch=0.;
562  massMatch=0.;
563  //constMatch=0.;
564  angulMatch=0.;
565  circMatch=0.;
566  lesubMatch=0.;
567  coronnaMatch =0.;
568 
569  }
570 
571 
572 
573  fShapesVar[8] = ptMatch;
574  fShapesVar[9] = ptDMatch;
575  fShapesVar[10] = massMatch;
576  fShapesVar[11] = angulMatch;
577  fShapesVar[12] = circMatch;
578  fShapesVar[13] = lesubMatch;
579  fShapesVar[14] = coronnaMatch;
580  fShapesVar[15] = kWeight;
581  //fShapesVar[16] = ntracksEvt;
582  fShapesVar[16] = rhoVal;
583  fShapesVar[17] = rhoMassVal;
584  fShapesVar[18] = jet1->Pt();
585 
586 
587  fTreeObservableTagging->Fill();
588 
589 
590 
591 
592 
593 
594  }
595 
596  }
597 
598  return kTRUE;
599 }
600 
601 //________________________________________________________________________
603  //calc subtracted jet mass
604  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
606  else return jet->GetShapeProperties()->GetSecondOrderSubtracted();
607  else
608  return jet->M();
609 }
610 
611 //________________________________________________________________________
613 
614  AliJetContainer *jetCont = GetJetContainer(jetContNb);
615  if (!jet->GetNumberOfTracks())
616  return 0;
617  Double_t den=0.;
618  Double_t num = 0.;
619  AliVParticle *vp1 = 0x0;
620  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
621  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
622 
623  if (!vp1){
624  Printf("AliVParticle associated to constituent not found");
625  continue;
626  }
627 
628  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
629  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
630  Double_t dr = TMath::Sqrt(dr2);
631  num=num+vp1->Pt()*dr;
632  den=den+vp1->Pt();
633  }
634  return num/den;
635 }
636 
637 //________________________________________________________________________
639 
640  if((fJetShapeSub==kDerivSub) && (jetContNb==0))
643  else
644  return Angularity(jet, jetContNb);
645 
646 }
647 
648 //____________________________________________________________________________
649 
651 
652  AliTrackContainer *PartCont = NULL;
653  AliParticleContainer *PartContMC = NULL;
654 
655 
656  if (fJetShapeSub==kConstSub){
658  else PartCont = GetTrackContainer(1);
659  }
660  else{
662  else PartCont = GetTrackContainer(0);
663  }
664 
665  TClonesArray *TracksArray = NULL;
666  TClonesArray *TracksArrayMC = NULL;
667 
668  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) TracksArrayMC = PartContMC->GetArray();
669  else TracksArray = PartCont->GetArray();
670 
672  if(!PartContMC || !TracksArrayMC) return -2;
673  }
674  else {
675  if(!PartCont || !TracksArray) return -2;
676  }
677 
678 
679  AliAODTrack *Track = 0x0;
680  Float_t sumpt=0;
681  Int_t NTracks=0;
682  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) NTracks = TracksArrayMC->GetEntriesFast();
683  else NTracks = TracksArray->GetEntriesFast();
684 
685  for(Int_t i=0; i < NTracks; i++){
687  if((Track = static_cast<AliAODTrack*>(PartContMC->GetAcceptParticle(i)))){
688  if (!Track) continue;
689  if(TMath::Abs(Track->Eta())>0.9) continue;
690  Double_t dphi = RelativePhi(Track->Phi(),jet->Phi());
691  Double_t dr2 = (Track->Eta()-jet->Eta())*(Track->Eta()-jet->Eta()) + dphi*dphi;
692  Double_t dr = TMath::Sqrt(dr2);
693  if((dr>=0.8) && (dr<1)) sumpt=sumpt+Track->Pt();
694  }
695  }
696  else{
697  if((Track = static_cast<AliAODTrack*>(PartCont->GetAcceptTrack(i)))){
698  if (!Track) continue;
699  if(TMath::Abs(Track->Eta())>0.9) continue;
700  if (Track->Pt()<0.15) continue;
701  Double_t dphi = RelativePhi(Track->Phi(),jet->Phi());
702  Double_t dr2 = (Track->Eta()-jet->Eta())*(Track->Eta()-jet->Eta()) + dphi*dphi;
703  Double_t dr = TMath::Sqrt(dr2);
704  if((dr>=0.8) && (dr<1)) sumpt=sumpt+Track->Pt();
705 
706  }
707  }
708  }
709 
710 
711 
712  return sumpt;
713 
714 
715 
716 
717 }
718 
719 //________________________________________________________________________
721 
722  if((fJetShapeSub==kDerivSub) && (jetContNb==0)) return -2;
723  else
724  return Coronna(jet, jetContNb);
725 
726 }
727 
728 
729 
730 
731 
732 //________________________________________________________________________
734 
735  AliJetContainer *jetCont = GetJetContainer(jetContNb);
736  if (!jet->GetNumberOfTracks())
737  return 0;
738  Double_t den=0.;
739  Double_t num = 0.;
740  AliVParticle *vp1 = 0x0;
741  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
742  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
743 
744  if (!vp1){
745  Printf("AliVParticle associated to constituent not found");
746  continue;
747  }
748 
749  num=num+vp1->Pt()*vp1->Pt();
750  den=den+vp1->Pt();
751  }
752  return TMath::Sqrt(num)/den;
753 }
754 
755 //________________________________________________________________________
757  //calc subtracted jet mass
758  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
760  else return jet->GetShapeProperties()->GetSecondOrderSubtractedpTD();
761  else
762  return PTD(jet, jetContNb);
763 
764 }
765 
766 //_____________________________________________________________________________
768 
769  AliJetContainer *jetCont = GetJetContainer(jetContNb);
770  if (!jet->GetNumberOfTracks())
771  return 0;
772  Double_t mxx = 0.;
773  Double_t myy = 0.;
774  Double_t mxy = 0.;
775  int nc = 0;
776  Double_t sump2 = 0.;
777  Double_t pxjet=jet->Px();
778  Double_t pyjet=jet->Py();
779  Double_t pzjet=jet->Pz();
780 
781 
782  //2 general normalized vectors perpendicular to the jet
783  TVector3 ppJ1(pxjet, pyjet, pzjet);
784  TVector3 ppJ3(- pxjet* pzjet, - pyjet * pzjet, pxjet * pxjet + pyjet * pyjet);
785  ppJ3.SetMag(1.);
786  TVector3 ppJ2(-pyjet, pxjet, 0);
787  ppJ2.SetMag(1.);
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  TVector3 pp(vp1->Px(), vp1->Py(), vp1->Pz());
798 
799  //local frame
800  TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
801  TVector3 pPerp = pp - pLong;
802  //projection onto the two perpendicular vectors defined above
803 
804  Float_t ppjX = pPerp.Dot(ppJ2);
805  Float_t ppjY = pPerp.Dot(ppJ3);
806  Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
807  if(ppjT<=0) return 0;
808 
809  mxx += (ppjX * ppjX / ppjT);
810  myy += (ppjY * ppjY / ppjT);
811  mxy += (ppjX * ppjY / ppjT);
812  nc++;
813  sump2 += ppjT;}
814 
815  if(nc<2) return 0;
816  if(sump2==0) return 0;
817  // Sphericity Matrix
818  Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
819  TMatrixDSym m0(2,ele);
820 
821  // Find eigenvectors
822  TMatrixDSymEigen m(m0);
823  TVectorD eval(2);
824  TMatrixD evecm = m.GetEigenVectors();
825  eval = m.GetEigenValues();
826  // Largest eigenvector
827  int jev = 0;
828  // cout<<eval[0]<<" "<<eval[1]<<endl;
829  if (eval[0] < eval[1]) jev = 1;
830  TVectorD evec0(2);
831  // Principle axis
832  evec0 = TMatrixDColumn(evecm, jev);
833  Double_t compx=evec0[0];
834  Double_t compy=evec0[1];
835  TVector2 evec(compx, compy);
836  Double_t circ=0;
837  if(jev==1) circ=2*eval[0];
838  if(jev==0) circ=2*eval[1];
839 
840  return circ;
841 
842 
843 
844 }
845 
846 
847 
848 
849 //________________________________________________________________________
851  //calc subtracted jet mass
852 
853  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
856  else
857  return Circularity(jet, jetContNb);
858 
859 }
860 
861 //________________________________________________________________________
863 
864  AliJetContainer *jetCont = GetJetContainer(jetContNb);
865  if (!jet->GetNumberOfTracks())
866  return 0;
867  Double_t den=0.;
868  Double_t num = 0.;
869  AliVParticle *vp1 = 0x0;
870  AliVParticle *vp2 = 0x0;
871  std::vector<int> ordindex;
872  ordindex=jet->GetPtSortedTrackConstituentIndexes(jetCont->GetParticleContainer()->GetArray());
873  //Printf("Nbof const = %d", jet->GetNumberOfTracks());
874  //Printf("ordindex[0] = %d, ordindex[1] = %d", ordindex[0], ordindex[1]);
875 
876  if(ordindex.size()<2) return -1;
877 
878  vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[0], jetCont->GetParticleContainer()->GetArray()));
879  if (!vp1){
880  Printf("AliVParticle associated to Leading constituent not found");
881  return -1;
882  }
883 
884  vp2 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[1], jetCont->GetParticleContainer()->GetArray()));
885  if (!vp2){
886  Printf("AliVParticle associated to Subleading constituent not found");
887  return -1;
888  }
889 
890 
891  num=vp1->Pt();
892  den=vp2->Pt();
893  //Printf("vp1->Pt() =%f, vp2->Pt() =%f", vp1->Pt(), vp2->Pt());
894 
895 return num-den;
896 }
897 
898 //________________________________________________________________________
900  //calc subtracted jet mass
901 
902  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
905  else
906  return LeSub(jet, jetContNb);
907 
908 }
909 
910 //________________________________________________________________________
912  //calc subtracted jet mass
913 
914  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
917  else
918  return jet->GetNumberOfTracks();
919 
920  }
921 
922 
923 //______________________________________________________________________________
925 
926  AliJetContainer *jetCont = GetJetContainer(jetContNb);
927  if (!jet->GetNumberOfTracks())
928  return 0;
929  Double_t mxx = 0.;
930  Double_t myy = 0.;
931  Double_t mxy = 0.;
932  int nc = 0;
933  Double_t sump2 = 0.;
934 
935  AliVParticle *vp1 = 0x0;
936  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
937  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
938 
939  if (!vp1){
940  Printf("AliVParticle associated to constituent not found");
941  continue;
942  }
943 
944  Double_t ppt=vp1->Pt();
945  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
946 
947  Double_t deta = vp1->Eta()-jet->Eta();
948  mxx += ppt*ppt*deta*deta;
949  myy += ppt*ppt*dphi*dphi;
950  mxy -= ppt*ppt*deta*TMath::Abs(dphi);
951  nc++;
952  sump2 += ppt*ppt;
953 
954  }
955  if(nc<2) return 0;
956  if(sump2==0) return 0;
957  // Sphericity Matrix
958  Double_t ele[4] = {mxx , mxy , mxy , myy };
959  TMatrixDSym m0(2,ele);
960 
961  // Find eigenvectors
962  TMatrixDSymEigen m(m0);
963  TVectorD eval(2);
964  TMatrixD evecm = m.GetEigenVectors();
965  eval = m.GetEigenValues();
966  // Largest eigenvector
967  int jev = 0;
968  // cout<<eval[0]<<" "<<eval[1]<<endl;
969  if (eval[0] < eval[1]) jev = 1;
970  TVectorD evec0(2);
971  // Principle axis
972  evec0 = TMatrixDColumn(evecm, jev);
973  Double_t compx=evec0[0];
974  Double_t compy=evec0[1];
975  TVector2 evec(compx, compy);
976  Double_t sig=0;
977  if(jev==1) sig=TMath::Sqrt(TMath::Abs(eval[0])/sump2);
978  if(jev==0) sig=TMath::Sqrt(TMath::Abs(eval[1])/sump2);
979 
980  return sig;
981 
982 }
983 
984 //________________________________________________________________________
986  //calc subtracted jet mass
987 
988  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
991  else
992  return Sigma2(jet, jetContNb);
993 
994 }
995 
996 //________________________________________________________________________
998 
999  AliTrackContainer *PartCont = NULL;
1000  AliParticleContainer *PartContMC = NULL;
1001 
1002 
1003  if (fJetShapeSub==kConstSub){
1005  else PartCont = GetTrackContainer(1);
1006  }
1007  else{
1009  else PartCont = GetTrackContainer(0);
1010  }
1011 
1012  TClonesArray *TracksArray = NULL;
1013  TClonesArray *TracksArrayMC = NULL;
1014 
1015  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) TracksArrayMC = PartContMC->GetArray();
1016  else TracksArray = PartCont->GetArray();
1017 
1019  if(!PartContMC || !TracksArrayMC) return -99999;
1020  }
1021  else {
1022  if(!PartCont || !TracksArray) return -99999;
1023  }
1024 
1025 
1026  AliAODTrack *Track = 0x0;
1027 
1028 
1029 
1030  TList *trackList = new TList();
1031  Int_t triggers[100];
1032  for (Int_t iTrigger=0; iTrigger<100; iTrigger++) triggers[iTrigger] = 0;
1033  Int_t iTT = 0;
1034  Int_t NTracks=0;
1035  if (fJetShapeType == AliAnalysisTaskEmcalQGTagging::kGenOnTheFly) NTracks = TracksArrayMC->GetEntriesFast();
1036  else NTracks = TracksArray->GetEntriesFast();
1037 
1038  for(Int_t i=0; i < NTracks; i++){
1040  if((Track = static_cast<AliAODTrack*>(PartContMC->GetAcceptParticle(i)))){
1041  if (!Track) continue;
1042  if(TMath::Abs(Track->Eta())>0.9) continue;
1043  if (Track->Pt()<0.15) continue;
1044  if ((Track->Pt() >= minpT) && (Track->Pt()< maxpT)) {
1045  triggers[iTT] = i;
1046  iTT++;
1047  }
1048  }
1049  }
1050  else{
1051  if((Track = static_cast<AliAODTrack*>(PartCont->GetAcceptTrack(i)))){
1052  if (!Track) continue;
1053  if(TMath::Abs(Track->Eta())>0.9) continue;
1054  if (Track->Pt()<0.15) continue;
1055  if ((Track->Pt() >= minpT) && (Track->Pt()< maxpT)) {
1056  triggers[iTT] = i;
1057  iTT++;
1058  }
1059  }
1060  }
1061  }
1062 
1063 
1064  if (iTT == 0) return -99999;
1065  Int_t nbRn = 0, index = 0 ;
1066  TRandom3* random = new TRandom3(0);
1067  nbRn = random->Integer(iTT);
1068  index = triggers[nbRn];
1069  //Printf("iTT Total= %d, nbRn = %d, Index = %d",iTT, nbRn, index );
1070  return index;
1071 
1072 }
1073 
1074 //__________________________________________________________________________________
1076 
1077  if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1078  else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1079  if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1080  else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1081  double dphi = mphi-vphi;
1082  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1083  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1084  return dphi;//dphi in [-Pi, Pi]
1085 }
1086 
1087 
1088 //________________________________________________________________________
1090  //
1091  // retrieve event objects
1092  //
1094  return kFALSE;
1095 
1096  return kTRUE;
1097 }
1098 
1099 //_______________________________________________________________________
1101 {
1102  // Called once at the end of the analysis.
1103 
1104  // fTreeObservableTagging = dynamic_cast<TTree*>(GetOutputData(1));
1105  // if (!fTreeObservableTagging){
1106  // Printf("ERROR: fTreeObservableTagging not available");
1107  // return;
1108  // }
1109 
1110 }
1111 
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:123
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:222
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:114
Double_t Py() const
Definition: AliEmcalJet.h:100
Double_t Phi() const
Definition: AliEmcalJet.h:110
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:117
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:153
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:132
Float_t GetJetCoronna(AliEmcalJet *jet, Int_t jetContNb)
TString kData
Declare data MC or deltaAOD.
Double_t Px() const
Definition: AliEmcalJet.h:99
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:102
Double_t GetRhoVal(Int_t i=0) const
Double_t GetFirstOrderSubtractedCircularity() const
AliEmcalList * fOutput
!output list
Float_t Sigma2(AliEmcalJet *jet, Int_t jetContNb)
Double_t RelativePhi(Double_t mphi, Double_t vphi)
const Int_t nVar
Double_t GetSecondOrderSubtractedCircularity() const
Float_t GetJetNumberOfConstituents(AliEmcalJet *jet, Int_t jetContNb)
AliTrackContainer * GetTrackContainer(Int_t i=0) const
Store some informaion about a Pythia eventThis class is used to store some information about a Pythia...
void SetMakeGeneralHistograms(Bool_t g)
virtual AliVTrack * GetAcceptTrack(Int_t i=-1) const
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
const AliEmcalPythiaInfo * GetPythiaInfo() const
Float_t GetJetpTD(AliEmcalJet *jet, Int_t jetContNb)
Double_t Pz() const
Definition: AliEmcalJet.h:101
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:256
Float_t GetJetLeSub(AliEmcalJet *jet, Int_t jetContNb)
Float_t GetPartonPt7() const
Double_t M() const
Definition: AliEmcalJet.h:113
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
AliEmcalJet * GetJet(Int_t i) const