AliPhysics  95775ff (95775ff)
 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  if(fJetMassVarType==kRatMPt) {
181  if(ptjetS>0. || ptjetS<0.) var = mjetS/ptjetS;
182  else var = -999.;
183  }
184 
185  //Fill histograms for all AA jets
186  fh2MSubPtRawAll[fCentBin]->Fill(var,ptjetS);
187 
188  Double_t fraction = 0.;
189  fMatch = 0;
190  fJet2Vec->SetPtEtaPhiM(0.,0.,0.,0.);
191  //Printf("EMB task: pT jet %d = %f", njet1, jet1->Pt());
192  if(fSingleTrackEmb) {
193  vpe = GetEmbeddedConstituent(jet1);
194  if(vpe){
195  Bool_t reject = kFALSE;
196  ptjetSMinusEmbTrpt -= vpe->Pt();
197  fhptjetSMinusSingleTrack->Fill(ptjetSMinusEmbTrpt);
198  //Printf("EMB task: pT jet matched = %f -> %f", jet1->Pt(), ptjetSMinusEmbTrpt);
199  if(fOverlap){
200  Int_t Njets = jetContO->GetNAcceptedJets();
201  fhNJetsSelEv->Fill(Njets);
202  jetContO->ResetCurrentID();
203  while((jetO = jetContO->GetNextAcceptJet())){
204  //print constituents of different jet containers
205  //jetO
206  //Printf("N particle %d",jetO->GetNumberOfTracks());
207  for(Int_t i=0; i<jetO->GetNumberOfTracks(); i++) {
208  AliVParticle* vp = static_cast<AliVParticle*>(jetO->TrackAt(i, jetContO->GetParticleContainer()->GetArray()));
209  // if (vp->TestBits(TObject::kBitMask) != (Int_t)(TObject::kBitMask) ) continue;
210  //Int_t lab = TMath::Abs(vp->GetLabel());
211  if(vp) fhpTTracksJetO -> Fill(vp->Pt());
212  }
213  Double_t deltaR = jetO->DeltaR(vpe);
214  fhRjetTrvspTj->Fill(deltaR, jetO->Pt());
215  fhJetEtaPhiOvl->Fill(jetO->Eta(), jetO->Phi());
216  if(deltaR < fRadius) {
217  reject = kTRUE;
218  break;
219  }
220  }
221  }
222  if(!reject) {
223  fJet2Vec->SetPxPyPzE(vpe->Px(),vpe->Py(),vpe->Pz(),vpe->E());
224  fMatch = 1;
225  }
226  }
227  } else {
228  jet2 = jet1->ClosestJet();
229  fraction = jetCont->GetFractionSharedPt(jet1);
230  fMatch = 1;
231  if(fMinFractionShared>0.) {
232  if(fraction>fMinFractionShared) {
233  fJet2Vec->SetPxPyPzE(jet2->Px(),jet2->Py(),jet2->Pz(),jet2->E());
234  fMatch = 1;
235 
236  //choose jet type for true axis of response matrix
238  jetR = jet2;
239  else if(fResponseReference==kPart)
240  jetR = jet2->GetTaggedJet();
241  } else
242  fMatch = 0;
243  }
244  }
245 
246  // if(fMatch==1 && jet2->GetTagStatus()>0) jet2T = jet2->GetTaggedJet();
247 
248  //Fill histograms for matched jets
249  fh2MSubMatch[fCentBin]->Fill(var,fMatch);
250 
251  if(fMatch==1) {
252  fhJetSubMatchEtaPhiPt->Fill(jetS->Eta(), jetS->Phi(), ptjetS);
253  Double_t drToLJ = -1.;
254  if(jetL) drToLJ = jet1->DeltaR(jetL);
255  if(fSingleTrackEmb && vpe)
256  drToLJ = jet1->DeltaR(vpe);
257  fh3MSubPtRawDRMatch[fCentBin]->Fill(var,ptjet1,drToLJ);
258  Double_t var2 = 0.;
259  //Double_t mJetR = 0.;
260  Double_t ptJetR = 0.;
261  if(jetR) {
262  //mJetR = jetR->M();
263  var2 = jetR->M();
264  ptJetR = jetR->Pt();
265  }
266  if(fSingleTrackEmb && vpe) {
267 
268  if(fFromTree){
269  Int_t exit = MatchEmbeddedConstituentWithParticleLevel(); //here is GetEntry
270 
271  if(exit>0) {
272  //mJetR = fVecP->M();
273  var2 = fVecP->M();
274  ptJetR = fVecP->Pt();
275 
276  }
277  } else{
278  //mJetR = vpe->M();
279  var2 = vpe->M();
280  ptJetR = vpe->Pt();
281  }
282  }
283 
284  if(fJetMassVarType==kRatMPt) {
285  if(ptJetR>0. || ptJetR<0.) var2 /= ptJetR;
286  }
287  fh3MSubPtTrueLeadPt[fCentBin]->Fill(var,ptJetR,jet1->MaxTrackPt());
288  fh3MTruePtTrueLeadPt[fCentBin]->Fill(var2,ptJetR,jet1->MaxTrackPt());
289  fh3PtTrueDeltaMLeadPt[fCentBin]->Fill(ptJetR,var-var2,jet1->MaxTrackPt());
290  if(var2>0.) fh3PtTrueDeltaMRelLeadPt[fCentBin]->Fill(ptJetR,(var-var2)/var2,jet1->MaxTrackPt());
291  //M sub;M true;#it{p}_{T,sub};#it{p}_{T,true};#it{p}_{T,lead trk}
292 
293  if(!fCreateTree){
294  if(fFromTree){
295  // Mass sub; Mass true;#it{p}_{T,sub};#it{p}_{T,true};%s (emb, det); #it{p}_{T,emb det}; #rho; #rho_{m};
296  Double_t varsp[10] = {var,var2,ptjetS,ptJetR, fVecD->M(), fVecD->Pt(), fRho, fRhoM, mUnsubjet1, ptUnsubjet1};
297  fhnMassResponse[fCentBin]->Fill(varsp);
298  } else {
299  Double_t varsp[9] = {var,var2,ptjetS,ptJetR,jetS->MaxTrackPt(), fRho, fRhoM, mUnsubjet1, ptUnsubjet1};
300  fhnMassResponse[fCentBin]->Fill(varsp);
301  }
302 
303  Double_t varsp1[8];
304  //#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}
305  varsp1[0] = var-var2;
306  varsp1[1] = ptjetS-ptJetR;
307  varsp1[2] = var;
308  varsp1[3] = var2;
309  varsp1[4] = ptjetS;
310  varsp1[5] = ptJetR;
311  varsp1[6] = ptjet1;
312  varsp1[7] = ptjet1 - ptJetR;
313 
314  //fhnDeltaMass[fCentBin]->Fill(varsp1);
315 
316  //#it{M}_{det} - #it{M}_{part}; #it{p}_{T,det} - #it{p}_{T,part}; #it{M}_{unsub} - #it{M}_{part}; #it{p}_{T,unsub} - #it{p}_{T,part}; #it{M}_{det}; #it{M}_{unsub}; #it{p}_{T,det}; #it{p}_{T,unsub}; #rho ; #rho_{m}
317 
318  Double_t varsp2[10] = {var-var2, ptjetS-ptJetR, mUnsubjet1 - var2, ptUnsubjet1 - ptJetR, var2, mUnsubjet1, ptjetS, ptUnsubjet1, fRho, fRhoM};
319  fhnDeltaMassAndBkgInfo->Fill(varsp2);
320 
321  Double_t varsp3[6] = {(var-var2)/var2, (ptjetS-ptJetR)/ptJetR, var2, ptJetR, fRho, fRhoM};
322  fhnResolution->Fill(varsp3);
323  }
324  fhNconstit->Fill(jet1->GetNumberOfConstituents());
325  fhAreaJet ->Fill(jet1->Area());
326  }
327 
328  if(fCreateTree) {
329  fJet1Vec->SetPxPyPzE(jet1->Px(),jet1->Py(),jet1->Pz(),jet1->E());
330  if(jetS->Pt()>0.) fJetSubVec->SetPtEtaPhiM(jetS->Pt(),jetS->Eta(),jetS->Phi(),jetS->M());
331  else fJetSubVec->SetPtEtaPhiM(0.,0.,0.,0.);
332  fArea = (Float_t)jet1->Area();
333  fAreaPhi = (Float_t)jet1->AreaPhi();
334  fAreaEta = (Float_t)jet1->AreaEta();
335  fRho = (Float_t)jetCont->GetRhoVal();
336  fRhoM = (Float_t)jetCont->GetRhoMassVal();
337  fNConst = (Int_t)jet1->GetNumberOfTracks();
338  fTreeJetBkg->Fill();
339  }
340  } //jet1 loop
341  }//jetCont
342 
343  fhJet1vsJetTag->Fill(njet1, ntagjet2);
344 
345  return kTRUE;
346 }
347 
348 
TH3F ** fh3PtTrueDeltaMLeadPt
! true jet pT vs (Msub - Mtrue) vs LeadPt for matched jets
Double_t Area() const
Definition: AliEmcalJet.h:123
Double_t GetRhoVal() const
TH1F * fhptjetSMinusSingleTrack
! pT distribution of jets subtracting the pT of the embedded track
const TString & GetRhoName() const
double Double_t
Definition: External.C:58
AliEmcalJet * GetTaggedJet() const
Definition: AliEmcalJet.h:229
Definition: External.C:236
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:219
AliJetContainer * GetJetContainer(Int_t i=0) const
Int_t GetTagStatus() const
Definition: AliEmcalJet.h:230
Double_t Eta() const
Definition: AliEmcalJet.h:114
Double_t Py() const
Definition: AliEmcalJet.h:100
Double_t Phi() const
Definition: AliEmcalJet.h:110
Int_t GetLabel() const
Definition: AliEmcalJet.h:117
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:112
TH2F * fhJet1vsJetTag
! N jet vs N jet tagged
UShort_t GetNumberOfConstituents() const
Definition: AliEmcalJet.h:133
Container for particles within the EMCAL framework.
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:132
TH3F ** fh3MSubPtTrueLeadPt
! subtracted jet mass vs true jet pT vs LeadPt for matched jets for matched jets
Double_t Px() const
Definition: AliEmcalJet.h:99
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:125
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
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:148
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
AliEmcalJet * GetNextAcceptJet()
Double_t DeltaR(const AliVParticle *part) const
Double_t AreaPhi() const
Definition: AliEmcalJet.h:126
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:102
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
THnSparse * fhnResolution
! Contains mass and pT resolution
Int_t fMatch
1: matched to MC jet; 0: no match
Short_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:153
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:101
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) ...
bool Bool_t
Definition: External.C:53
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:113
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