]>
Commit | Line | Data |
---|---|---|
cfa7b4df | 1 | // |
2 | // Jet QG tagging analysis task. | |
3 | // | |
4 | // Author: D. Caffarri, L. Cunqueiro | |
5 | ||
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 "TMatrixD.h" | |
20 | #include "TMatrixDSym.h" | |
21 | #include "TMatrixDSymEigen.h" | |
22 | #include "TVector3.h" | |
23 | #include "TVector2.h" | |
cfa7b4df | 24 | #include "AliVCluster.h" |
25 | #include "AliVTrack.h" | |
26 | #include "AliEmcalJet.h" | |
27 | #include "AliRhoParameter.h" | |
28 | #include "AliLog.h" | |
29 | #include "AliEmcalParticle.h" | |
30 | #include "AliMCEvent.h" | |
31 | #include "AliGenPythiaEventHeader.h" | |
32 | #include "AliAODMCHeader.h" | |
33 | #include "AliMCEvent.h" | |
34 | #include "AliAnalysisManager.h" | |
35 | #include "AliJetContainer.h" | |
36 | #include "AliParticleContainer.h" | |
b929303d | 37 | #include "AliPythiaInfo.h" |
364835fd | 38 | #include "TRandom3.h" |
cfa7b4df | 39 | |
40 | ||
41 | #include "AliAODEvent.h" | |
42 | ||
43 | #include "AliAnalysisTaskEmcalQGTagging.h" | |
44 | ||
e8c28d9e | 45 | using std::cout; |
46 | using std::endl; | |
47 | ||
cfa7b4df | 48 | ClassImp(AliAnalysisTaskEmcalQGTagging) |
49 | ||
50 | //________________________________________________________________________ | |
51 | AliAnalysisTaskEmcalQGTagging::AliAnalysisTaskEmcalQGTagging() : | |
52 | AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalQGTagging", kTRUE), | |
53 | fContainer(0), | |
54 | fMinFractionShared(0), | |
a68b1989 | 55 | fJetShapeType(kData), |
56 | fJetShapeSub(kNoSub), | |
364835fd | 57 | fJetSelection(kInclusive), |
cfa7b4df | 58 | fShapesVar(0), |
34b1e79e | 59 | fPtThreshold(-9999.), |
cfa7b4df | 60 | fRMatching(0.3), |
364835fd | 61 | fminpTTrig(20.), |
62 | fmaxpTTrig(50.), | |
63 | fangWindowRecoil(0.6), | |
67a0e4a8 | 64 | fh2ResponseUW(0x0), |
65 | fh2ResponseW(0x0), | |
cfa7b4df | 66 | fPhiJetCorr6(0x0), |
67 | fPhiJetCorr7(0x0), | |
68 | fEtaJetCorr6(0x0), | |
69 | fEtaJetCorr7(0x0), | |
70 | fPtJetCorr(0x0), | |
71 | fPtJet(0x0), | |
364835fd | 72 | fhpTjetpT(0x0), |
73 | fhPt(0x0), | |
74 | fhPhi(0x0), | |
cfa7b4df | 75 | fTreeObservableTagging(0) |
34b1e79e | 76 | |
cfa7b4df | 77 | { |
78 | SetMakeGeneralHistograms(kTRUE); | |
79 | } | |
80 | ||
81 | //________________________________________________________________________ | |
82 | AliAnalysisTaskEmcalQGTagging::AliAnalysisTaskEmcalQGTagging(const char *name) : | |
83 | AliAnalysisTaskEmcalJet(name, kTRUE), | |
84 | fContainer(0), | |
85 | fMinFractionShared(0), | |
a68b1989 | 86 | fJetShapeType(kData), |
87 | fJetShapeSub(kNoSub), | |
364835fd | 88 | fJetSelection(kInclusive), |
cfa7b4df | 89 | fShapesVar(0), |
34b1e79e | 90 | fPtThreshold(-9999.), |
cfa7b4df | 91 | fRMatching(0.3), |
364835fd | 92 | fminpTTrig(20.), |
93 | fmaxpTTrig(50.), | |
94 | fangWindowRecoil(0.6), | |
67a0e4a8 | 95 | fh2ResponseUW(0x0), |
96 | fh2ResponseW(0x0), | |
cfa7b4df | 97 | fPhiJetCorr6(0x0), |
98 | fPhiJetCorr7(0x0), | |
99 | fEtaJetCorr6(0x0), | |
100 | fEtaJetCorr7(0x0), | |
101 | fPtJetCorr(0x0), | |
102 | fPtJet(0x0), | |
364835fd | 103 | fhpTjetpT(0x0), |
104 | fhPt(0x0), | |
105 | fhPhi(0x0), | |
cfa7b4df | 106 | fTreeObservableTagging(0) |
34b1e79e | 107 | |
cfa7b4df | 108 | { |
109 | // Standard constructor. | |
110 | ||
111 | SetMakeGeneralHistograms(kTRUE); | |
112 | ||
113 | DefineOutput(1, TTree::Class()); | |
114 | ||
115 | } | |
116 | ||
117 | //________________________________________________________________________ | |
118 | AliAnalysisTaskEmcalQGTagging::~AliAnalysisTaskEmcalQGTagging() | |
119 | { | |
120 | // Destructor. | |
121 | } | |
122 | ||
123 | //________________________________________________________________________ | |
124 | void AliAnalysisTaskEmcalQGTagging::UserCreateOutputObjects() | |
125 | { | |
126 | // Create user output. | |
127 | ||
128 | AliAnalysisTaskEmcalJet::UserCreateOutputObjects(); | |
129 | ||
130 | Bool_t oldStatus = TH1::AddDirectoryStatus(); | |
131 | TH1::AddDirectory(kFALSE); | |
132 | ||
133 | fTreeObservableTagging = new TTree("fTreeJetShape", "fTreeJetShape"); | |
6e3203fb | 134 | Int_t nVar = 18; |
cfa7b4df | 135 | fShapesVar = new Float_t [nVar]; |
136 | TString *fShapesVarNames = new TString [nVar]; | |
137 | ||
138 | fShapesVarNames[0] = "partonCode"; | |
139 | fShapesVarNames[1] = "ptJet"; | |
140 | fShapesVarNames[2] = "ptDJet"; | |
141 | fShapesVarNames[3] = "mJet"; | |
142 | fShapesVarNames[4] = "nbOfConst"; | |
143 | fShapesVarNames[5] = "angularity"; | |
144 | fShapesVarNames[6] = "circularity"; | |
145 | fShapesVarNames[7] = "lesub"; | |
146 | fShapesVarNames[8] = "sigma2"; | |
a68b1989 | 147 | |
148 | fShapesVarNames[9] = "ptJetMatch"; | |
149 | fShapesVarNames[10] = "ptDJetMatch"; | |
150 | fShapesVarNames[11] = "mJetMatch"; | |
151 | fShapesVarNames[12] = "nbOfConstMatch"; | |
152 | fShapesVarNames[13] = "angularityMatch"; | |
153 | fShapesVarNames[14] = "circularityMatch"; | |
154 | fShapesVarNames[15] = "lesubMatch"; | |
155 | fShapesVarNames[16] = "sigma2Match"; | |
6e3203fb | 156 | fShapesVarNames[17]="weightPythia"; |
157 | ||
cfa7b4df | 158 | for(Int_t ivar=0; ivar < nVar; ivar++){ |
159 | cout<<"looping over variables"<<endl; | |
160 | fTreeObservableTagging->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/F", fShapesVarNames[ivar].Data())); | |
161 | ||
162 | //if( ivar == 4 ) fTreeObservableTagging->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/I", fShapesVarNames[ivar].Data())); | |
163 | ||
164 | } | |
67a0e4a8 | 165 | |
166 | fh2ResponseUW= new TH2F("fh2ResponseUW", "fh2ResponseUW", 100, 0, 200, 100, 0, 200); | |
167 | fOutput->Add(fh2ResponseUW); | |
364835fd | 168 | fh2ResponseW= new TH2F("fh2ResponseW", "fh2ResponseW", 100, 0, 200, 100, 0, 200); |
169 | fOutput->Add(fh2ResponseW); | |
cfa7b4df | 170 | fPhiJetCorr6= new TH2F("fPhiJetCorr6", "fPhiJetCorr6", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi()); |
171 | fOutput->Add(fPhiJetCorr6); | |
172 | fEtaJetCorr6= new TH2F("fEtaJetCorr6", "fEtaJetCorr6", 50, -1.5, 1.5, 50, -1.5, 1.5); | |
173 | fOutput->Add(fEtaJetCorr6); | |
174 | ||
175 | fPhiJetCorr7= new TH2F("fPhiJetCorr7", "fPhiJetCorr7", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi()); | |
176 | fOutput->Add(fPhiJetCorr7); | |
177 | fEtaJetCorr7= new TH2F("fEtaJetCorr7", "fEtaJetCorr7", 50, -1.5, 1.5, 50, -1.5, 1.5); | |
178 | fOutput->Add(fEtaJetCorr7); | |
179 | ||
180 | fPtJetCorr= new TH2F("fPtJetCorr", "fPtJetCorr", 100, 0, 200, 100, 0, 200); | |
181 | fOutput->Add(fPtJetCorr); | |
182 | fPtJet= new TH1F("fPtJet", "fPtJet", 100, 0, 200); | |
364835fd | 183 | fOutput->Add(fPtJet); |
184 | ||
185 | fhpTjetpT= new TH2F("fhpTjetpT", "fhpTjetpT", 100, 0, 200, 100, 0, 200); | |
186 | fOutput->Add(fhpTjetpT); | |
187 | fhPt= new TH1F("fhPt", "fhPt", 100, 0, 200); | |
188 | fOutput->Add(fhPt); | |
189 | fhPhi= new TH1F("fhPhi", "fhPhi", 100, -TMath::Pi(), TMath::Pi()); | |
190 | fOutput->Add(fhPhi); | |
191 | ||
cfa7b4df | 192 | fOutput->Add(fTreeObservableTagging); |
193 | TH1::AddDirectory(oldStatus); | |
194 | PostData(1, fOutput); // Post data for ALL output slots > 0 here. | |
195 | ||
196 | } | |
197 | ||
198 | //________________________________________________________________________ | |
199 | Bool_t AliAnalysisTaskEmcalQGTagging::Run() | |
200 | { | |
201 | // Run analysis code here, if needed. It will be executed before FillHistograms(). | |
202 | ||
203 | return kTRUE; | |
204 | } | |
205 | ||
206 | //________________________________________________________________________ | |
207 | Bool_t AliAnalysisTaskEmcalQGTagging::FillHistograms() | |
208 | { | |
209 | // Fill histograms. | |
210 | //cout<<"base container"<<endl; | |
211 | AliEmcalJet* jet1 = NULL; | |
212 | AliJetContainer *jetCont = GetJetContainer(0); | |
b929303d | 213 | Float_t kWeight=1; |
2daefd7c | 214 | if(fCent>10) return 0; |
364835fd | 215 | |
216 | AliAODTrack *triggerHadron = 0x0; | |
217 | ||
218 | if (fJetSelection == kRecoil) { | |
219 | //Printf("Recoil jets!!!, fminpTTrig = %f, fmaxpTTrig = %f", fminpTTrig, fmaxpTTrig); | |
220 | Int_t triggerHadronLabel = SelectTrigger(fminpTTrig, fmaxpTTrig); | |
221 | if (triggerHadronLabel==-9999) { | |
222 | Printf ("Trigger Hadron not found, return"); | |
223 | return 0; | |
224 | } | |
225 | AliParticleContainer *partContAn = GetParticleContainer(0); | |
226 | TClonesArray *trackArrayAn = partContAn->GetArray(); | |
227 | triggerHadron = static_cast<AliAODTrack*>(trackArrayAn->At(triggerHadronLabel)); | |
228 | if (!triggerHadron) { | |
229 | Printf("No Trigger hadron with the found label!!"); | |
230 | return 0; | |
231 | } | |
232 | } | |
233 | ||
cfa7b4df | 234 | if(jetCont) { |
235 | jetCont->ResetCurrentID(); | |
236 | while((jet1 = jetCont->GetNextAcceptJet())) { | |
237 | if (!jet1) continue; | |
a68b1989 | 238 | AliEmcalJet* jet2 = 0x0; |
cfa7b4df | 239 | fPtJet->Fill(jet1->Pt()); |
5f1a8018 | 240 | AliEmcalJet *jetUS = NULL; |
241 | Int_t ifound=0; | |
242 | Int_t ilab=-1; | |
364835fd | 243 | |
244 | if (!(fJetShapeType == kData)) { | |
245 | AliPythiaInfo *partonsInfo = 0x0; | |
246 | ||
247 | if((fJetShapeType == kTrueDet) || (fJetShapeType == kDetEmb)){ | |
248 | AliJetContainer *jetContTrue = GetJetContainer(1); | |
5f1a8018 | 249 | AliJetContainer *jetContUS = GetJetContainer(2); |
250 | if(fJetShapeSub==kConstSub){ | |
364835fd | 251 | for(Int_t i = 0; i<jetContUS->GetNJets(); i++) { |
252 | jetUS = jetContUS->GetJet(i); | |
253 | if(jetUS->GetLabel()==jet1->GetLabel()) { | |
254 | ifound++; | |
255 | if(ifound==1) ilab = i; | |
256 | } | |
257 | } | |
258 | if(ilab==-1) continue; | |
259 | jetUS=jetContUS->GetJet(ilab); | |
260 | ||
261 | jet2=jetUS->ClosestJet(); | |
262 | } | |
263 | if(!(fJetShapeSub==kConstSub)) jet2 = jet1->ClosestJet(); | |
264 | if (!jet2) { | |
265 | Printf("jet2 not exists, returning"); | |
266 | continue; | |
267 | } | |
268 | ||
269 | fh2ResponseUW->Fill(jet1->Pt(),jet2->Pt()); | |
270 | ||
5f1a8018 | 271 | Double_t fraction=0; |
364835fd | 272 | if(!(fJetShapeSub==kConstSub)) fraction = jetCont->GetFractionSharedPt(jet1); |
5f1a8018 | 273 | if(fJetShapeSub==kConstSub) fraction = jetContUS->GetFractionSharedPt(jetUS); |
2daefd7c | 274 | cout<<"hey a jet"<<fraction<<" "<<jet1->Pt()<<" "<<jet2->Pt()<<" "<<fCent<<endl; |
275 | ||
364835fd | 276 | if(fraction<fMinFractionShared) continue; |
277 | partonsInfo = (AliPythiaInfo*) jetContTrue->GetPythiaInfo(); | |
278 | if(!partonsInfo) return 0; | |
b929303d | 279 | |
364835fd | 280 | } |
281 | else { | |
282 | partonsInfo = (AliPythiaInfo*) jetCont->GetPythiaInfo(); | |
283 | jet2=jet1; | |
be58fab7 | 284 | if(!partonsInfo) return 0; |
364835fd | 285 | } |
286 | ||
287 | Double_t jp1=(jet2->Phi())-(partonsInfo->GetPartonPhi6()); | |
288 | Double_t detap1=(jet2->Eta())-(partonsInfo->GetPartonEta6()); | |
289 | kWeight=partonsInfo->GetPythiaEventWeight(); | |
67a0e4a8 | 290 | fh2ResponseW->Fill(jet1->Pt(),jet2->Pt(),kWeight); |
364835fd | 291 | if (jp1< -1*TMath::Pi()) jp1 = (-2*TMath::Pi())-jp1; |
292 | else if (jp1 > TMath::Pi()) jp1 = (2*TMath::Pi())-jp1; | |
293 | Float_t dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1); | |
294 | fEtaJetCorr6->Fill(jet2->Eta(), partonsInfo->GetPartonEta6()); | |
295 | fPhiJetCorr6->Fill(jet2->Phi(), partonsInfo->GetPartonPhi6()); | |
296 | if(dRp1 < fRMatching) { | |
297 | fShapesVar[0] = partonsInfo->GetPartonFlag6(); | |
298 | fPtJetCorr ->Fill(partonsInfo->GetPartonPt6(), jet2->Pt()); | |
299 | } | |
300 | else { | |
301 | jp1=(jet2->Phi())-(partonsInfo->GetPartonPhi7()); | |
302 | detap1=(jet2->Eta())-(partonsInfo->GetPartonEta7()); | |
303 | if (jp1< -1*TMath::Pi()) jp1= (-2*TMath::Pi())-jp1; | |
304 | else if (jp1 > TMath::Pi()) jp1 = (2*TMath::Pi())-jp1; | |
305 | dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1); | |
306 | fEtaJetCorr7->Fill(jet2->Eta(), partonsInfo->GetPartonEta7()); | |
307 | fPhiJetCorr7->Fill(jet2->Phi(), partonsInfo->GetPartonPhi7()); | |
308 | if(dRp1 < fRMatching) { | |
309 | fShapesVar[0] = partonsInfo->GetPartonFlag7(); | |
310 | fPtJetCorr->Fill(partonsInfo->GetPartonPt7(), jet2->Pt()); | |
311 | } | |
312 | else continue; | |
313 | } | |
cfa7b4df | 314 | } |
315 | else | |
364835fd | 316 | fShapesVar[0] = 0.; |
317 | ||
318 | Double_t ptSubtracted = 0; | |
319 | ||
320 | Float_t dphiRecoil = 0.; | |
321 | if (fJetSelection == kRecoil){ | |
322 | dphiRecoil = RelativePhi(triggerHadron->Phi(), jet1->Phi()); | |
323 | if (TMath::Abs(dphiRecoil) < TMath::Pi() - fangWindowRecoil) { | |
324 | // Printf("Recoil jets back to back not found! continuing"); | |
325 | continue; | |
326 | } | |
327 | //Printf("Recoil jets back to back found! filling histos!"); | |
328 | fhpTjetpT->Fill(triggerHadron->Pt(), jet1->Pt()); | |
329 | //Printf("triggerHadron =%f, jet1 = %f", triggerHadron->Pt(), jet1->Pt()); | |
330 | fhPt->Fill(triggerHadron->Pt()); | |
331 | fhPhi->Fill(RelativePhi(triggerHadron->Phi(), jet1->Phi())); | |
332 | } | |
333 | ||
a68b1989 | 334 | if (((fJetShapeType == kData) || (fJetShapeType == kDetEmb)) && (fJetShapeSub == kConstSub)) |
364835fd | 335 | ptSubtracted = jet1->Pt(); |
a68b1989 | 336 | else ptSubtracted = jet1->Pt() - GetRhoVal(0)*jet1->Area(); |
364835fd | 337 | |
338 | if ((fJetShapeType == kData || fJetShapeType== kDetEmb)) | |
339 | if (ptSubtracted < fPtThreshold) continue; | |
a68b1989 | 340 | |
34b1e79e | 341 | fShapesVar[1] = ptSubtracted; |
a68b1989 | 342 | fShapesVar[2] = GetJetpTD(jet1,0); |
e5aca8ab | 343 | fShapesVar[3] = GetJetMass(jet1,0); |
344 | fShapesVar[4] = 1.*GetJetNumberOfConstituents(jet1,0); | |
a68b1989 | 345 | fShapesVar[5] = GetJetAngularity(jet1,0); |
346 | fShapesVar[6] = GetJetCircularity(jet1,0); | |
347 | fShapesVar[7] = GetJetLeSub(jet1,0); | |
348 | fShapesVar[8] = GetSigma2(jet1,0); | |
349 | ||
350 | Float_t ptMatch=0., ptDMatch=0., massMatch=0., constMatch=0.,angulMatch=0.,circMatch=0., lesubMatch=0., sigma2Match=0.; | |
351 | Int_t kMatched = 0; | |
352 | if (fJetShapeType == kTrueDet || fJetShapeType == kDetEmb) { | |
364835fd | 353 | kMatched = 1; |
354 | ptMatch=jet2->Pt(); | |
355 | ptDMatch=GetJetpTD(jet2, kMatched); | |
356 | massMatch=GetJetMass(jet2,kMatched); | |
357 | constMatch=1.*GetJetNumberOfConstituents(jet2,kMatched); | |
358 | angulMatch=GetJetAngularity(jet2, kMatched); | |
359 | circMatch=GetJetCircularity(jet2, kMatched); | |
360 | lesubMatch=GetJetLeSub(jet2, kMatched); | |
361 | sigma2Match = GetSigma2(jet2, kMatched); | |
6e3203fb | 362 | |
a68b1989 | 363 | } |
6e3203fb | 364 | |
364835fd | 365 | if (fJetShapeType == kTrue || fJetShapeType == kData) { |
366 | kMatched = 0; | |
367 | ptMatch=0.; | |
368 | ptDMatch=0.; | |
369 | massMatch=0.; | |
370 | constMatch=0.; | |
371 | angulMatch=0.; | |
372 | circMatch=0.; | |
373 | lesubMatch=0.; | |
374 | sigma2Match =0.; | |
375 | ||
a68b1989 | 376 | } |
377 | ||
378 | fShapesVar[9] = ptMatch; | |
379 | fShapesVar[10] = ptDMatch; | |
380 | fShapesVar[11] = massMatch; | |
381 | fShapesVar[12] = constMatch; | |
382 | fShapesVar[13] = angulMatch; | |
383 | fShapesVar[14] = circMatch; | |
384 | fShapesVar[15] = lesubMatch; | |
385 | fShapesVar[16] = sigma2Match; | |
6e3203fb | 386 | fShapesVar[17] = kWeight; |
cfa7b4df | 387 | fTreeObservableTagging->Fill(); |
364835fd | 388 | |
cfa7b4df | 389 | } |
390 | ||
364835fd | 391 | } |
cfa7b4df | 392 | |
393 | return kTRUE; | |
394 | } | |
395 | ||
396 | //________________________________________________________________________ | |
e5aca8ab | 397 | Float_t AliAnalysisTaskEmcalQGTagging::GetJetMass(AliEmcalJet *jet,Int_t jetContNb=0) { |
cfa7b4df | 398 | //calc subtracted jet mass |
e5aca8ab | 399 | if((fJetShapeSub==kDerivSub)&&(jetContNb==0)) |
cfa7b4df | 400 | return jet->GetSecondOrderSubtracted(); |
401 | else | |
402 | return jet->M(); | |
403 | } | |
404 | ||
405 | //________________________________________________________________________ | |
a68b1989 | 406 | Float_t AliAnalysisTaskEmcalQGTagging::Angularity(AliEmcalJet *jet, Int_t jetContNb = 0){ |
cfa7b4df | 407 | |
a68b1989 | 408 | AliJetContainer *jetCont = GetJetContainer(jetContNb); |
cfa7b4df | 409 | if (!jet->GetNumberOfTracks()) |
410 | return 0; | |
411 | Double_t den=0.; | |
412 | Double_t num = 0.; | |
413 | AliVParticle *vp1 = 0x0; | |
414 | for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) { | |
415 | vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray())); | |
416 | Double_t dphi = vp1->Phi()-jet->Phi(); | |
417 | if(dphi<-1.*TMath::Pi()) dphi+=TMath::TwoPi(); | |
418 | if(dphi>TMath::Pi()) dphi-=TMath::TwoPi(); | |
419 | Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi; | |
420 | Double_t dr = TMath::Sqrt(dr2); | |
421 | num=num+vp1->Pt()*dr; | |
422 | den=den+vp1->Pt(); | |
423 | } | |
424 | return num/den; | |
425 | } | |
426 | ||
427 | //________________________________________________________________________ | |
a68b1989 | 428 | Float_t AliAnalysisTaskEmcalQGTagging::GetJetAngularity(AliEmcalJet *jet, Int_t jetContNb = 0) { |
cfa7b4df | 429 | |
e5aca8ab | 430 | if((fJetShapeSub==kDerivSub) && (jetContNb==0)) |
cfa7b4df | 431 | return jet->GetSecondOrderSubtractedAngularity(); |
432 | else | |
a68b1989 | 433 | return Angularity(jet, jetContNb); |
cfa7b4df | 434 | |
435 | } | |
436 | ||
437 | ||
438 | //________________________________________________________________________ | |
a68b1989 | 439 | Float_t AliAnalysisTaskEmcalQGTagging::PTD(AliEmcalJet *jet, Int_t jetContNb = 0){ |
cfa7b4df | 440 | |
a68b1989 | 441 | AliJetContainer *jetCont = GetJetContainer(jetContNb); |
cfa7b4df | 442 | if (!jet->GetNumberOfTracks()) |
443 | return 0; | |
444 | Double_t den=0.; | |
445 | Double_t num = 0.; | |
446 | AliVParticle *vp1 = 0x0; | |
447 | for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) { | |
448 | vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray())); | |
449 | num=num+vp1->Pt()*vp1->Pt(); | |
450 | den=den+vp1->Pt(); | |
451 | } | |
452 | return TMath::Sqrt(num)/den; | |
453 | } | |
454 | ||
455 | //________________________________________________________________________ | |
a68b1989 | 456 | Float_t AliAnalysisTaskEmcalQGTagging::GetJetpTD(AliEmcalJet *jet, Int_t jetContNb = 0) { |
cfa7b4df | 457 | //calc subtracted jet mass |
e5aca8ab | 458 | if((fJetShapeSub==kDerivSub)&&(jetContNb==0)) |
cfa7b4df | 459 | return jet->GetSecondOrderSubtractedpTD(); |
460 | else | |
a68b1989 | 461 | return PTD(jet, jetContNb); |
cfa7b4df | 462 | |
463 | } | |
464 | ||
465 | //_____________________________________________________________________________ | |
a68b1989 | 466 | Float_t AliAnalysisTaskEmcalQGTagging::Circularity(AliEmcalJet *jet, Int_t jetContNb = 0){ |
cfa7b4df | 467 | |
a68b1989 | 468 | AliJetContainer *jetCont = GetJetContainer(jetContNb); |
cfa7b4df | 469 | if (!jet->GetNumberOfTracks()) |
470 | return 0; | |
471 | Double_t mxx = 0.; | |
472 | Double_t myy = 0.; | |
473 | Double_t mxy = 0.; | |
474 | int nc = 0; | |
475 | Double_t sump2 = 0.; | |
476 | Double_t pxjet=jet->Px(); | |
477 | Double_t pyjet=jet->Py(); | |
478 | Double_t pzjet=jet->Pz(); | |
479 | ||
480 | ||
481 | //2 general normalized vectors perpendicular to the jet | |
482 | TVector3 ppJ1(pxjet, pyjet, pzjet); | |
483 | TVector3 ppJ3(- pxjet* pzjet, - pyjet * pzjet, pxjet * pxjet + pyjet * pyjet); | |
484 | ppJ3.SetMag(1.); | |
485 | TVector3 ppJ2(-pyjet, pxjet, 0); | |
486 | ppJ2.SetMag(1.); | |
487 | AliVParticle *vp1 = 0x0; | |
488 | for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) { | |
489 | vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray())); | |
490 | ||
491 | ||
492 | TVector3 pp(vp1->Px(), vp1->Py(), vp1->Pz()); | |
493 | ||
494 | //local frame | |
495 | TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1; | |
496 | TVector3 pPerp = pp - pLong; | |
497 | //projection onto the two perpendicular vectors defined above | |
498 | ||
499 | Float_t ppjX = pPerp.Dot(ppJ2); | |
500 | Float_t ppjY = pPerp.Dot(ppJ3); | |
501 | Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY); | |
502 | if(ppjT<=0) return 0; | |
503 | ||
504 | mxx += (ppjX * ppjX / ppjT); | |
505 | myy += (ppjY * ppjY / ppjT); | |
506 | mxy += (ppjX * ppjY / ppjT); | |
507 | nc++; | |
508 | sump2 += ppjT;} | |
509 | ||
510 | if(nc<2) return 0; | |
511 | if(sump2==0) return 0; | |
512 | // Sphericity Matrix | |
513 | Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2}; | |
514 | TMatrixDSym m0(2,ele); | |
515 | ||
516 | // Find eigenvectors | |
517 | TMatrixDSymEigen m(m0); | |
518 | TVectorD eval(2); | |
519 | TMatrixD evecm = m.GetEigenVectors(); | |
520 | eval = m.GetEigenValues(); | |
521 | // Largest eigenvector | |
522 | int jev = 0; | |
523 | // cout<<eval[0]<<" "<<eval[1]<<endl; | |
524 | if (eval[0] < eval[1]) jev = 1; | |
525 | TVectorD evec0(2); | |
526 | // Principle axis | |
527 | evec0 = TMatrixDColumn(evecm, jev); | |
528 | Double_t compx=evec0[0]; | |
529 | Double_t compy=evec0[1]; | |
530 | TVector2 evec(compx, compy); | |
531 | Double_t circ=0; | |
532 | if(jev==1) circ=2*eval[0]; | |
533 | if(jev==0) circ=2*eval[1]; | |
534 | ||
535 | return circ; | |
536 | ||
537 | ||
538 | ||
539 | } | |
540 | ||
541 | ||
542 | ||
543 | ||
544 | //________________________________________________________________________ | |
a68b1989 | 545 | Float_t AliAnalysisTaskEmcalQGTagging::GetJetCircularity(AliEmcalJet *jet, Int_t jetContNb =0 ) { |
cfa7b4df | 546 | //calc subtracted jet mass |
547 | ||
e5aca8ab | 548 | if((fJetShapeSub==kDerivSub)&&(jetContNb==0)) |
cfa7b4df | 549 | return jet->GetSecondOrderSubtractedCircularity(); |
550 | else | |
a68b1989 | 551 | return Circularity(jet, jetContNb); |
cfa7b4df | 552 | |
553 | } | |
554 | ||
555 | //________________________________________________________________________ | |
a68b1989 | 556 | Float_t AliAnalysisTaskEmcalQGTagging::LeSub(AliEmcalJet *jet, Int_t jetContNb =0 ){ |
cfa7b4df | 557 | |
a68b1989 | 558 | AliJetContainer *jetCont = GetJetContainer(jetContNb); |
cfa7b4df | 559 | if (!jet->GetNumberOfTracks()) |
560 | return 0; | |
561 | Double_t den=0.; | |
562 | Double_t num = 0.; | |
563 | AliVParticle *vp1 = 0x0; | |
564 | AliVParticle *vp2 = 0x0; | |
565 | std::vector<int> ordindex; | |
566 | ordindex=jet->SortConstituentsPt(jetCont->GetParticleContainer()->GetArray()); | |
567 | ||
568 | vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[0], jetCont->GetParticleContainer()->GetArray())); | |
569 | vp2 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[1], jetCont->GetParticleContainer()->GetArray())); | |
570 | ||
571 | num=vp1->Pt(); | |
572 | den=vp2->Pt(); | |
573 | ||
574 | return num-den; | |
575 | } | |
576 | ||
577 | //________________________________________________________________________ | |
a68b1989 | 578 | Float_t AliAnalysisTaskEmcalQGTagging::GetJetLeSub(AliEmcalJet *jet, Int_t jetContNb =0) { |
cfa7b4df | 579 | //calc subtracted jet mass |
580 | ||
e5aca8ab | 581 | if((fJetShapeSub==kDerivSub)&&(jetContNb==0)) |
cfa7b4df | 582 | return jet->GetSecondOrderSubtractedLeSub(); |
583 | else | |
a68b1989 | 584 | return LeSub(jet, jetContNb); |
cfa7b4df | 585 | |
586 | } | |
587 | ||
588 | //________________________________________________________________________ | |
e5aca8ab | 589 | Float_t AliAnalysisTaskEmcalQGTagging::GetJetNumberOfConstituents(AliEmcalJet *jet,Int_t jetContNb=0) { |
cfa7b4df | 590 | //calc subtracted jet mass |
591 | ||
e5aca8ab | 592 | if((fJetShapeSub==kDerivSub)&&(jetContNb==0)) |
cfa7b4df | 593 | return jet->GetSecondOrderSubtractedConstituent(); |
594 | else | |
595 | return jet->GetNumberOfTracks(); | |
596 | ||
597 | } | |
598 | ||
599 | ||
600 | //______________________________________________________________________________ | |
a68b1989 | 601 | Float_t AliAnalysisTaskEmcalQGTagging::Sigma2(AliEmcalJet *jet, Int_t jetContNb=0){ |
cfa7b4df | 602 | |
a68b1989 | 603 | AliJetContainer *jetCont = GetJetContainer(jetContNb); |
cfa7b4df | 604 | if (!jet->GetNumberOfTracks()) |
605 | return 0; | |
606 | Double_t mxx = 0.; | |
607 | Double_t myy = 0.; | |
608 | Double_t mxy = 0.; | |
609 | int nc = 0; | |
610 | Double_t sump2 = 0.; | |
611 | ||
612 | AliVParticle *vp1 = 0x0; | |
613 | for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) { | |
614 | vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray())); | |
615 | Double_t ppt=vp1->Pt(); | |
616 | Double_t dphi = vp1->Phi()-jet->Phi(); | |
617 | if(dphi<-1.*TMath::Pi()) dphi+=TMath::TwoPi(); | |
618 | if(dphi>TMath::Pi()) dphi-=TMath::TwoPi(); | |
619 | Double_t deta = vp1->Eta()-jet->Eta(); | |
620 | mxx += ppt*ppt*deta*deta; | |
621 | myy += ppt*ppt*dphi*dphi; | |
622 | mxy -= ppt*ppt*deta*dphi; | |
623 | nc++; | |
624 | sump2 += ppt*ppt; | |
625 | ||
626 | } | |
627 | if(nc<2) return 0; | |
628 | if(sump2==0) return 0; | |
629 | // Sphericity Matrix | |
630 | Double_t ele[4] = {mxx , mxy , mxy , myy }; | |
631 | TMatrixDSym m0(2,ele); | |
632 | ||
633 | // Find eigenvectors | |
634 | TMatrixDSymEigen m(m0); | |
635 | TVectorD eval(2); | |
636 | TMatrixD evecm = m.GetEigenVectors(); | |
637 | eval = m.GetEigenValues(); | |
638 | // Largest eigenvector | |
639 | int jev = 0; | |
640 | // cout<<eval[0]<<" "<<eval[1]<<endl; | |
641 | if (eval[0] < eval[1]) jev = 1; | |
642 | TVectorD evec0(2); | |
643 | // Principle axis | |
644 | evec0 = TMatrixDColumn(evecm, jev); | |
645 | Double_t compx=evec0[0]; | |
646 | Double_t compy=evec0[1]; | |
647 | TVector2 evec(compx, compy); | |
648 | Double_t sig=0; | |
649 | if(jev==1) sig=TMath::Sqrt(TMath::Abs(eval[0])/sump2); | |
650 | if(jev==0) sig=TMath::Sqrt(TMath::Abs(eval[1])/sump2); | |
651 | ||
652 | return sig; | |
653 | ||
654 | } | |
655 | ||
656 | //________________________________________________________________________ | |
a68b1989 | 657 | Float_t AliAnalysisTaskEmcalQGTagging::GetSigma2(AliEmcalJet *jet, Int_t jetContNb=0) { |
cfa7b4df | 658 | //calc subtracted jet mass |
659 | ||
e5aca8ab | 660 | if((fJetShapeSub==kDerivSub)&&(jetContNb==0)) |
cfa7b4df | 661 | return jet->GetSecondOrderSubtractedSigma2(); |
662 | else | |
a68b1989 | 663 | return Sigma2(jet, jetContNb); |
cfa7b4df | 664 | |
665 | } | |
666 | ||
364835fd | 667 | //________________________________________________________________________ |
668 | Int_t AliAnalysisTaskEmcalQGTagging::SelectTrigger(Float_t minpT, Float_t maxpT){ | |
669 | ||
670 | AliParticleContainer *partCont = GetParticleContainer(0); | |
671 | TClonesArray *tracksArray = partCont->GetArray(); | |
672 | ||
673 | if(!partCont || !tracksArray) return 0; | |
674 | AliAODTrack *track = 0x0; | |
675 | ||
676 | TList *trackList = new TList(); | |
677 | Int_t triggers[100]; | |
678 | for (Int_t iTrigger=0; iTrigger<100; iTrigger++) triggers[iTrigger] = 0; | |
679 | Int_t iTT = 0; | |
680 | ||
681 | for(Int_t iTrack=0; iTrack <= tracksArray->GetEntriesFast(); iTrack++){ | |
682 | track = static_cast<AliAODTrack*>(tracksArray->At(iTrack)); | |
683 | if (!track) continue; | |
684 | ||
685 | if(TMath::Abs(track->Eta())>0.9) continue; | |
686 | if (track->Pt()<0.15) continue; | |
687 | if (!(track->TestFilterBit(768))) continue; | |
688 | ||
689 | if ((track->Pt() >= minpT) && (track->Pt()< maxpT)) { | |
690 | trackList->Add(track); | |
691 | triggers[iTT] = iTrack; | |
692 | iTT++; | |
693 | } | |
694 | } | |
695 | ||
696 | if (iTT == 0) return -9999; | |
697 | Int_t nbRn = 0, index = 0 ; | |
698 | TRandom3* random = new TRandom3(0); | |
699 | nbRn = random->Integer(iTT); | |
700 | ||
701 | index = triggers[nbRn]; | |
702 | //Printf("iTT= %d, nbRn = %d, Index = %d",iTT, nbRn, index ); | |
703 | return index; | |
704 | ||
705 | } | |
706 | ||
707 | //__________________________________________________________________________________ | |
708 | Double_t AliAnalysisTaskEmcalQGTagging::RelativePhi(Double_t mphi,Double_t vphi){ | |
709 | ||
710 | if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi()); | |
711 | else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi()); | |
712 | if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi()); | |
713 | else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi()); | |
714 | double dphi = mphi-vphi; | |
715 | if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi()); | |
716 | else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi()); | |
717 | return dphi;//dphi in [-Pi, Pi] | |
718 | } | |
719 | ||
720 | ||
cfa7b4df | 721 | //________________________________________________________________________ |
722 | Bool_t AliAnalysisTaskEmcalQGTagging::RetrieveEventObjects() { | |
723 | // | |
724 | // retrieve event objects | |
725 | // | |
726 | if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects()) | |
727 | return kFALSE; | |
728 | ||
729 | return kTRUE; | |
730 | } | |
731 | ||
732 | //_______________________________________________________________________ | |
733 | void AliAnalysisTaskEmcalQGTagging::Terminate(Option_t *) | |
734 | { | |
735 | // Called once at the end of the analysis. | |
736 | ||
737 | // fTreeObservableTagging = dynamic_cast<TTree*>(GetOutputData(1)); | |
738 | // if (!fTreeObservableTagging){ | |
739 | // Printf("ERROR: fTreeObservableTagging not available"); | |
740 | // return; | |
741 | // } | |
742 | ||
743 | } | |
744 |