AliPhysics  e59a9ba (e59a9ba)
 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 
130  Int_t njet1 = 0, ntagjet2 = 0;
131 
132  AliDebug(11,Form("NJets Incl: %d Csub: %d",jetCont->GetNJets(),jetContS->GetNJets()));
133  if(jetCont) {
134  jetCont->ResetCurrentID();
135  while((jet1 = jetCont->GetNextAcceptJet())) {
136  njet1++;
137  jet2 = NULL;
138  jetS = NULL;
139  if(jet1->GetTagStatus()<1 || !jet1->GetTaggedJet())
140  continue;
141 
142  ntagjet2++;
143  //print constituents of different jet containers
144  //jet1
145 
146  for(Int_t i=0; i<jet1->GetNumberOfTracks(); i++) {
147  AliVParticle *vp = static_cast<AliVParticle*>(jet1->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
148  // if (vp->TestBits(TObject::kBitMask) != (Int_t)(TObject::kBitMask) ) continue;
149  //Int_t lab = TMath::Abs(vp->GetLabel());
150  if(vp) fhpTTracksJet1 -> Fill(vp->Pt());
151  }
152 
153  //Get constituent subtacted version of jet
154 
155  Int_t ifound = 0;
156  Int_t ilab = -1;
157  for(Int_t i = 0; i<jetContS->GetNJets(); i++) {
158  //if(ifound==1) continue;
159  jetS = jetContS->GetJet(i);
160  if(!jetS) continue;
161  if(jetS->GetLabel()==jet1->GetLabel()) {
162  ifound++;
163  if(ifound==1) ilab = i;
164  }
165  }
166  if(ifound>1) AliDebug(2,Form("Found %d partners",ifound));
167  if(ifound==0) jetS = 0x0;
168  else jetS = jetContS->GetJet(ilab);
169  if(!jetS) continue;
170 
171  Double_t mjetS = jetS->M();
172  Double_t mUnsubjet1 = jet1->M();
173  Double_t ptjet1 = jet1->Pt()-jetCont->GetRhoVal()*jet1->Area();
174  Double_t ptjetS = jetS->Pt();
175  Double_t ptUnsubjet1 = jet1->Pt();
176  Double_t var = mjetS;
177  //Double_t ptjetSMinusEmbTrpt = ptjetS;
178  Double_t ptjetSMinusEmbTrpt = ptUnsubjet1;
179 
180  Double_t rhofactor = GetRhoFactor(ptjet1, fRho); // the pt used in the mass task is the one subtracted with area based method
181 
182  if(fJetMassVarType==kRatMPt) {
183  if(ptjetS>0. || ptjetS<0.) var = mjetS/ptjetS;
184  else var = -999.;
185  }
186 
187  //Fill histograms for all AA jets
188  fh2MSubPtRawAll[fCentBin]->Fill(var,ptjetS);
189 
190  Double_t fraction = 0.;
191  fMatch = 0;
192  fJet2Vec->SetPtEtaPhiM(0.,0.,0.,0.);
193  //Printf("EMB task: pT jet %d = %f", njet1, jet1->Pt());
194  if(fSingleTrackEmb) {
195  vpe = GetEmbeddedConstituent(jet1);
196  if(vpe){
197  Bool_t reject = kFALSE;
198  ptjetSMinusEmbTrpt -= vpe->Pt();
199  fhptjetSMinusSingleTrack->Fill(ptjetSMinusEmbTrpt);
200  //Printf("EMB task: pT jet matched = %f -> %f", jet1->Pt(), ptjetSMinusEmbTrpt);
201  if(fOverlap){
202  Int_t Njets = jetContO->GetNAcceptedJets();
203  fhNJetsSelEv->Fill(Njets);
204  jetContO->ResetCurrentID();
205  while((jetO = jetContO->GetNextAcceptJet())){
206  //print constituents of different jet containers
207  //jetO
208  //Printf("N particle %d",jetO->GetNumberOfTracks());
209  for(Int_t i=0; i<jetO->GetNumberOfTracks(); i++) {
210  AliVParticle* vp = static_cast<AliVParticle*>(jetO->TrackAt(i, jetContO->GetParticleContainer()->GetArray()));
211  // if (vp->TestBits(TObject::kBitMask) != (Int_t)(TObject::kBitMask) ) continue;
212  //Int_t lab = TMath::Abs(vp->GetLabel());
213  if(vp) fhpTTracksJetO -> Fill(vp->Pt());
214  }
215  Double_t deltaR = jetO->DeltaR(vpe);
216  fhRjetTrvspTj->Fill(deltaR, jetO->Pt());
217  fhJetEtaPhiOvl->Fill(jetO->Eta(), jetO->Phi());
218  if(deltaR < fRadius) {
219  reject = kTRUE;
220  break;
221  }
222  }
223  }
224  if(!reject) {
225  fJet2Vec->SetPxPyPzE(vpe->Px(),vpe->Py(),vpe->Pz(),vpe->E());
226  fMatch = 1;
227  }
228  }
229  } else {
230  jet2 = jet1->ClosestJet();
231  fraction = jetCont->GetFractionSharedPt(jet1);
232  fMatch = 1;
233  if(fMinFractionShared>0.) {
234  if(fraction>fMinFractionShared) {
235  fJet2Vec->SetPxPyPzE(jet2->Px(),jet2->Py(),jet2->Pz(),jet2->E());
236  fMatch = 1;
237 
238  //choose jet type for true axis of response matrix
240  jetR = jet2;
241  else if(fResponseReference==kPart)
242  jetR = jet2->GetTaggedJet();
243  } else
244  fMatch = 0;
245  }
246  }
247 
248  // if(fMatch==1 && jet2->GetTagStatus()>0) jet2T = jet2->GetTaggedJet();
249 
250  //Fill histograms for matched jets
251  fh2MSubMatch[fCentBin]->Fill(var,fMatch);
252 
253  if(fMatch==1) {
254  fhJetSubMatchEtaPhiPt->Fill(jetS->Eta(), jetS->Phi(), ptjetS);
255  Double_t drToLJ = -1.;
256  if(jetL) drToLJ = jet1->DeltaR(jetL);
257  if(fSingleTrackEmb && vpe)
258  drToLJ = jet1->DeltaR(vpe);
259  fh3MSubPtRawDRMatch[fCentBin]->Fill(var,ptjet1,drToLJ);
260  Double_t var2 = 0.;
261  //Double_t mJetR = 0.;
262  Double_t ptJetR = 0.;
263  if(jetR) {
264  //mJetR = jetR->M();
265  var2 = jetR->M();
266  ptJetR = jetR->Pt();
267  }
268  if(fSingleTrackEmb && vpe) {
269 
270  if(fFromTree){
271  Int_t exit = MatchEmbeddedConstituentWithParticleLevel(); //here is GetEntry
272 
273  if(exit>0) {
274  //mJetR = fVecP->M();
275  var2 = fVecP->M();
276  ptJetR = fVecP->Pt();
277 
278  }
279  } else{
280  //mJetR = vpe->M();
281  var2 = vpe->M();
282  ptJetR = vpe->Pt();
283  }
284  }
285 
286  if(fJetMassVarType==kRatMPt) {
287  if(ptJetR>0. || ptJetR<0.) var2 /= ptJetR;
288  }
289  fh3MSubPtTrueLeadPt[fCentBin]->Fill(var,ptJetR,jet1->MaxTrackPt());
290  fh3MTruePtTrueLeadPt[fCentBin]->Fill(var2,ptJetR,jet1->MaxTrackPt());
291  fh3PtTrueDeltaMLeadPt[fCentBin]->Fill(ptJetR,var-var2,jet1->MaxTrackPt());
292  if(var2>0.) fh3PtTrueDeltaMRelLeadPt[fCentBin]->Fill(ptJetR,(var-var2)/var2,jet1->MaxTrackPt());
293  //M sub;M true;#it{p}_{T,sub};#it{p}_{T,true};#it{p}_{T,lead trk}
294  //if(fRho > 4 && ptjet1 > 40) Printf("Bin rho is (%f) -> %f", fRho, rhofactor);
295  fRhoFactorQA->Fill(rhofactor);
296  if(fFromTree){
297  // Mass sub; Mass true;#it{p}_{T,sub};#it{p}_{T,true};%s (emb, det); #it{p}_{T,emb det}; #rho; #rho_{m};
298  Double_t varsp[8] = {var,var2,ptjetS,ptJetR, fVecD->M(), fVecD->Pt(), fRho, fRhoM};
299  fhnMassResponse[fCentBin]->Fill(varsp, rhofactor);
300  } else {
301  Double_t varsp[7] = {var,var2,ptjetS,ptJetR,jetS->MaxTrackPt(), fRho, fRhoM};
302  fhnMassResponse[fCentBin]->Fill(varsp, rhofactor);
303  }
304  Double_t varsp1[8];
305  //#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}
306  varsp1[0] = var-var2;
307  varsp1[1] = ptjetS-ptJetR;
308  varsp1[2] = var;
309  varsp1[3] = var2;
310  varsp1[4] = ptjetS;
311  varsp1[5] = ptJetR;
312  varsp1[6] = ptjet1;
313  varsp1[7] = ptjet1 - ptJetR;
314 
315  fhnDeltaMass[fCentBin]->Fill(varsp1);
316 
317  //#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}
318  Double_t varsp2[8] = {var-var2, ptjetS-ptJetR, var, mUnsubjet1, ptjetS, ptUnsubjet1, fRho, fRhoM};
319  fhnDeltaMassAndBkgInfo->Fill(varsp2);
320 
321  fhNconstit->Fill(jet1->GetNumberOfConstituents());
322  fhAreaJet ->Fill(jet1->Area());
323  }
324 
325  if(fCreateTree) {
326  fJet1Vec->SetPxPyPzE(jet1->Px(),jet1->Py(),jet1->Pz(),jet1->E());
327  if(jetS->Pt()>0.) fJetSubVec->SetPtEtaPhiM(jetS->Pt(),jetS->Eta(),jetS->Phi(),jetS->M());
328  else fJetSubVec->SetPtEtaPhiM(0.,0.,0.,0.);
329  fArea = (Float_t)jet1->Area();
330  fAreaPhi = (Float_t)jet1->AreaPhi();
331  fAreaEta = (Float_t)jet1->AreaEta();
332  fRho = (Float_t)jetCont->GetRhoVal();
333  fRhoM = (Float_t)jetCont->GetRhoMassVal();
334  fNConst = (Int_t)jet1->GetNumberOfTracks();
335  fTreeJetBkg->Fill();
336  }
337  } //jet1 loop
338  }//jetCont
339 
340  fhJet1vsJetTag->Fill(njet1, ntagjet2);
341 
342  return kTRUE;
343 }
344 
345 
TH3F ** fh3PtTrueDeltaMLeadPt
! true jet pT vs (Msub - Mtrue) vs LeadPt for matched jets
TH1D * fRhoFactorQA
! actual rho weighting factors used
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
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
Double_t GetRhoFactor(Double_t ptdet, Double_t rhoCurrent)
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
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)
AliEmcalList * fOutput
!output list
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
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