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