AliPhysics  a9863a5 (a9863a5)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskJetShapeConst.cxx
Go to the documentation of this file.
1 #include <TH1F.h>
2 #include <TH2F.h>
3 #include <TH3F.h>
4 #include <THnSparse.h>
5 #include <TLorentzVector.h>
6 #include <TTree.h>
7 
8 #include <AliParticleContainer.h>
9 #include <AliJetContainer.h>
10 #include <AliRhoParameter.h>
11 #include <AliEmcalJet.h>
12 #include <AliVParticle.h>
13 
15 
17 
18 //________________________________________________________________________
21  fhptjetSMinusSingleTrack(0x0),
22  fhJet1vsJetTag(0x0),
23  fhNconstit(0x0),
24  fhAreaJet(0x0)
25 {
26  // Default constructor.
27 
28 
29 }
30 
31 //________________________________________________________________________
34  fhptjetSMinusSingleTrack(0x0),
35  fhJet1vsJetTag(0x0),
36  fhNconstit(0x0),
37  fhAreaJet(0x0)
38 {
39  // Standard constructor.
40 
41 }
42 
43 //________________________________________________________________________
45 {
46  // Destructor.
47 }
48 
49 //________________________________________________________________________
51 {
52  // Create user output.
53 
55  Bool_t oldStatus = TH1::AddDirectoryStatus();
56  TH1::AddDirectory(kFALSE);
57 
58  fhptjetSMinusSingleTrack = new TH1F("fhptjetSMinusSingleTrack", "Subtraction of single track #it{p}_{T}; #it{p}_{T, jet} - #it{p}_{T, Emb Track};Entries", 500,-10.,110.);
60 
61  fhJet1vsJetTag = new TH2F("fhJet1vsJetTag", "Number of jets vs tagged jets; #it{N}_{jet,Tot}; #it{N}_{jet,Tag}", 30, 1., 30., 30, 1., 30.);
62  fOutput->Add(fhJet1vsJetTag);
63 
64  fhNconstit = new TH1F("fhNconstit", "Number of constituents (matched jets); #it{N}_{constituents}", 21, 0., 20.);
65  fOutput->Add(fhNconstit);
66 
67  fhAreaJet = new TH1F("fhAreaJet", "Area (matched jets); Area", 400., 0., 4);
68  fOutput->Add(fhAreaJet);
69 
70  TH1::AddDirectory(oldStatus);
71 
72  PostData(1, fOutput);
73 }
74 
75 //________________________________________________________________________
76 //Bool_t AliAnalysisTaskJetShapeConst::Run()
77 //{
78 // // Run analysis code here, if needed. It will be executed before FillHistograms().
79 //
80 // return kTRUE;
81 //}
82 
83 //________________________________________________________________________
85 {
86  // Fill histograms.
87 
88  AliEmcalJet* jet1 = NULL; //AA jet
89  AliEmcalJet *jet2 = NULL; //Embedded Pythia jet
90  AliEmcalJet *jetR = NULL; //true jet for response matrix
91  AliEmcalJet *jetS = NULL; //subtracted jet
92  AliEmcalJet *jetO = NULL; //hard-ish jet to avoid overlap of single track with
93  AliVParticle *vpe = NULL; //embedded particle
94 
98 
99  AliParticleContainer *trackCont = GetParticleContainer(0); //track used for the jet finder jet1
100  for(Int_t itr=0; itr < trackCont->GetNAcceptedParticles(); itr++){
101  fhpTTracksCont->Fill(trackCont->GetParticle(itr)->Pt());
102 
103  }
104  //Get leading jet in Pb-Pb event without embedded objects
106  AliEmcalJet *jetL = NULL;
107  if(jetContNoEmb) jetL = jetContNoEmb->GetLeadingJet("rho");
108 
109  if(fOverlap && !jetContO){
110  Printf("Jet Container %d not found, return", fContainerOverlap);
111  return kFALSE;
112  }
113 
114  //rho
115  AliRhoParameter* rhoParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(jetContS->GetRhoName()));
116  fRho = 0;
117  if (!rhoParam) {
118  AliError(Form("%s: Could not retrieve rho %s (some histograms will be filled with zero)!", GetName(), jetContS->GetRhoName().Data()));
119 
120  } else fRho = rhoParam->GetVal();
121  //rhom
122  AliRhoParameter* rhomParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(jetContS->GetRhoMassName()));
123  fRhoM = 0;
124  if (!rhomParam) {
125  AliError(Form("%s: Could not retrieve rho_m %s (some histograms will be filled with zero)!", GetName(), jetContS->GetRhoMassName().Data()));
126 
127  } else fRhoM = rhomParam->GetVal();
128 
129  Int_t njet1 = 0, ntagjet2 = 0;
130 
131  AliDebug(11,Form("NJets Incl: %d Csub: %d",jetCont->GetNJets(),jetContS->GetNJets()));
132  if(jetCont) {
133  jetCont->ResetCurrentID();
134  while((jet1 = jetCont->GetNextAcceptJet())) {
135  njet1++;
136  jet2 = NULL;
137  jetS = NULL;
138  if(jet1->GetTagStatus()<1 || !jet1->GetTaggedJet())
139  continue;
140 
141  ntagjet2++;
142  //print constituents of different jet containers
143  //jet1
144 
145  for(Int_t i=0; i<jet1->GetNumberOfTracks(); i++) {
146  AliVParticle *vp = static_cast<AliVParticle*>(jet1->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
147  // if (vp->TestBits(TObject::kBitMask) != (Int_t)(TObject::kBitMask) ) continue;
148  //Int_t lab = TMath::Abs(vp->GetLabel());
149  if(vp) fhpTTracksJet1 -> Fill(vp->Pt());
150  }
151 
152  //Get constituent subtacted version of jet
153 
154  Int_t ifound = 0;
155  Int_t ilab = -1;
156  for(Int_t i = 0; i<jetContS->GetNJets(); i++) {
157  //if(ifound==1) continue;
158  jetS = jetContS->GetJet(i);
159  if(!jetS) continue;
160  if(jetS->GetLabel()==jet1->GetLabel()) {
161  ifound++;
162  if(ifound==1) ilab = i;
163  }
164  }
165  if(ifound>1) AliDebug(2,Form("Found %d partners",ifound));
166  if(ifound==0) jetS = 0x0;
167  else jetS = jetContS->GetJet(ilab);
168  if(!jetS) continue;
169 
170  Double_t mjetS = jetS->M();
171  Double_t mUnsubjet1 = jet1->M();
172  Double_t ptjet1 = jet1->Pt()-jetCont->GetRhoVal()*jet1->Area();
173  Double_t ptjetS = jetS->Pt();
174  Double_t ptUnsubjet1 = jet1->Pt();
175  Double_t var = mjetS;
176  //Double_t ptjetSMinusEmbTrpt = ptjetS;
177  Double_t ptjetSMinusEmbTrpt = ptUnsubjet1;
178  if(fJetMassVarType==kRatMPt) {
179  if(ptjetS>0. || ptjetS<0.) var = mjetS/ptjetS;
180  else var = -999.;
181  }
182 
183  //Fill histograms for all AA jets
184  fh2MSubPtRawAll[fCentBin]->Fill(var,ptjetS);
185 
186  Double_t fraction = 0.;
187  fMatch = 0;
188  fJet2Vec->SetPtEtaPhiM(0.,0.,0.,0.);
189  //Printf("EMB task: pT jet %d = %f", njet1, jet1->Pt());
190  if(fSingleTrackEmb) {
191  vpe = GetEmbeddedConstituent(jet1);
192  if(vpe){
193  Bool_t reject = kFALSE;
194  ptjetSMinusEmbTrpt -= vpe->Pt();
195  fhptjetSMinusSingleTrack->Fill(ptjetSMinusEmbTrpt);
196  //Printf("EMB task: pT jet matched = %f -> %f", jet1->Pt(), ptjetSMinusEmbTrpt);
197  if(fOverlap){
198  Int_t Njets = jetContO->GetNAcceptedJets();
199  fhNJetsSelEv->Fill(Njets);
200  jetContO->ResetCurrentID();
201  while((jetO = jetContO->GetNextAcceptJet())){
202  //print constituents of different jet containers
203  //jetO
204  //Printf("N particle %d",jetO->GetNumberOfTracks());
205  for(Int_t i=0; i<jetO->GetNumberOfTracks(); i++) {
206  AliVParticle* vp = static_cast<AliVParticle*>(jetO->TrackAt(i, jetContO->GetParticleContainer()->GetArray()));
207  // if (vp->TestBits(TObject::kBitMask) != (Int_t)(TObject::kBitMask) ) continue;
208  //Int_t lab = TMath::Abs(vp->GetLabel());
209  if(vp) fhpTTracksJetO -> Fill(vp->Pt());
210  }
211  Double_t deltaR = jetO->DeltaR(vpe);
212  fhRjetTrvspTj->Fill(deltaR, jetO->Pt());
213  fhJetEtaPhiOvl->Fill(jetO->Eta(), jetO->Phi());
214  if(deltaR < fRadius) {
215  reject = kTRUE;
216  break;
217  }
218  }
219  }
220  if(!reject) {
221  fJet2Vec->SetPxPyPzE(vpe->Px(),vpe->Py(),vpe->Pz(),vpe->E());
222  fMatch = 1;
223  }
224  }
225  } else {
226  jet2 = jet1->ClosestJet();
227  fraction = jetCont->GetFractionSharedPt(jet1);
228  fMatch = 1;
229  if(fMinFractionShared>0.) {
230  if(fraction>fMinFractionShared) {
231  fJet2Vec->SetPxPyPzE(jet2->Px(),jet2->Py(),jet2->Pz(),jet2->E());
232  fMatch = 1;
233 
234  //choose jet type for true axis of response matrix
236  jetR = jet2;
237  else if(fResponseReference==kPart)
238  jetR = jet2->GetTaggedJet();
239  } else
240  fMatch = 0;
241  }
242  }
243 
244  // if(fMatch==1 && jet2->GetTagStatus()>0) jet2T = jet2->GetTaggedJet();
245 
246  //Fill histograms for matched jets
247  fh2MSubMatch[fCentBin]->Fill(var,fMatch);
248 
249  if(fMatch==1) {
250  fhJetSubMatchEtaPhiPt->Fill(jetS->Eta(), jetS->Phi(), ptjetS);
251  Double_t drToLJ = -1.;
252  if(jetL) drToLJ = jet1->DeltaR(jetL);
253  if(fSingleTrackEmb && vpe)
254  drToLJ = jet1->DeltaR(vpe);
255  fh3MSubPtRawDRMatch[fCentBin]->Fill(var,ptjet1,drToLJ);
256  Double_t var2 = 0.;
257  //Double_t mJetR = 0.;
258  Double_t ptJetR = 0.;
259  if(jetR) {
260  //mJetR = jetR->M();
261  var2 = jetR->M();
262  ptJetR = jetR->Pt();
263  }
264  if(fSingleTrackEmb && vpe) {
265 
266  if(fFromTree){
267  Int_t exit = MatchEmbeddedConstituentWithParticleLevel(); //here is GetEntry
268 
269  if(exit>0) {
270  //mJetR = fVecP->M();
271  var2 = fVecP->M();
272  ptJetR = fVecP->Pt();
273 
274  }
275  } else{
276  //mJetR = vpe->M();
277  var2 = vpe->M();
278  ptJetR = vpe->Pt();
279  }
280  }
281 
282  if(fJetMassVarType==kRatMPt) {
283  if(ptJetR>0. || ptJetR<0.) var2 /= ptJetR;
284  }
285  fh3MSubPtTrueLeadPt[fCentBin]->Fill(var,ptJetR,jet1->MaxTrackPt());
286  fh3MTruePtTrueLeadPt[fCentBin]->Fill(var2,ptJetR,jet1->MaxTrackPt());
287  fh3PtTrueDeltaMLeadPt[fCentBin]->Fill(ptJetR,var-var2,jet1->MaxTrackPt());
288  if(var2>0.) fh3PtTrueDeltaMRelLeadPt[fCentBin]->Fill(ptJetR,(var-var2)/var2,jet1->MaxTrackPt());
289  //M sub;M true;#it{p}_{T,sub};#it{p}_{T,true};#it{p}_{T,lead trk}
290  if(fFromTree){
291  Double_t varsp[6] = {var,var2,ptjetS,ptJetR, fVecD->M(), fVecD->Pt()};
292  fhnMassResponse[fCentBin]->Fill(varsp);
293  } else {
294  Double_t varsp[5] = {var,var2,ptjetS,ptJetR,jetS->MaxTrackPt()};
295  fhnMassResponse[fCentBin]->Fill(varsp);
296  }
297  Double_t varsp1[8];
298  //#it{M}_{det,Const} - #it{M}_{part}; #it{p}_{T,det,Const} - #it{p}_{T,part}; #it{M}_{det,Const}; #it{M}_{part}; #it{p}_{T,det,Const}; #it{p}_{T,part}; #it{p}_{T,det,A}
299  varsp1[0] = var-var2;
300  varsp1[1] = ptjetS-ptJetR;
301  varsp1[2] = var;
302  varsp1[3] = var2;
303  varsp1[4] = ptjetS;
304  varsp1[5] = ptJetR;
305  varsp1[6] = ptjet1;
306  varsp1[7] = ptjet1 - ptJetR;
307 
308  fhnDeltaMass[fCentBin]->Fill(varsp1);
309 
310  //#it{M}_{det} - #it{M}_{part}; #it{p}_{T,det} - #it{p}_{T,part}; #it{M}_{det}; #it{M}_{unsub}; #it{p}_{T,det}; #it{p}_{T,unsub}; #rho ; #rho_{m}
311  Double_t varsp2[8] = {var-var2, ptjetS-ptJetR, var, mUnsubjet1, ptjetS, ptUnsubjet1, fRho, fRhoM};
312  fhnDeltaMassAndBkgInfo->Fill(varsp2);
313 
314  fhNconstit->Fill(jet1->GetNumberOfConstituents());
315  fhAreaJet ->Fill(jet1->Area());
316  }
317 
318  if(fCreateTree) {
319  fJet1Vec->SetPxPyPzE(jet1->Px(),jet1->Py(),jet1->Pz(),jet1->E());
320  if(jetS->Pt()>0.) fJetSubVec->SetPtEtaPhiM(jetS->Pt(),jetS->Eta(),jetS->Phi(),jetS->M());
321  else fJetSubVec->SetPtEtaPhiM(0.,0.,0.,0.);
322  fArea = (Float_t)jet1->Area();
323  fAreaPhi = (Float_t)jet1->AreaPhi();
324  fAreaEta = (Float_t)jet1->AreaEta();
325  fRho = (Float_t)jetCont->GetRhoVal();
326  fRhoM = (Float_t)jetCont->GetRhoMassVal();
327  fNConst = (Int_t)jet1->GetNumberOfTracks();
328  fTreeJetBkg->Fill();
329  }
330  } //jet1 loop
331  }//jetCont
332 
333  fhJet1vsJetTag->Fill(njet1, ntagjet2);
334 
335  return kTRUE;
336 }
337 
338 
TH3F ** fh3PtTrueDeltaMLeadPt
! true jet pT vs (Msub - Mtrue) vs LeadPt for matched jets
Double_t Area() const
Definition: AliEmcalJet.h:97
Double_t GetRhoVal() const
TH1F * fhptjetSMinusSingleTrack
! pT distribution of jets subtracting the pT of the embedded track
const TString & GetRhoName() const
AliEmcalJet * GetTaggedJet() const
Definition: AliEmcalJet.h:203
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:193
AliJetContainer * GetJetContainer(Int_t i=0) const
Int_t GetTagStatus() const
Definition: AliEmcalJet.h:204
Double_t Eta() const
Definition: AliEmcalJet.h:88
Double_t Py() const
Definition: AliEmcalJet.h:74
Double_t Phi() const
Definition: AliEmcalJet.h:84
Int_t GetLabel() const
Definition: AliEmcalJet.h:91
TH1F * fhAreaJet
! area of the matched jet
THnSparse * fhnDeltaMassAndBkgInfo
! DeltaM, DeltapT bkg-unsubtracted M and pT, rho and rhom
Int_t fCentBin
!event centrality bin
ClassImp(AliAnalysisTaskJetShapeConst) AliAnalysisTaskJetShapeConst
TH3F * fhJetSubMatchEtaPhiPt
! eta, phi, pt distribution of jet subtracted and matched
Double_t E() const
Definition: AliEmcalJet.h:86
TH2F * fhJet1vsJetTag
! N jet vs N jet tagged
TList * fOutput
!output list
UShort_t GetNumberOfConstituents() const
Definition: AliEmcalJet.h:107
Container for particles within the EMCAL framework.
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:106
TH3F ** fh3MSubPtTrueLeadPt
! subtracted jet mass vs true jet pT vs LeadPt for matched jets for matched jets
Double_t Px() const
Definition: AliEmcalJet.h:73
TH2F * fhJetEtaPhiOvl
! eta-phi distribution of the selected signal jets
AliParticleContainer * GetParticleContainer(Int_t i=0) const
const TString & GetRhoMassName() const
Int_t fContainerOverlap
jets (jetO) with a pT cut selection to reject overlapping embedded single track (used only in single ...
AliParticleContainer * GetParticleContainer() const
AliEmcalJet * GetLeadingJet(const char *opt="")
TLorentzVector * fVecD
! vector with detector level jet
Double_t AreaEta() const
Definition: AliEmcalJet.h:99
TH1F * fhNconstit
! number of constituents of the matched jets
Double_t fRadius
Radius that define overlap.
TLorentzVector * fJet1Vec
jet1(AA) vector
virtual AliVParticle * GetParticle(Int_t i=-1) const
Int_t GetNJets() const
Double_t MaxTrackPt() const
Definition: AliEmcalJet.h:122
TClonesArray * GetArray() const
ResponseReference fResponseReference
true axis of response matrix
Double_t GetRhoMassVal() const
Background fluctuation studies: dMdpT spectrum for PYTHIA and single track embedding for Constituent ...
TH3F ** fh3MSubPtRawDRMatch
! subtracted jet mass vs subtracted jet pT vs distance to leading Pb-Pb jet
THnSparse ** fhnDeltaMass
! deltaM vs deltapT
AliEmcalJet * GetNextAcceptJet()
Double_t DeltaR(const AliVParticle *part) const
Double_t AreaPhi() const
Definition: AliEmcalJet.h:100
THnSparse ** fhnMassResponse
! Msub vs Mtrue vs PtCorr vs PtTrue vs DR
Double_t fMinFractionShared
only fill histos for jets if shared fraction larger than X
Double_t Pt() const
Definition: AliEmcalJet.h:76
TH1F * fhNJetsSelEv
! number of selected signal jets per event
Bool_t fOverlap
activate the check on overlap between single particle embedded and jetO (jet with a pT of at least 5 ...
TH3F ** fh3MTruePtTrueLeadPt
! true jet mass vs true jet pT vs LeadPt for matched jets for matched jets
JetMassVarType fJetMassVarType
observable to use
TH2F ** fh2MSubMatch
! subtracted jet mass vs match index (0: no match; 1:match)
Int_t fMatch
1: matched to MC jet; 0: no match
Short_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:127
TTree * fTreeJetBkg
! tree with jet and bkg variables
TH2F ** fh2MSubPtRawAll
! subtracted jet mass vs subtracted jet pT
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
AliVParticle * GetEmbeddedConstituent(AliEmcalJet *jet)
Double_t Pz() const
Definition: AliEmcalJet.h:75
Int_t fNConst
N constituents in jet1.
Double_t GetFractionSharedPt(const AliEmcalJet *jet, AliParticleContainer *cont2=0x0) const
TH2F * fhRjetTrvspTj
! distance in R between each jetO and embedded single track (those below fRadius are rejected) ...
Int_t GetNAcceptedParticles() const
Int_t fContainerSub
subtracted jets to be analyzed
Bool_t fFromTree
Input embedding from tree.
TLorentzVector * fJetSubVec
subtracted AA jet vector
TLorentzVector * fJet2Vec
jet2(probe) vector
TLorentzVector * fVecP
! vector with particle level jet
Bool_t fCreateTree
create output tree
Double_t M() const
Definition: AliEmcalJet.h:87
void ResetCurrentID(Int_t i=-1)
Int_t fContainerNoEmb
subtracted jets from Pb-Pb only events
Bool_t fSingleTrackEmb
single track embedding
Background fluctuation studies: dMdpT spectrum for PYTHIA and single track embedding.
Container for jet within the EMCAL jet framework.
Int_t fContainerBase
jets to be analyzed
TH3F ** fh3PtTrueDeltaMRelLeadPt
! true jet pT vs (Msub - Mtrue)/Mtrue vs LeadPt for matched jets
Bool_t reject
AliEmcalJet * GetJet(Int_t i) const