AliPhysics  1168478 (1168478)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskSEDStarJets.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 //
16 //
17 // Base class for DStar in Jets Analysis
18 //
19 //-----------------------------------------------------------------------
20 // Author A.Grelli
21 // ERC-QGP Utrecht University - a.grelli@uu.nl
22 //-----------------------------------------------------------------------
23 
24 #include <TDatabasePDG.h>
25 #include <TParticle.h>
26 #include <TVector3.h>
27 #include "TROOT.h"
28 
31 #include "AliMCEvent.h"
32 #include "AliAODMCHeader.h"
33 #include "AliAODHandler.h"
34 #include "AliAnalysisManager.h"
35 #include "AliLog.h"
36 #include "AliAODVertex.h"
37 #include "AliAODJet.h"
38 #include "AliAODRecoDecay.h"
39 #include "AliAODRecoCascadeHF.h"
41 #include "AliESDtrack.h"
42 #include "AliAODMCParticle.h"
43 
47 
48 //__________________________________________________________________________
51  fEvents(0),
52  fchargeFrCorr(0),
53  fUseMCInfo(kTRUE),
54  fRequireNormalization(kTRUE),
55  fOutput(0),
56  fCuts(0),
57  ftrigger(0),
58  fPtPion(0),
59  fInvMass(0),
60  fRECOPtDStar(0),
61  fRECOPtBkg(0),
62  fDStar(0),
63  fDiff(0),
64  fDiffSideBand(0),
65  fDStarMass(0),
66  fPhi(0),
67  fPhiBkg(0),
68  fTrueDiff(0),
69  fResZ(0),
70  fResZBkg(0),
71  fEjet(0),
72  fPhijet(0),
73  fEtaJet(0),
74  theMCFF(0),
75  fDphiD0Dstar(0),
76  fPtJet(0)
77 {
78  //
80  //
81 }
82 //___________________________________________________________________________
84  AliAnalysisTaskSE(name),
85  fEvents(0),
86  fchargeFrCorr(0),
87  fUseMCInfo(kTRUE),
88  fRequireNormalization(kTRUE),
89  fOutput(0),
90  fCuts(0),
91  ftrigger(0),
92  fPtPion(0),
93  fInvMass(0),
94  fRECOPtDStar(0),
95  fRECOPtBkg(0),
96  fDStar(0),
97  fDiff(0),
98  fDiffSideBand(0),
99  fDStarMass(0),
100  fPhi(0),
101  fPhiBkg(0),
102  fTrueDiff(0),
103  fResZ(0),
104  fResZBkg(0),
105  fEjet(0),
106  fPhijet(0),
107  fEtaJet(0),
108  theMCFF(0),
109  fDphiD0Dstar(0),
110  fPtJet(0)
111 {
112  //
114  //
115  fCuts=cuts;
116  Info("AliAnalysisTaskSEDStarJets","Calling Constructor");
117 
118  DefineOutput(1,TList::Class()); // histos
119  DefineOutput(2,AliRDHFCutsDStartoKpipi::Class()); // my cuts
120 }
121 //___________________________________________________________________________
123  //
125  //
126 
127  Info("~AliAnalysisTaskSEDStarJets","Calling Destructor");
128  if (fOutput) {
129  delete fOutput;
130  fOutput = 0;
131  }
132 
133  if (fCuts) {
134  delete fCuts;
135  fCuts = 0;
136  }
137 
138  if (ftrigger) { delete ftrigger; ftrigger = 0;}
139  if (fPtPion) { delete fPtPion; fPtPion = 0;}
140  if (fInvMass) { delete fInvMass; fInvMass = 0;}
141  if (fRECOPtDStar) { delete fRECOPtDStar; fRECOPtDStar = 0;}
142  if (fRECOPtBkg) { delete fRECOPtBkg; fRECOPtBkg = 0;}
143  if (fDStar) { delete fDStar; fDStar = 0;}
144  if (fDiff) { delete fDiff; fDiff = 0;}
145  if (fDiffSideBand) { delete fDiffSideBand; fDiffSideBand = 0;}
146  if (fDStarMass) { delete fDStarMass; fDStarMass = 0;}
147  if (fPhi) { delete fPhi; fPhi = 0;}
148  if (fPhiBkg) { delete fPhiBkg; fPhiBkg = 0;}
149  if (fTrueDiff){ delete fTrueDiff; fTrueDiff = 0;}
150  if (fResZ) { delete fResZ; fResZ = 0;}
151  if (fResZBkg) { delete fResZBkg; fResZBkg = 0;}
152  if (fEjet) { delete fEjet; fEjet = 0;}
153  if (fPhijet) { delete fPhijet; fPhijet = 0;}
154  if (fEtaJet) { delete fEtaJet; fEtaJet = 0;}
155  if (theMCFF) { delete theMCFF; theMCFF = 0;}
156  if (fDphiD0Dstar) { delete fDphiD0Dstar; fDphiD0Dstar = 0;}
157  if (fPtJet) { delete fPtJet; fPtJet = 0;}
158 }
159 
160 //___________________________________________________________
162  //
164  //
165  if(fDebug > 1) printf("AnalysisTaskSEDStarJets::Init() \n");
167  // Post the cuts
168  PostData(2,copyfCuts);
169 
170  return;
171 }
172 
173 //_________________________________________________
175 {
177  if (!fInputEvent) {
178  Error("UserExec","NO EVENT FOUND!");
179  return;
180  }
181 
182  fEvents++;
183  AliInfo(Form("Event %d",fEvents));
184  if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents));
185 
186  // Load the event
187  AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
188 
189  TClonesArray *arrayDStartoD0pi=0;
190 
191  if(!aodEvent && AODEvent() && IsStandardAOD()) {
192  // In case there is an AOD handler writing a standard AOD, use the AOD
193  // event in memory rather than the input (ESD) event.
194  aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
195  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
196  // have to taken from the AOD event hold by the AliAODExtension
197  AliAODHandler* aodHandler = (AliAODHandler*)
198  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
199  if(aodHandler->GetExtensions()) {
200  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
201  AliAODEvent *aodFromExt = ext->GetAOD();
202  arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
203  }
204  } else {
205  arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
206  }
207 
208  if (!arrayDStartoD0pi){
209  AliInfo("Could not find array of HF vertices, skipping the event");
210  return;
211  }else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
212 
213  // fix for temporary bug in ESDfilter
214  // the AODs with null vertex pointer didn't pass the PhysSel
215  if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
216 
217  // Simulate a jet triggered sample
218  TClonesArray *arrayofJets = (TClonesArray*)aodEvent->GetJets();
219  if(aodEvent->GetNJets()<=0) return;
220 
221  // counters for efficiencies
222  Int_t icountReco = 0;
223 
224  // Normalization factor
226  ftrigger->Fill(1);
227  }
228 
229  //D* and D0 prongs needed to MatchToMC method
230  Int_t pdgDgDStartoD0pi[2]={421,211};
231  Int_t pdgDgD0toKpi[2]={321,211};
232 
233  Double_t max =0;
234  Double_t ejet = 0;
235  Double_t phiJet = 0;
236  Double_t etaJet = 0;
237  Double_t ptjet = 0;
238 
239  //loop over jets and consider only the leading jey in the event
240  for (Int_t iJets = 0; iJets<arrayofJets->GetEntriesFast(); iJets++) {
241  AliAODJet* jet = (AliAODJet*)arrayofJets->At(iJets);
242 
243  //jets variables
244  ejet = jet->E();
245 
246  if(ejet>max){
247  max = jet->E();
248  phiJet = jet->Phi();
249  etaJet = jet->Eta();
250  ptjet = jet->Pt();
251  }
252 
253  // fill energy, eta and phi of the jet
254  fEjet ->Fill(ejet);
255  fPhijet ->Fill(phiJet);
256  fEtaJet ->Fill(etaJet);
257  fPtJet->Fill(ptjet);
258  }
259 
260  //loop over D* candidates
261  for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
262 
263  // D* candidates
264  AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
265  AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
266 
267  Double_t finvM =-999;
268  Double_t finvMDStar =-999;
269  Double_t dPhi =-999;
270  Bool_t isDStar =0;
271  Int_t mcLabel = -9999;
272 
273  // find associated MC particle for D* ->D0toKpi
274  if(fUseMCInfo){
275  TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
276  if(!mcArray) {
277  printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
278  return;
279  }
280  mcLabel = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
281 
282  if(mcLabel>=0) isDStar = 1;
283  if(mcLabel>0){
284  Double_t zMC =-999;
285  AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel));
286  //fragmentation function in mc
287  zMC= FillMCFF(mcPart,mcArray,mcLabel);
288  if(zMC>0) theMCFF->Fill(zMC);
289  }
290  }
291 
292  // soft pion
293  AliAODTrack *track2 = (AliAODTrack*)dstarD0pi->GetBachelor();
294  Double_t pt = dstarD0pi->Pt();
295 
296  // track quality cuts
297  Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks
298  if(!isTkSelected) continue;
299 
300  // region of interest + topological cuts + PID
301  if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;
302  Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected
303  if (!isSelected) continue;
304 
305  // fill histos
306  finvM = dstarD0pi->InvMassD0();
307  fInvMass->Fill(finvM);
308 
309  //DStar invariant mass
310  finvMDStar = dstarD0pi->InvMassDstarKpipi();
311 
312  Double_t EGjet = 0;
313  Double_t dStarMom = dstarD0pi->P();
314  Double_t phiDStar = dstarD0pi->Phi();
315  Double_t phiD0 = theD0particle->Phi();
316  //check suggested by Federico
317  Double_t dPhiD0Dstar = phiD0 - phiDStar;
318 
319  dPhi = phiJet - phiDStar;
320 
321  //plot right dphi
322  if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
323  if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
324 
325  Double_t corrFactorCharge = (ejet/100)*20;
326  EGjet = ejet + corrFactorCharge;
327 
328  // fill D* candidates
329  Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
330  if(finvM >= (mPDGD0-0.05) && finvM <=(mPDGD0+0.05)){ // ~3 sigma (sigma=17MeV, conservative)
331 
332  if(isDStar == 1) {
333  fDphiD0Dstar->Fill(dPhiD0Dstar);
334  fDStarMass->Fill(finvMDStar);
335  fTrueDiff->Fill(finvMDStar-finvM);
336  }
337  if(isDStar == 0) fDphiD0Dstar->Fill(dPhiD0Dstar); // angle between D0 and D*
338 
339  fDStar->Fill(finvMDStar);
340  fDiff ->Fill(finvMDStar-finvM);
341 
342  Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
343  Double_t invmassDelta = dstarD0pi->DeltaInvMass();
344 
345  // now the dphi signal and the fragmentation function
346  if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<0.0019){
347  //fill candidates D* and soft pion reco pt
348 
349  fRECOPtDStar->Fill(pt);
350  fPtPion->Fill(track2->Pt());
351 
352  fPhi ->Fill(dPhi);
353 
354  Double_t jetCone = 0.4;
355  if(dPhi>=-jetCone && dPhi<=jetCone){ // evaluate in the near side inside UA1 radius
356  Double_t zFrag = (TMath::Cos(dPhi)*dStarMom)/EGjet;
357  fResZ->Fill(TMath::Abs(zFrag));
358  }
359  }
360  }
361  // evaluate side band background
362  SideBandBackground(finvM, finvMDStar, dStarMom, EGjet, dPhi);
363 
364  } // D* cand
365 
366 AliDebug(2, Form("Found %i Reco particles that are D*!!",icountReco));
367 
368 PostData(1,fOutput);
369 }
370 
371 //________________________________________ terminate ___________________________
373 {
377 
378  Info("Terminate"," terminate");
379  AliAnalysisTaskSE::Terminate();
380 
381  fOutput = dynamic_cast<TList*> (GetOutputData(1));
382  if (!fOutput) {
383  printf("ERROR: fOutput not available\n");
384  return;
385  }
386 
387  fDStarMass = dynamic_cast<TH1F*>(fOutput->FindObject("fDStarMass"));
388  fTrueDiff = dynamic_cast<TH1F*>(fOutput->FindObject("fTrueDiff"));
389  fInvMass = dynamic_cast<TH1F*>(fOutput->FindObject("fInvMass"));
390  fPtPion = dynamic_cast<TH1F*>(fOutput->FindObject("fPtPion "));
391  fDStar = dynamic_cast<TH1F*>(fOutput->FindObject("fDStar"));
392  fDiff = dynamic_cast<TH1F*>(fOutput->FindObject("fDiff"));
393  fDiffSideBand = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffSideBand"));
394  ftrigger = dynamic_cast<TH1F*>(fOutput->FindObject("ftrigger"));
395  fRECOPtDStar = dynamic_cast<TH1F*>(fOutput->FindObject("fRECOPtDStar"));
396  fRECOPtBkg = dynamic_cast<TH1F*>(fOutput->FindObject("fRECOPtBkg"));
397  fEjet = dynamic_cast<TH1F*>(fOutput->FindObject("fEjet"));
398  fPhijet = dynamic_cast<TH1F*>(fOutput->FindObject("fPhijet"));
399  fEtaJet = dynamic_cast<TH1F*>(fOutput->FindObject("fEtaJet"));
400  fPhi = dynamic_cast<TH1F*>(fOutput->FindObject("fPhi"));
401  fResZ = dynamic_cast<TH1F*>(fOutput->FindObject("fResZ"));
402  fResZBkg = dynamic_cast<TH1F*>(fOutput->FindObject("fResZBkg"));
403  fPhiBkg = dynamic_cast<TH1F*>(fOutput->FindObject("fPhiBkg"));
404  theMCFF = dynamic_cast<TH1F*>(fOutput->FindObject("theMCFF"));
405  fDphiD0Dstar = dynamic_cast<TH1F*>(fOutput->FindObject("fDphiD0Dstar"));
406  fPtJet = dynamic_cast<TH1F*>(fOutput->FindObject("fPtJet"));
407 
408 }
409 //___________________________________________________________________________
410 
413  Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
414 
415  //slot #1
416  OpenFile(1);
417  fOutput = new TList();
418  fOutput->SetOwner();
419  // define histograms
421 
422  return;
423 }
424 //___________________________________ hiostograms _______________________________________
425 
427 
429  fInvMass = new TH1F("invMass","Kpi invariant mass distribution",1500,.5,3.5);
430  fInvMass->SetStats(kTRUE);
431  fInvMass->GetXaxis()->SetTitle("GeV/c");
432  fInvMass->GetYaxis()->SetTitle("Entries");
433  fOutput->Add(fInvMass);
434 
435  fDStar = new TH1F("invMassDStar","DStar invariant mass after D0 cuts ",600,1.8,2.4);
436  fDStar->SetStats(kTRUE);
437  fDStar->GetXaxis()->SetTitle("GeV/c");
438  fDStar->GetYaxis()->SetTitle("Entries");
439  fOutput->Add(fDStar);
440 
441  fDiff = new TH1F("Diff","M(kpipi)-M(kpi)",750,0.1,0.2);
442  fDiff->SetStats(kTRUE);
443  fDiff->GetXaxis()->SetTitle("M(kpipi)-M(kpi) GeV/c^2");
444  fDiff->GetYaxis()->SetTitle("Entries");
445  fOutput->Add(fDiff);
446 
447  fDiffSideBand = new TH1F("DiffSide","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
448  fDiffSideBand->SetStats(kTRUE);
449  fDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
450  fDiffSideBand->GetYaxis()->SetTitle("Entries");
451  fOutput->Add(fDiffSideBand);
452 
453  fDStarMass = new TH1F("RECODStar2","RECO DStar invariant mass distribution",750,1.5,2.5);
454  fOutput->Add(fDStarMass);
455 
456  fTrueDiff = new TH1F("dstar","True Reco diff",750,0,0.2);
457  fOutput->Add(fTrueDiff);
458 
459  // trigger normalization
460  ftrigger = new TH1F("Normalization","Normalization factor for correlations",1,0,10);
461  ftrigger->SetStats(kTRUE);
462  fOutput->Add(ftrigger);
463 
464  //correlation fistograms
465  fPhi = new TH1F("phi","Delta phi between Jet axis and DStar ",25,-1.57,4.72);
466  fPhi->SetStats(kTRUE);
467  fPhi->GetXaxis()->SetTitle("#Delta #phi (rad)");
468  fPhi->GetYaxis()->SetTitle("Entries");
469  fOutput->Add(fPhi);
470 
471  fDphiD0Dstar = new TH1F("phiD0Dstar","Delta phi between D0 and DStar ",1000,-6.5,6.5);
472  fOutput->Add(fDphiD0Dstar);
473 
474  fPhiBkg = new TH1F("phiBkg","Delta phi between Jet axis and DStar background ",25,-1.57,4.72);
475  fPhiBkg->SetStats(kTRUE);
476  fPhiBkg->GetXaxis()->SetTitle("#Delta #phi (rad)");
477  fPhiBkg->GetYaxis()->SetTitle("Entries");
478  fOutput->Add(fPhiBkg);
479 
480  fRECOPtDStar = new TH1F("RECODStar1","RECO DStar pt distribution",30,0,30);
481  fRECOPtDStar->SetStats(kTRUE);
482  fRECOPtDStar->SetLineColor(2);
483  fRECOPtDStar->GetXaxis()->SetTitle("GeV/c");
484  fRECOPtDStar->GetYaxis()->SetTitle("Entries");
485  fOutput->Add(fRECOPtDStar);
486 
487  fRECOPtBkg = new TH1F("RECOptBkg","RECO pt distribution side bands",30,0,30);
488  fRECOPtBkg->SetStats(kTRUE);
489  fRECOPtBkg->SetLineColor(2);
490  fRECOPtBkg->GetXaxis()->SetTitle("GeV/c");
491  fRECOPtBkg->GetYaxis()->SetTitle("Entries");
492  fOutput->Add(fRECOPtBkg);
493 
494  fPtPion = new TH1F("pionpt","Primary pions candidates pt ",500,0,10);
495  fPtPion->SetStats(kTRUE);
496  fPtPion->GetXaxis()->SetTitle("GeV/c");
497  fPtPion->GetYaxis()->SetTitle("Entries");
498  fOutput->Add(fPtPion);
499 
500  // jet related fistograms
501  fEjet = new TH1F("ejet", "UA1 algorithm jet energy distribution",1000,0,500);
502  fPhijet = new TH1F("Phijet","UA1 algorithm jet #phi distribution", 200,-7,7);
503  fEtaJet = new TH1F("Etajet","UA1 algorithm jet #eta distribution", 200,-7,7);
504  fPtJet = new TH1F("PtJet", "UA1 algorithm jet Pt distribution",1000,0,500);
505  fOutput->Add(fEjet);
506  fOutput->Add(fPhijet);
507  fOutput->Add(fEtaJet);
508  fOutput->Add(fPtJet);
509 
510  theMCFF = new TH1F("FragFuncMC","Fragmentation function in MC for FC ",100,0,10);
511  fResZ = new TH1F("FragFunc","Fragmentation function ",50,0,1);
512  fResZBkg = new TH1F("FragFuncBkg","Fragmentation function background",50,0,1);
513  fOutput->Add(theMCFF);
514  fOutput->Add(fResZ);
515  fOutput->Add(fResZBkg);
516 
517  return kTRUE;
518 }
519 
520 //______________________________ side band background for D*___________________________________
521 
523 
527 
528  if((invM>=1.7 && invM<=1.8) || (invM>=1.92 && invM<=2.19)){
529  fDiffSideBand->Fill(invMDStar-invM); // M(Kpipi)-M(Kpi) side band background
530  if ((invMDStar-invM)<=0.14732 && (invMDStar-invM)>=0.14352) {
531  fPhiBkg->Fill(dPhi);
532  fRECOPtBkg->Fill(dStarMomBkg);
533  if(dPhi>=-0.4 && dPhi<=0.4){ // evaluate in the near side
534  Double_t zFragBkg = (TMath::Cos(dPhi)*dStarMomBkg)/EGjet;
535  fResZBkg->Fill(TMath::Abs(zFragBkg));
536  }
537  }
538  }
539 }
540 
541 //_____________________________________________________________________________________________-
542 double AliAnalysisTaskSEDStarJets::FillMCFF(AliAODMCParticle* mcPart, TClonesArray* mcArray, Int_t mcLabel){
547  Double_t zMC2 =-999;
548 
549  Double_t leading =0;
550  Double_t PartE = 0;
551  Double_t PhiLeading = -999;
552  Double_t EtaLeading = -999;
553  Double_t PtLeading = -999;
554  Int_t counter =-999;
555 
556  if (!mcPart) return zMC2;
557  if (!mcArray) return zMC2;
558 
559  //find leading particle
560  for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
561  AliAODMCParticle* Part = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
562  if (!Part) {
563  AliWarning("MC Particle not found in tree, skipping");
564  continue;
565  }
566 
567  // remove quarks and the leading particle (it will be counted later)
568  if(iPart == mcLabel) continue;
569  if(iPart <= 8) continue;
570 
571  //remove resonances not directly detected in detector
572  Int_t PDGCode = Part->GetPdgCode();
573 
574  // be sure the particle reach the detector
575  Double_t x = Part->Xv();
576  Double_t y = Part->Yv();
577  Double_t z = Part->Zv();
578 
579  if(TMath::Abs(PDGCode)== 2212 && x<3 && y<3) continue;
580  if(TMath::Abs(x)>30 || TMath::Abs(y)>30 || TMath::Abs(z)>30 ) continue;
581  if(TMath::Abs(PDGCode)!=211 && TMath::Abs(PDGCode)!=321 && TMath::Abs(PDGCode)!=11 && TMath::Abs(PDGCode)!=13 && TMath::Abs(PDGCode)!=2212) continue;
582 
583  Int_t daug0 = -999;
584  Double_t xd =-999;
585  Double_t yd =-999;
586  Double_t zd =-999;
587 
588  daug0 = Part->GetDaughter(0);
589 
590  if(daug0>=0){
591  AliAODMCParticle* tdaug = dynamic_cast<AliAODMCParticle*>(mcArray->At(daug0));
592  if(tdaug){
593  xd = tdaug->Xv();
594  yd = tdaug->Yv();
595  zd = tdaug->Zv();
596  }
597  }
598  if(TMath::Abs(xd)<3 || TMath::Abs(yd)<3) continue;
599 
600  Bool_t AliceAcc = (TMath::Abs(Part->Eta()) <= 0.9);
601  if(!AliceAcc) continue;
602 
603  PartE = Part->E();
604 
605  if(PartE>leading){
606  leading = Part->E();
607  PhiLeading = Part->Phi();
608  EtaLeading = Part->Eta();
609  PtLeading = Part->Pt();
610  counter = iPart;
611  }
612 
613  }
614 
615  Double_t jetEnergy = 0;
616 
617  //reconstruct the jet
618  for (Int_t iiPart=0; iiPart<mcArray->GetEntriesFast(); iiPart++) {
619  AliAODMCParticle* tPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iiPart));
620  if (!tPart) {
621  AliWarning("MC Particle not found in tree, skipping");
622  continue;
623  }
624  // remove quarks and the leading particle (it will be counted later)
625  if(iiPart == counter) continue; // do not count again the leading particle
626  if(iiPart == mcLabel) continue;
627  if(iiPart <= 8) continue;
628 
629  //remove resonances not directly detected in detector
630  Int_t PDGCode = tPart->GetPdgCode();
631 
632  // be sure the particle reach the detector
633  Double_t x = tPart->Xv();
634  Double_t y = tPart->Yv();
635  Double_t z = tPart->Zv();
636 
637  if(TMath::Abs(PDGCode)== 2212 && (x<3 && y<3)) continue;
638  if(TMath::Abs(x)>30 || TMath::Abs(y)>30 || TMath::Abs(z)>30 ) continue; // has to be generated at least in the silicon tracker or beam pipe
639  if(TMath::Abs(PDGCode)!=211 && TMath::Abs(PDGCode)!=321 && TMath::Abs(PDGCode)!=11 && TMath::Abs(PDGCode)!=13 && TMath::Abs(PDGCode)!=2212) continue;
640 
641 
642  Int_t daug0 = -999;
643  Double_t xd =-999;
644  Double_t yd =-999;
645  Double_t zd =-999;
646 
647  daug0 = tPart->GetDaughter(0);
648 
649  if(daug0>=0){
650  AliAODMCParticle* tdaug = dynamic_cast<AliAODMCParticle*>(mcArray->At(daug0));
651  if(tdaug){
652  xd = tdaug->Xv();
653  yd = tdaug->Yv();
654  zd = tdaug->Zv();
655  }
656  }
657  if(TMath::Abs(xd)<3 && TMath::Abs(yd)<3) continue;
658  //remove particles not in ALICE acceptance
659  if(tPart->Pt()<0.07) continue;
660  Bool_t AliceAcc = (TMath::Abs(tPart->Eta()) <= 0.9);
661  if(!AliceAcc) continue;
662 
663  Double_t EtaMCp = tPart->Eta();
664  Double_t PhiMCp = tPart->Phi();
665 
666  Double_t DphiMClead = PhiLeading-PhiMCp;
667 
668  if(DphiMClead<=-(TMath::Pi())/2) DphiMClead = DphiMClead+2*(TMath::Pi());
669  if(DphiMClead>(3*(TMath::Pi()))/2) DphiMClead = DphiMClead-2*(TMath::Pi());
670 
671  Double_t deta = (EtaLeading-EtaMCp);
672  //cone radius
673  Double_t radius = TMath::Sqrt((DphiMClead*DphiMClead)+(deta*deta));
674 
675  if(radius>0.4) continue; // in the jet cone
676  if(tPart->Charge()==0) continue; // only charged fraction
677 
678  jetEnergy = jetEnergy+(tPart->E());
679  }
680 
681  jetEnergy = jetEnergy + leading;
682 
683  // delta phi D*, jet axis
684  Double_t dPhi = PhiLeading - (mcPart->Phi());
685 
686  //plot right dphi
687  if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
688  if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
689 
690  zMC2 = (TMath::Cos(dPhi)*(mcPart->P()))/jetEnergy;
691 
692  return zMC2;
693 }
double Double_t
Definition: External.C:58
Double_t DeltaInvMass() const
virtual void UserExec(Option_t *option)
Double_t YDstar() const
Int_t MatchToMC(Int_t pdgabs, Int_t pdgabs2prong, Int_t *pdgDg, Int_t *pdgDg2prong, TClonesArray *mcArray, Bool_t isV0=kFALSE) const
AliRDHFCutsDStartoKpipi * fCuts
char Char_t
Definition: External.C:18
Bool_t fUseMCInfo
Charge fraction correction UA1 algorithm.
Bool_t DefineHistoFroAnalysis()
inizializations
int Int_t
Definition: External.C:63
AliAODTrack * GetBachelor() const
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel, AliAODEvent *aod)
double FillMCFF(AliAODMCParticle *mcPart, TClonesArray *mcArray, Int_t mcLabel)
MC FF.
virtual Bool_t IsInFiducialAcceptance(Double_t pt, Double_t y) const
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
Double_t InvMassD0() const
const char Option_t
Definition: External.C:48
Double_t InvMassDstarKpipi() const
bool Bool_t
Definition: External.C:53
AliAODRecoDecayHF2Prong * Get2Prong() const
void SideBandBackground(Double_t finvM, Double_t finvMDStar, Double_t dStarMomBkg, Double_t fejet, Double_t ejet)
side band background eval
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65