New task for D* pt-dep analysis (A. Grelli)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSEDStarSpectra.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 Analysis
18 //
19 //  Side Band and like sign background are implemented in the macro
20 //
21 //  Utrech Optimized cuts used and PID is on request (flag in che .C) 
22 //  The D* spectra study is done in pt bins:
23 //
24 //  [0,1] [1,2] [2,3] [3,5] [5,8] [8,14]
25 //
26 //-----------------------------------------------------------------------
27 //
28 //                         Author A.Grelli 
29 //              ERC-QGP Utrecht University - a.grelli@uu.nl
30 //           Contributors: C.Ivan  - Utrecht University
31 //
32 //-----------------------------------------------------------------------
33
34 #include <TSystem.h>
35 #include <TParticle.h>
36 #include <TH1I.h>
37 #include "TROOT.h"
38
39 #include "AliPID.h"
40 #include "AliTPCPIDResponse.h"
41 #include "AliStack.h"
42 #include "AliMCEvent.h"
43 #include "AliAnalysisManager.h"
44 #include "AliAODMCHeader.h"
45 #include "AliAODHandler.h"
46 #include "AliLog.h"
47 #include "AliAODVertex.h"
48 #include "AliAODJet.h"
49 #include "AliAODRecoDecay.h"
50 #include "AliAODRecoDecayHF.h"
51 #include "AliAODRecoCascadeHF.h"
52 #include "AliAODRecoDecayHF2Prong.h"
53 #include "AliAnalysisVertexingHF.h"
54 #include "AliESDtrack.h"
55 #include "AliAODMCParticle.h"
56 #include "AliAnalysisTaskSEDStarSpectra.h"
57
58 ClassImp(AliAnalysisTaskSEDStarSpectra)
59
60 //__________________________________________________________________________
61 AliAnalysisTaskSEDStarSpectra::AliAnalysisTaskSEDStarSpectra():
62   AliAnalysisTaskSE(),
63   fEvents(0),
64   fVHF(0),  
65   fMinITSClusters(0),
66   fMinITSClustersSoft(0),
67   fUseMCInfo(kTRUE), 
68   fOutput(0),
69   fNSigma(3),
70   fPID(kTRUE),
71   fAODTrack(0),
72   fMCDStarPt(0),    
73   fCEvents(0),     
74   fDStarMass(0),  
75   fTrueDiff(0),   
76   fTrueDiff2(0),  
77   fInvMass(0),    
78   fInvMass1(0),   
79   fInvMass2(0),  
80   fInvMass3(0),    
81   fInvMass4(0),   
82   fInvMass5(0),   
83   fPtDStar(0),     
84   fDStar(0),   
85   fDiff(0),   
86   fDiff1(0),    
87   fDiff2(0),    
88   fDiff3(0),  
89   fDiff4(0),  
90   fDiff5(0), 
91   fDiffSideBand(0),
92   fDiffSideBand1(0),
93   fDiffSideBand2(0),
94   fDiffSideBand3(0),
95   fDiffSideBand4(0),
96   fDiffSideBand5(0),
97   fDiffWrongSign(0),
98   fDiffWrongSign1(0),
99   fDiffWrongSign2(0),
100   fDiffWrongSign3(0),
101   fDiffWrongSign4(0),
102   fDiffWrongSign5(0)
103
104 {
105   //
106   // Default ctor
107   //
108 }
109 //___________________________________________________________________________
110 AliAnalysisTaskSEDStarSpectra::AliAnalysisTaskSEDStarSpectra(const Char_t* name) :
111   AliAnalysisTaskSE(name),
112   fEvents(0),
113   fVHF(0),  
114   fMinITSClusters(0),
115   fMinITSClustersSoft(0),
116   fUseMCInfo(kTRUE),
117   fOutput(0),
118   fNSigma(3),
119   fPID(kTRUE),
120   fAODTrack(0),
121   fMCDStarPt(0),    
122   fCEvents(0),     
123   fDStarMass(0),  
124   fTrueDiff(0),   
125   fTrueDiff2(0),  
126   fInvMass(0),    
127   fInvMass1(0),   
128   fInvMass2(0),  
129   fInvMass3(0),    
130   fInvMass4(0),   
131   fInvMass5(0),   
132   fPtDStar(0),    
133   fDStar(0),   
134   fDiff(0),   
135   fDiff1(0),    
136   fDiff2(0),    
137   fDiff3(0),  
138   fDiff4(0),  
139   fDiff5(0), 
140   fDiffSideBand(0),
141   fDiffSideBand1(0),
142   fDiffSideBand2(0),
143   fDiffSideBand3(0),
144   fDiffSideBand4(0),
145   fDiffSideBand5(0),
146   fDiffWrongSign(0),
147   fDiffWrongSign1(0),
148   fDiffWrongSign2(0),
149   fDiffWrongSign3(0),
150   fDiffWrongSign4(0),
151   fDiffWrongSign5(0)
152 {
153   //
154   // Constructor. Initialization of Inputs and Outputs
155   //
156   Info("AliAnalysisTaskSEDStarSpectra","Calling Constructor");
157   
158   DefineOutput(1,TList::Class());
159 }
160
161 //___________________________________________________________________________
162 AliAnalysisTaskSEDStarSpectra& AliAnalysisTaskSEDStarSpectra::operator=(const AliAnalysisTaskSEDStarSpectra& c) 
163 {
164   //
165   // Assignment operator
166   //
167   if (this!=&c) {
168     AliAnalysisTaskSE::operator=(c) ;
169   }
170   return *this;
171 }
172
173 //___________________________________________________________________________
174 AliAnalysisTaskSEDStarSpectra::AliAnalysisTaskSEDStarSpectra(const AliAnalysisTaskSEDStarSpectra& c) :
175   AliAnalysisTaskSE(c),
176   fEvents(c.fEvents),
177   fVHF(c.fVHF),  
178   fMinITSClusters(c.fMinITSClusters),
179   fMinITSClustersSoft(c.fMinITSClustersSoft),
180   fUseMCInfo(c.fUseMCInfo),
181   fOutput(c.fOutput),
182   fNSigma(c.fNSigma),
183   fPID(c.fPID),
184   fAODTrack(c.fAODTrack),
185   fMCDStarPt(c.fMCDStarPt),    
186   fCEvents(c.fCEvents),     
187   fDStarMass(c.fDStarMass),  
188   fTrueDiff(c.fTrueDiff),   
189   fTrueDiff2(c.fTrueDiff2),  
190   fInvMass(c.fInvMass),    
191   fInvMass1(c.fInvMass1),   
192   fInvMass2(c.fInvMass2),  
193   fInvMass3(c.fInvMass3),    
194   fInvMass4(c.fInvMass4),   
195   fInvMass5(c.fInvMass5),   
196   fPtDStar(c.fPtDStar),    
197   fDStar(c.fDStar),   
198   fDiff(c.fDiff),   
199   fDiff1(c.fDiff1),    
200   fDiff2(c.fDiff2),    
201   fDiff3(c.fDiff3),  
202   fDiff4(c.fDiff4),  
203   fDiff5(c.fDiff5), 
204   fDiffSideBand(c.fDiffSideBand),
205   fDiffSideBand1(c.fDiffSideBand1),
206   fDiffSideBand2(c.fDiffSideBand2),
207   fDiffSideBand3(c.fDiffSideBand3),
208   fDiffSideBand4(c.fDiffSideBand4),
209   fDiffSideBand5(c.fDiffSideBand5),
210   fDiffWrongSign(c.fDiffWrongSign),
211   fDiffWrongSign1(c.fDiffWrongSign1),
212   fDiffWrongSign2(c.fDiffWrongSign2),
213   fDiffWrongSign3(c.fDiffWrongSign3),
214   fDiffWrongSign4(c.fDiffWrongSign4),
215   fDiffWrongSign5(c.fDiffWrongSign5)
216 {
217   //
218   // Copy Constructor
219   //
220 }
221
222 //___________________________________________________________________________
223 AliAnalysisTaskSEDStarSpectra::~AliAnalysisTaskSEDStarSpectra() {
224   //
225   // destructor
226   //
227   Info("~AliAnalysisTaskSEDStarSpectra","Calling Destructor");
228   
229   if (fOutput) {
230     delete fOutput;
231     fOutput = 0;
232   }
233   if (fVHF) {
234     delete fVHF;
235     fVHF = 0;
236   } 
237 }
238 //_________________________________________________
239 void AliAnalysisTaskSEDStarSpectra::Init(){
240   //
241   // Initialization
242   //
243
244   if(fDebug > 1) printf("AnalysisTaskSEDStarSpectra::Init() \n");
245
246   gROOT->LoadMacro("$ALICE_ROOT/PWG3/vertexingHF/ConfigVertexingHF.C");
247   fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
248   //fVHF->PrintStatus();
249
250   return;
251 }
252
253 //_________________________________________________
254 void AliAnalysisTaskSEDStarSpectra::UserExec(Option_t *)
255 {
256   // user exec
257   if (!fInputEvent) {
258     Error("UserExec","NO EVENT FOUND!");
259     return;
260   }
261   
262   fCEvents->Fill(1);
263   // Load the event
264   fEvents++;
265   AliInfo(Form("Event %d",fEvents));
266   if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents));
267   AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
268   TClonesArray *arrayDStartoD0pi=0;
269   Init();
270   if(!aodEvent && AODEvent() && IsStandardAOD()) {
271     // In case there is an AOD handler writing a standard AOD, use the AOD 
272     // event in memory rather than the input (ESD) event.    
273     aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
274     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
275     // have to taken from the AOD event hold by the AliAODExtension
276     AliAODHandler* aodHandler = (AliAODHandler*) 
277       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
278     if(aodHandler->GetExtensions()) {
279       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
280       AliAODEvent *aodFromExt = ext->GetAOD();
281       arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
282     }
283   } else {
284     arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
285   }
286  
287  // AOD primary vertex
288   AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
289
290   // counters for efficiencies
291   Int_t icountReco = 0;
292   
293   //D* and D0 prongs needed to MatchToMC method
294   Int_t pdgDgDStartoD0pi[2]={421,211};
295   Int_t pdgDgD0toKpi[2]={321,211};
296   
297   if (!arrayDStartoD0pi){
298     AliInfo("Could not find array of HF vertices, skipping the event");
299     return;
300   }else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast())); 
301     
302   // loop over the tracks to search for candidates soft pion
303   
304   for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
305     
306     // D* candidates
307     AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
308     
309     // D0 from the reco cascade
310     AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
311     Bool_t unsetvtx=kFALSE;
312     
313     Double_t finvM =0;
314     Double_t finvMDStar = 0;
315     // needed for pointing angle
316     if(!theD0particle->GetOwnPrimaryVtx()) {
317       theD0particle->SetOwnPrimaryVtx(vtx1);
318       unsetvtx=kTRUE;
319     }    
320     Bool_t isDStar = 0;
321
322     // mc analysis 
323     if(fUseMCInfo){
324     //MC array need for maching
325       TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
326       if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
327       // find associated MC particle for D* ->D0toKpi
328       Int_t mcLabel = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
329       if(mcLabel>=0) isDStar = 1;
330     }
331
332     // soft pion
333     AliAODTrack *track2 = (AliAODTrack*)dstarD0pi->GetBachelor(); 
334
335     //D0tokpi
336     AliAODTrack *track0 = (AliAODTrack*)theD0particle->GetDaughter(0);
337     AliAODTrack *track1 = (AliAODTrack*)theD0particle->GetDaughter(1);
338     
339     Double_t pt = dstarD0pi->Pt();
340     
341     //kaon PID on low pt D*
342     if(pt<3 && fPID){
343       if (dstarD0pi->Charge()>0){
344         if(!SelectPID(track1,fNSigma)) continue;
345       }else{
346         if(!SelectPID(track0,fNSigma)) continue;
347       }
348     }
349
350     // reft in ITS for soft pion
351     //if((!(track2->GetStatus()&AliESDtrack::kITSrefit))) continue;
352        
353     // cut in acceptance for the soft pion and for the D0 daughters      
354     Bool_t acceptanceProng0 = (TMath::Abs(theD0particle->EtaProng(0))<= 0.9 && theD0particle->PtProng(0) >= 0.1);
355     Bool_t acceptanceProng1 = (TMath::Abs(theD0particle->EtaProng(1))<= 0.9 && theD0particle->PtProng(1) >= 0.1);
356     // soft pion acceptance ... is it fine 0.9?????
357     Bool_t acceptanceProng2 = (TMath::Abs(track2->Eta())<= 1.0 && track2->Pt() >= 0.05);
358     
359     if (acceptanceProng0 && acceptanceProng1 && acceptanceProng2) {
360       AliDebug(2,"D* reco daughters in acceptance");
361
362       // cut on the min n. of clusters in ITS for the D0 and soft pion
363       Int_t ncls0=0,ncls1=0,ncls2=0;
364       for(Int_t l=0;l<6;l++) {
365         if(TESTBIT(track0->GetITSClusterMap(),l)) ncls0++;
366         if(TESTBIT(track1->GetITSClusterMap(),l)) ncls1++;
367         if(TESTBIT(track2->GetITSClusterMap(),l)) ncls2++;
368       }       
369       // see AddTask for soft pion and D0 prongs ITS clusters request
370       if (ncls0 >= fMinITSClusters && ncls1 >= fMinITSClusters && ncls2>=fMinITSClustersSoft) { 
371  
372         // tag the D* decay do not double count background
373         Int_t decayTag = track2->Charge(); 
374         // D0 pt needed for the cuts       
375         Double_t ptD0 = theD0particle->Pt();
376        
377         Int_t okD0 = 0;
378         Int_t okD0bar = 0;
379         Int_t okD0WrongSign =0;
380         Int_t okD0barWrongSign =0;
381
382         // optimized D* cuts from UU
383         Bool_t tCutDStar = SetUtrechtSelections(ptD0,fVHF);  
384         Bool_t tCutOk = kFALSE;
385         if(tCutDStar) tCutOk = theD0particle->SelectD0(fVHF->GetD0toKpiCuts(),okD0,okD0bar);
386
387         // flags for wrong sign
388         okD0WrongSign=okD0;
389         okD0barWrongSign=okD0bar;
390
391         // search for D*
392         if(tCutOk){     
393           // correct decay 
394           if(decayTag>0 ? (okD0bar = 0) : (okD0 = 0));
395
396           finvM = dstarD0pi->InvMassD0();
397
398           if(okD0 == 1) fInvMass->Fill(finvM); 
399           if(okD0bar == 1) fInvMass->Fill(finvM); 
400
401           //DStar invariant mass
402           finvMDStar = dstarD0pi->InvMassDstarKpipi();
403           
404           if(finvM >= 1.829 && finvM <= 1.901){ // ~3 sigma cut on D0 mass
405            
406             if(isDStar == 1) { 
407               fDStarMass->Fill(finvMDStar); 
408               fTrueDiff ->Fill(dstarD0pi->DeltaInvMass());
409               fTrueDiff2->Fill(pt,dstarD0pi->DeltaInvMass());
410               fMCDStarPt->Fill(pt);  
411             }
412
413             // D* candidates - pt integrated
414             if(okD0==1){
415               fDStar->Fill(finvMDStar);
416               fDiff->Fill(dstarD0pi->DeltaInvMass()); // M(Kpipi)-M(Kpi) 
417             }else if(okD0bar==1){
418               fDStar->Fill(finvMDStar);
419               fDiff->Fill(dstarD0pi->DeltaInvMass()); // M(Kpipi)-M(Kpi) 
420             }
421
422             //D* candidates - pt bins
423             if(pt>1 && pt<=2){ // [1-2]
424               if(okD0==1){
425                 fDiff1->Fill(dstarD0pi->DeltaInvMass());
426                 fInvMass1->Fill(finvM);
427               }else if(okD0bar==1){
428                 fDiff1->Fill(dstarD0pi->DeltaInvMass());
429                 fInvMass1->Fill(finvM);
430               }
431             }
432             if( pt>2 && pt<=3){  // [2-3]
433               if(okD0==1){
434                 fDiff2->Fill(dstarD0pi->DeltaInvMass());
435                 fInvMass2->Fill(finvM);
436               }else if(okD0bar==1){
437                 fDiff2->Fill(dstarD0pi->DeltaInvMass());
438                 fInvMass2->Fill(finvM);
439               }
440             }
441             if(pt>3  && pt<=5){ // [3-5]
442               if(okD0==1){
443                 fDiff3->Fill(dstarD0pi->DeltaInvMass());
444                 fInvMass3->Fill(finvM);
445               }else if(okD0bar==1){
446                 fDiff3->Fill(dstarD0pi->DeltaInvMass());
447                 fInvMass3->Fill(finvM);
448               }
449             }
450             if( pt>5 && pt<8){ // [5-8]
451               if(okD0==1){
452                 fDiff4->Fill(dstarD0pi->DeltaInvMass());
453                 fInvMass4->Fill(finvM);
454               }else if(okD0bar==1){
455                 fDiff4->Fill(dstarD0pi->DeltaInvMass());
456                 fInvMass4->Fill(finvM); 
457               }
458             }       
459             if(pt>=8){ // [>8]
460               if(okD0==1){
461                 fDiff5->Fill(dstarD0pi->DeltaInvMass());
462                 fInvMass5->Fill(finvM);
463               }else if(okD0bar==1){
464                 fDiff5->Fill(dstarD0pi->DeltaInvMass());
465                 fInvMass5->Fill(finvM);
466               }
467             }
468             
469             // D* pt distribution
470             if((dstarD0pi->DeltaInvMass())>=0.143920 && (dstarD0pi->DeltaInvMass())<=0.14702){
471               if(okD0==1) fPtDStar->Fill(pt);
472               if(okD0bar==1) fPtDStar->Fill(pt);
473             }
474  
475           }     
476           // evaluate side band background
477           SideBandBackground(finvM,finvMDStar,pt,okD0,okD0bar);
478           
479           finvM      = 0;      
480           finvMDStar = 0;          
481           
482         } // end cuts
483         
484         // wrong sign background
485         if(tCutOk){     
486           // correct decay 
487           if(decayTag>0 ? ( okD0WrongSign = 0) : ( okD0barWrongSign = 0));
488
489           //wrong D0 inv mass
490           if(okD0WrongSign==1){
491             finvM = theD0particle->InvMassD0();
492           }else if(okD0WrongSign==0){
493             finvM = theD0particle->InvMassD0bar();
494           }
495           // wrong D* inv mass    
496           Double_t px[3],py[3],pz[3];
497           UInt_t pdg[3]={321,211,211};
498           pdg[0] = (decayTag>0 ? 321 : 211); // positive daughter of D0
499           px[0] = theD0particle->PxProng(0);
500           py[0] = theD0particle->PyProng(0);
501           pz[0] = theD0particle->PzProng(0);
502           pdg[1] = (decayTag>0 ? 211 : 321); // negative daughter of D0
503           px[1] = theD0particle->PxProng(1);
504           py[1] = theD0particle->PyProng(1);
505           pz[1] = theD0particle->PzProng(1);
506           pdg[2] = 211; // soft pion
507           px[2] = track2->Px();
508           py[2] = track2->Py();
509           pz[2] = track2->Pz();
510
511           Short_t dummycharge=0;
512           Double_t dummyd0[3]={0,0,0};
513           AliAODRecoDecay *rd = new AliAODRecoDecay(0x0,3,dummycharge,px,py,pz,dummyd0);
514           
515           finvMDStar = rd->InvMass(3,pdg);
516           
517           delete rd; rd=NULL;
518           
519           if(finvM >= 1.829 && finvM <= 1.901){ // ~3 sigma cut on D0 mass
520             WrongSignForDStar(finvM,finvMDStar,pt,okD0WrongSign,okD0barWrongSign);
521           }
522         }// end of wrong sign
523       }
524     }
525   }
526   
527   AliDebug(2, Form("Found %i Reco particles that are D*!!",icountReco));
528   
529   PostData(1,fOutput);
530 }
531 //________________________________________ terminate ___________________________
532 void AliAnalysisTaskSEDStarSpectra::Terminate(Option_t*)
533 {    
534   // The Terminate() function is the last function to be called during
535   // a query. It always runs on the client, it can be used to present
536   // the results graphically or save the results to file.
537   
538   Info("Terminate","");
539   AliAnalysisTaskSE::Terminate();
540   
541   fOutput = dynamic_cast<TList*> (GetOutputData(1));
542   if (!fOutput) {     
543     printf("ERROR: fOutput not available\n");
544     return;
545   }
546   
547   fMCDStarPt      = dynamic_cast<TH1F*>(fOutput->FindObject("fMCDStarPt"));
548   fCEvents        = dynamic_cast<TH1F*>(fOutput->FindObject("fCEvents"));
549   fDStarMass      = dynamic_cast<TH1F*>(fOutput->FindObject("fDStarMass"));
550   fTrueDiff       = dynamic_cast<TH1F*>(fOutput->FindObject("fTrueDiff"));
551   fTrueDiff2      = dynamic_cast<TH2F*>(fOutput->FindObject("fTrueDiff2"));
552   fInvMass        = dynamic_cast<TH1F*>(fOutput->FindObject("fInvMass"));
553   fInvMass1       = dynamic_cast<TH1F*>(fOutput->FindObject("fInvMass1"));
554   fInvMass2       = dynamic_cast<TH1F*>(fOutput->FindObject("fInvMass2"));
555   fInvMass3       = dynamic_cast<TH1F*>(fOutput->FindObject("fInvMass3"));
556   fInvMass4       = dynamic_cast<TH1F*>(fOutput->FindObject("fInvMass4"));
557   fInvMass5       = dynamic_cast<TH1F*>(fOutput->FindObject("fInvMass5"));
558   fPtDStar        = dynamic_cast<TH1F*>(fOutput->FindObject("fPtDStar "));
559   fDStar          = dynamic_cast<TH1F*>(fOutput->FindObject("fDStar"));
560   fDiff           = dynamic_cast<TH1F*>(fOutput->FindObject("fDiff"));
561   fDiff1          = dynamic_cast<TH1F*>(fOutput->FindObject("fDiff1"));
562   fDiff2          = dynamic_cast<TH1F*>(fOutput->FindObject("fDiff2"));
563   fDiff3          = dynamic_cast<TH1F*>(fOutput->FindObject("fDiff3"));
564   fDiff4          = dynamic_cast<TH1F*>(fOutput->FindObject("fDiff4"));
565   fDiff5          = dynamic_cast<TH1F*>(fOutput->FindObject("fDiff5"));
566   fDiffSideBand   = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffSideBand"));
567   fDiffSideBand1  = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffSideBand1"));
568   fDiffSideBand2  = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffSideBand2"));
569   fDiffSideBand3  = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffSideBand3"));
570   fDiffSideBand4  = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffSideBand4"));
571   fDiffSideBand5  = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffSideBand5"));
572   fDiffWrongSign  = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffWrongSign"));
573   fDiffWrongSign1 = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffWrongSign1"));
574   fDiffWrongSign2 = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffWrongSign2"));
575   fDiffWrongSign3 = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffWrongSign3"));
576   fDiffWrongSign4 = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffWrongSign4"));
577   fDiffWrongSign5 = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffWrongSign5"));
578
579 }
580 //___________________________________________________________________________
581 void AliAnalysisTaskSEDStarSpectra::UserCreateOutputObjects() { 
582  // output
583   Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
584   
585   //slot #1  
586   OpenFile(1);
587   fOutput = new TList();
588   fOutput->SetOwner();
589   // define histograms
590   DefineHistoFroAnalysis();
591   
592   return;
593 }
594
595 //___________________________________ hiostograms _______________________________________
596 Bool_t  AliAnalysisTaskSEDStarSpectra::DefineHistoFroAnalysis(){
597
598   // Invariant mass related histograms
599   fInvMass = new TH1F("invMass","K#pi invariant mass distribution",1500,.5,3.5);
600   fInvMass->SetStats(kTRUE);
601   fInvMass->GetXaxis()->SetTitle("M(K#pi) GeV/c");
602   fInvMass->GetYaxis()->SetTitle("Entries");
603   
604   fOutput->Add(fInvMass);
605   
606   fInvMass1 = new TH1F("invMass1","K#pi invariant mass distribution",3500,1.5,2.2);
607   fInvMass1->SetStats(kTRUE);
608   fInvMass1->GetXaxis()->SetTitle("M(K#pi) GeV/c");
609   fInvMass1->GetYaxis()->SetTitle("Entries");
610   
611   fOutput->Add(fInvMass1);
612   
613   fInvMass2 = new TH1F("invMass2","K#pi invariant mass distribution",350,1.5,2.2);
614   fInvMass2->SetStats(kTRUE);
615   fInvMass2->GetXaxis()->SetTitle("M(K#pi) GeV/c");
616   fInvMass2->GetYaxis()->SetTitle("Entries");
617   
618   fOutput->Add(fInvMass2);
619   
620   fInvMass3 = new TH1F("invMass3","K#pi invariant mass distribution",350,1.5,2.2);
621   fInvMass3->SetStats(kTRUE);
622   fInvMass3->GetXaxis()->SetTitle("M(K#pi) GeV/c");
623   fInvMass3->GetYaxis()->SetTitle("Entries");
624   
625   fOutput->Add(fInvMass3);
626   
627   fInvMass4 = new TH1F("invMass4","K#pi invariant mass distribution",350,1.5,2.2);
628   fInvMass4->SetStats(kTRUE);
629   fInvMass4->GetXaxis()->SetTitle("M(K#pi) GeV/c");
630   fInvMass4->GetYaxis()->SetTitle("Entries");
631
632   fOutput->Add(fInvMass4);
633   
634   fInvMass5 = new TH1F("invMass5","K#pi invariant mass distribution",350,1.5,2.2);
635   fInvMass5->SetStats(kTRUE);
636   fInvMass5->GetXaxis()->SetTitle("M(K#pi) GeV/c");
637   fInvMass5->GetYaxis()->SetTitle("Entries");
638
639   fOutput->Add(fInvMass5);
640
641   fDStar = new TH1F("invMassDStar","DStar invariant mass after D0 cuts ",600,1.8,2.4);
642   fDStar->SetStats(kTRUE);
643   fDStar->GetXaxis()->SetTitle("GeV/c");
644   fDStar->GetYaxis()->SetTitle("Entries");
645
646   fOutput->Add(fDStar);
647
648   fCEvents = new TH1F("fCEvents","conter",10,0,10);
649   fCEvents->SetStats(kTRUE);
650   fCEvents->GetXaxis()->SetTitle("1");
651   fCEvents->GetYaxis()->SetTitle("counts");
652
653   fOutput->Add(fCEvents);
654
655   fDiff = new TH1F("Diff","M(K#pi#pi)-M(K#pi)",750,0.1,0.2);
656   fDiff->SetStats(kTRUE);
657   fDiff->GetXaxis()->SetTitle("M(K#pi#pi)-M(K#pi) GeV/c^2");
658   fDiff->GetYaxis()->SetTitle("Entries");
659
660   fOutput->Add(fDiff);
661
662   fDiff1 = new TH1F("Diff1","M(K#pi#pi)-M(K#pi)",750,0.1,0.2);
663   fDiff1->SetStats(kTRUE);
664   fDiff1->GetXaxis()->SetTitle("M(K#pi#pi)-M(K#pi) GeV/c^2");
665   fDiff1->GetYaxis()->SetTitle("Entries");
666
667   fOutput->Add(fDiff1);
668
669   fDiff2 = new TH1F("Diff2","M(K#pi#pi)-M(K#pi)",750,0.1,0.2);
670   fDiff2->SetStats(kTRUE);
671   fDiff2->GetXaxis()->SetTitle("M(K#pi#pi)-M(K#pi) GeV/c^2");
672   fDiff2->GetYaxis()->SetTitle("Entries");
673
674   fOutput->Add(fDiff2);
675
676   fDiff3 = new TH1F("Diff3","M(K#pi#pi)-M(K#pi)",750,0.1,0.2);
677   fDiff3->SetStats(kTRUE);
678   fDiff3->GetXaxis()->SetTitle("M(K#pi#pi)-M(K#pi) GeV/c^2");
679   fDiff3->GetYaxis()->SetTitle("Entries");
680
681   fOutput->Add(fDiff3);
682
683   fDiff4 = new TH1F("Diff4","M(K#pi#pi)-M(K#pi)",750,0.1,0.2);
684   fDiff4->SetStats(kTRUE);
685   fDiff4->GetXaxis()->SetTitle("M(K#pi#pi)-M(K#pi) GeV/c^2");
686   fDiff4->GetYaxis()->SetTitle("Entries");
687
688   fOutput->Add(fDiff4);
689   
690   fDiff5 = new TH1F("Diff5","M(K#pi#pi)-M(K#pi)",750,0.1,0.2);
691   fDiff5->SetStats(kTRUE);
692   fDiff5->GetXaxis()->SetTitle("M(K#pi#pi)-M(K#pi) GeV/c^2");
693   fDiff5->GetYaxis()->SetTitle("Entries");
694
695   fOutput->Add(fDiff5);
696
697   fDiffSideBand = new TH1F("DiffSide","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
698   fDiffSideBand->SetStats(kTRUE);
699   fDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
700   fDiffSideBand->GetYaxis()->SetTitle("Entries");
701
702   fOutput->Add(fDiffSideBand); 
703
704   fDiffSideBand1 = new TH1F("DiffSide1","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
705   fDiffSideBand1->SetStats(kTRUE);
706   fDiffSideBand1->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
707   fDiffSideBand1->GetYaxis()->SetTitle("Entries");
708
709   fOutput->Add(fDiffSideBand1); 
710
711   fDiffSideBand2 = new TH1F("DiffSide2","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
712   fDiffSideBand2->SetStats(kTRUE);
713   fDiffSideBand2->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
714   fDiffSideBand2->GetYaxis()->SetTitle("Entries");
715
716   fOutput->Add(fDiffSideBand2); 
717
718   fDiffSideBand3 = new TH1F("DiffSide3","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
719   fDiffSideBand3->SetStats(kTRUE);
720   fDiffSideBand3->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
721   fDiffSideBand3->GetYaxis()->SetTitle("Entries");
722
723   fOutput->Add(fDiffSideBand3); 
724
725   fDiffSideBand4 = new TH1F("DiffSide4","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
726   fDiffSideBand4->SetStats(kTRUE);
727   fDiffSideBand4->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
728   fDiffSideBand4->GetYaxis()->SetTitle("Entries");
729
730   fOutput->Add(fDiffSideBand4); 
731
732   fDiffSideBand5 = new TH1F("DiffSide5","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
733   fDiffSideBand5->SetStats(kTRUE);
734   fDiffSideBand5->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
735   fDiffSideBand5->GetYaxis()->SetTitle("Entries");
736
737   fOutput->Add(fDiffSideBand5); 
738
739   fDiffWrongSign = new TH1F("DiffWrong","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
740   fDiffWrongSign->SetStats(kTRUE);
741   fDiffWrongSign->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
742   fDiffWrongSign->GetYaxis()->SetTitle("Entries");
743
744   fOutput->Add(fDiffWrongSign); 
745
746   fDiffWrongSign1 = new TH1F("DiffWrong1","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
747   fDiffWrongSign1->SetStats(kTRUE);
748   fDiffWrongSign1->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
749   fDiffWrongSign1->GetYaxis()->SetTitle("Entries");
750
751   fOutput->Add(fDiffWrongSign1); 
752
753   fDiffWrongSign2 = new TH1F("DiffWrong2","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
754   fDiffWrongSign2->SetStats(kTRUE);
755   fDiffWrongSign2->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
756   fDiffWrongSign2->GetYaxis()->SetTitle("Entries");
757
758   fOutput->Add(fDiffWrongSign2); 
759
760   fDiffWrongSign3 = new TH1F("DiffWrong3","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
761   fDiffWrongSign3->SetStats(kTRUE);
762   fDiffWrongSign3->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
763   fDiffWrongSign3->GetYaxis()->SetTitle("Entries");
764
765   fOutput->Add(fDiffWrongSign3); 
766
767   fDiffWrongSign4 = new TH1F("DiffWrong4","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
768   fDiffWrongSign4->SetStats(kTRUE);
769   fDiffWrongSign4->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
770   fDiffWrongSign4->GetYaxis()->SetTitle("Entries");
771
772   fOutput->Add(fDiffWrongSign4); 
773
774   fDiffWrongSign5 = new TH1F("DiffWrong5","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
775   fDiffWrongSign5->SetStats(kTRUE);
776   fDiffWrongSign5->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
777   fDiffWrongSign5->GetYaxis()->SetTitle("Entries");
778
779   fOutput->Add(fDiffWrongSign5); 
780  
781   fDStarMass = new TH1F("RECODStar2","RECO DStar invariant mass distribution",750,1.5,2.5);
782
783   fOutput->Add(fDStarMass);
784  
785   fTrueDiff  = new TH1F("dstar","True Reco diff",750,0,0.2);
786
787   fOutput->Add(fTrueDiff);
788   fTrueDiff2 = new TH2F("DiffDstar_pt","True Reco diff vs pt",100,0,15,900,0,0.3);
789
790   fOutput->Add(fTrueDiff2);
791
792   fPtDStar = new TH1F("DStarpt","Reconstructed D* candidates momentum ",500,0,25);
793   fPtDStar->SetStats(kTRUE);
794   fPtDStar->GetXaxis()->SetTitle("GeV/c");
795   fPtDStar->GetYaxis()->SetTitle("Entries");
796
797   fOutput->Add(fPtDStar);
798
799   fMCDStarPt = new TH1F("DStarMC2","Reconstructed D* momentum MC tagged ",500,0,25);
800   fMCDStarPt->SetStats(kTRUE);
801   fMCDStarPt->GetXaxis()->SetTitle("GeV/c");
802   fMCDStarPt->GetYaxis()->SetTitle("Entries");
803
804   fOutput->Add(fMCDStarPt);
805
806   return kTRUE;
807   
808 }
809 //______________________________ side band background for D*___________________________________
810 void AliAnalysisTaskSEDStarSpectra::SideBandBackground(Double_t finvM, Double_t finvMDStar,  Double_t pt, Int_t okD0, Int_t okD0bar ){
811
812   //  D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas 
813   // (expected detector resolution) on the left and right frm the D0 mass. Each band
814   //  has a width of ~5 sigmas. Two band needed  for opening angle considerations   
815   
816   if((finvM>=1.740 && finvM<=1.829) || (finvM>=1.901 && finvM<=1.990)){
817     
818     if(okD0==1)       fDiffSideBand->Fill(finvMDStar-finvM); // M(Kpipi)-M(Kpi) side band background
819     if(okD0bar==1)    fDiffSideBand->Fill(finvMDStar-finvM); // M(Kpipi)-M(Kpi) side band background
820
821     // pt bins
822     if( pt <= 1){
823       if(okD0==1)     fDiffSideBand1->Fill(finvMDStar-finvM);
824       if(okD0bar==1)  fDiffSideBand1->Fill(finvMDStar-finvM);
825     }
826     if( pt > 1 && pt <=3){
827       if(okD0==1 )    fDiffSideBand2->Fill(finvMDStar-finvM);
828       if(okD0bar==1 ) fDiffSideBand2->Fill(finvMDStar-finvM);
829     }
830     if( pt >3  && pt <=5){
831       if(okD0==1)     fDiffSideBand3->Fill(finvMDStar-finvM);
832       if(okD0bar==1)  fDiffSideBand3->Fill(finvMDStar-finvM);
833     }
834     if( pt > 5  && pt <8){
835       if(okD0==1)     fDiffSideBand4->Fill(finvMDStar-finvM);
836       if(okD0bar==1)  fDiffSideBand4->Fill(finvMDStar-finvM);
837     }  
838     if( pt >= 8){
839       if(okD0==1)     fDiffSideBand5->Fill(finvMDStar-finvM);
840       if(okD0bar==1)  fDiffSideBand5->Fill(finvMDStar-finvM);
841     }
842     
843   }
844 }
845
846 //________________________________________________________________________________________________________________
847 void AliAnalysisTaskSEDStarSpectra::WrongSignForDStar(Double_t finvM, Double_t finvMDStar,  Double_t pt, Int_t okD0, Int_t okD0bar){
848   //
849   // assign the wrong charge to the soft pion to create background
850   //
851   if(okD0==1 )       fDiffWrongSign->Fill(finvMDStar-finvM); // M(Kpipi)-M(Kpi) side band background
852   if(okD0bar==1 )     fDiffWrongSign->Fill(finvMDStar-finvM); // M(Kpipi)-M(Kpi) side band background
853   
854   // pt bins
855   if( pt <= 1){
856     if(okD0==1)      fDiffWrongSign1->Fill(finvMDStar-finvM);
857     if(okD0bar==1)   fDiffWrongSign1->Fill(finvMDStar-finvM);
858   }
859   if( pt > 1 && pt <=3){
860     if(okD0==1)      fDiffWrongSign2->Fill(finvMDStar-finvM);
861     if(okD0bar==1)   fDiffWrongSign2->Fill(finvMDStar-finvM);
862   }
863   if( pt >3  && pt <=5){
864     if(okD0==1)      fDiffWrongSign3->Fill(finvMDStar-finvM);
865     if(okD0bar==1)   fDiffWrongSign3->Fill(finvMDStar-finvM);
866   }
867   if( pt > 5  && pt <8){
868     if(okD0==1)      fDiffWrongSign4->Fill(finvMDStar-finvM);
869     if(okD0bar==1)   fDiffWrongSign4->Fill(finvMDStar-finvM);
870   }  
871   if( pt >= 8){
872     if(okD0==1)      fDiffWrongSign5->Fill(finvMDStar-finvM);
873     if(okD0bar==1)   fDiffWrongSign5->Fill(finvMDStar-finvM);
874   }
875   
876 }
877 //_____________________________________________________________________________________________
878 Bool_t AliAnalysisTaskSEDStarSpectra::SetUtrechtSelections(Double_t ptD0, AliAnalysisVertexingHF *fVHF){
879
880   //cuts[0] = inv. mass half width [GeV]
881   //cuts[1] = dca [cm]
882   //cuts[2] = cosThetaStar
883   //cuts[3] = pTK [GeV/c]
884   //cuts[4] = pTPi [GeV/c]
885   //cuts[5] = d0K [cm]   upper limit!
886   //cuts[6] = d0Pi [cm]  upper limit!
887   //cuts[7] = d0d0 [cm^2]
888   //cuts[8] = cosThetaPoint
889
890   if(ptD0<1.){
891     fVHF->SetD0toKpiCuts(0.450,0.04,0.8,0.21,0.21,0.021,0.021,-0.0002,0.9);
892   }
893   else if(ptD0 >=1. && ptD0 <=2.){
894     fVHF->SetD0toKpiCuts(0.450,0.02,0.7,0.8,0.8,0.021,0.021,-0.0002,0.9);
895   }  
896   else if(ptD0 >2. && ptD0 <=3.){
897     fVHF->SetD0toKpiCuts(0.450,0.04,0.8,0.8,0.8,0.035,0.042,-0.000085,0.9);
898   } 
899   else if(ptD0 >3. && ptD0 <=5.){
900     fVHF->SetD0toKpiCuts(0.450,0.016,0.8,1.2,1.2,0.042,0.056,-0.000085,0.9);
901   } 
902   else if(ptD0 >5.){
903     fVHF->SetD0toKpiCuts(0.450,0.08,1.0,1.2,1.2,0.07,0.07,0.0001,0.9);
904   }  
905   return kTRUE;
906 }
907 //_____________________________ pid _______________________________________-
908 Bool_t AliAnalysisTaskSEDStarSpectra::SelectPID(AliAODTrack *track, Double_t nsig){//pid - K = 3
909   //TPC
910   Bool_t isKaon=kTRUE;
911   if ((track->GetStatus()&AliESDtrack::kTPCpid )==0) return isKaon;
912   AliAODPid *pid = track->GetDetPid();
913   static AliTPCPIDResponse theTPCpid;
914   Double_t nsigma = theTPCpid.GetNumberOfSigmas(track->P(),pid->GetTPCsignal(),track->GetTPCClusterMap().CountBits(), (AliPID::kKaon));
915   if (TMath::Abs(nsigma)>nsig) isKaon=kFALSE;
916   //ITS
917   //
918   //
919   //
920
921   return isKaon;
922 }
923