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