Update and addition of LS analysis (Renu, Giacomo, Francesco)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSEDplus.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2008, 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 //
18 // AliAnalysisTaskSE for the extraction of signal(e.g D+) of heavy flavor
19 // decay candidates with the MC truth.
20 // Authors: Renu Bala, bala@to.infn.it
21 // F. Prino, prino@to.infn.it
22 // G. Ortona, ortona@to.infn.it
23 /////////////////////////////////////////////////////////////
24
25 #include <TClonesArray.h>
26 #include <TNtuple.h>
27 #include <TList.h>
28 #include <TString.h>
29 #include <TH1F.h>
30 #include <TDatabasePDG.h>
31
32 #include "AliAnalysisManager.h"
33 #include "AliAODHandler.h"
34 #include "AliAODEvent.h"
35 #include "AliAODVertex.h"
36 #include "AliAODTrack.h"
37 #include "AliAODMCHeader.h"
38 #include "AliAODMCParticle.h"
39 #include "AliAODRecoDecayHF3Prong.h"
40 #include "AliAnalysisVertexingHF.h"
41 #include "AliAnalysisTaskSE.h"
42 #include "AliAnalysisTaskSEDplus.h"
43
44 ClassImp(AliAnalysisTaskSEDplus)
45
46
47 //________________________________________________________________________
48 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus():
49 AliAnalysisTaskSE(),
50 fOutput(0), 
51 fHistNEvents(0),
52 fNtupleDplus(0),
53 fHistLS(0),
54 fHistLStrip(0),
55 fHistLStripcut(0),
56 fHistOS(0),
57 fHistOSbkg(0),
58 fHistLSnoweight(0),
59 fHistLSsinglecut(0),
60 fFillNtuple(kFALSE),
61 fReadMC(kFALSE),
62 fDoLS(kFALSE),
63 fVHF(0)
64 {
65   // Default constructor
66 }
67
68 //________________________________________________________________________
69 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name,Bool_t fillNtuple):
70 AliAnalysisTaskSE(name),
71 fOutput(0),
72 fHistNEvents(0),
73 fNtupleDplus(0),
74 fHistLS(0),
75 fHistLStrip(0),
76 fHistLStripcut(0),
77 fHistOS(0),
78 fHistOSbkg(0),
79 fHistLSnoweight(0),
80 fHistLSsinglecut(0),
81 fFillNtuple(fillNtuple),
82 fReadMC(kFALSE),
83 fDoLS(kFALSE),
84 fVHF(0)
85 {
86   // Default constructor
87
88   // Output slot #1 writes into a TList container
89   DefineOutput(1,TList::Class());  //My private output
90
91   if(fFillNtuple){
92     // Output slot #2 writes into a TNtuple container
93     DefineOutput(2,TNtuple::Class());  //My private output
94   }
95 }
96
97 //________________________________________________________________________
98 AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus()
99 {
100   // Destructor
101   if (fOutput) {
102     delete fOutput;
103     fOutput = 0;
104   }
105   if(fNtupleDplus){
106     delete fNtupleDplus;
107     fNtupleDplus=0;
108   }
109   if (fVHF) {
110     delete fVHF;
111     fVHF = 0;
112   }
113
114 }  
115
116 //________________________________________________________________________
117 void AliAnalysisTaskSEDplus::Init()
118 {
119   // Initialization
120
121   if(fDebug > 1) printf("AnalysisTaskSEDplus::Init() \n");
122
123   gROOT->LoadMacro("ConfigVertexingHF.C");
124
125   fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");  
126   fVHF->PrintStatus();
127
128   return;
129 }
130
131 //________________________________________________________________________
132 void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
133 {
134   // Create the output container
135   //
136   if(fDebug > 1) printf("AnalysisTaskSEDplus::UserCreateOutputObjects() \n");
137
138   // Several histograms are more conveniently managed in a TList
139   fOutput = new TList();
140   fOutput->SetOwner();
141   fOutput->SetName("OutputHistos");
142
143   Int_t nPtBins=4;
144   TString hisname;
145   for(Int_t i=0;i<nPtBins;i++){
146     hisname.Form("hMassPt%d",i);
147     TH1F* hm=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
148     hm->Sumw2();
149     hisname.Form("hSigPt%d",i);
150     TH1F* hs=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
151     hs->Sumw2();
152     hisname.Form("hBkgPt%d",i);
153     TH1F* hb=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
154     hb->Sumw2();
155
156     hisname.Form("hMassPt%dTC",i);
157     TH1F* hmtc=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
158     hmtc->Sumw2();
159     hisname.Form("hSigPt%dTC",i);
160     TH1F* hstc=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
161     hstc->Sumw2();
162     hisname.Form("hBkgPt%dTC",i);
163     TH1F* hbtc=new TH1F(hisname.Data(),hisname.Data(),100,1.765,1.965);
164     hbtc->Sumw2();
165
166
167     fOutput->Add(hm);
168     fOutput->Add(hs);
169     fOutput->Add(hb);
170     fOutput->Add(hmtc);
171     fOutput->Add(hstc);
172     fOutput->Add(hbtc);
173   }
174   fHistNEvents = new TH1F("fHistNEvents", "Number of processed events; ; Events",3,-1.5,1.5);
175   fHistNEvents->Sumw2();
176   fHistNEvents->SetMinimum(0);
177   fOutput->Add(fHistNEvents);
178   if(fDoLS){
179     Double_t massDplus=TDatabasePDG::Instance()->GetParticle(411)->Mass();
180   
181     fHistLS = new TH1F("fHistLS","fHistLS",100,massDplus-0.2,massDplus+0.2);
182     fHistLS->Sumw2();
183     fOutput->Add(fHistLS);
184
185     fHistLStrip = new TH1F("fHistLStrip","fHistLStrip",100,massDplus-0.2,massDplus+0.2);
186     fHistLStrip->Sumw2();
187     fOutput->Add(fHistLStrip);
188
189     fHistLStripcut = new TH1F("fHistLStripcut","fHistLStripcut",100,massDplus-0.2,massDplus+0.2);
190     fHistLStripcut->Sumw2();
191     fOutput->Add(fHistLStripcut);
192
193     fHistOS = new TH1F("fHistOS","fHistOS",100,massDplus-0.2,massDplus+0.2);
194     fHistOS->Sumw2();
195     fOutput->Add(fHistOS);
196
197     fHistOSbkg = new TH1F("fHistOSbkg","fHistOSbkg",100,massDplus-0.2,massDplus+0.2);
198     fHistOSbkg->Sumw2();
199     fOutput->Add(fHistOSbkg);
200
201     fHistLSnoweight = new TH1F("fHistLSnoweight","fHistLSnoweight",100,massDplus-0.2,massDplus+0.2);
202     fHistLSnoweight->Sumw2();
203     fOutput->Add(fHistLSnoweight);
204
205     fHistLSsinglecut = new TH1F("fHistLSsinglecut","fHistLSsinglecut",100,massDplus-0.2,massDplus+0.2);
206     fHistLSsinglecut->Sumw2();
207     fOutput->Add(fHistLSsinglecut);
208   }
209
210   if(fFillNtuple){
211     OpenFile(2); // 2 is the slot number of the ntuple
212     fNtupleDplus = new TNtuple("fNtupleDplus","D +","pdg:Px:Py:Pz:PtTrue:VxTrue:VyTrue:VzTrue:Ptpi:PtK:Ptpi2:PtRec:PointingAngle:DecLeng:VxRec:VyRec:VzRec:InvMass:sigvert:d0Pi:d0K:d0Pi2");
213   }
214
215   return;
216 }
217
218 //________________________________________________________________________
219 void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
220 {
221   // Execute analysis for current event:
222   // heavy flavor candidates association to MC truth
223
224   
225   
226   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
227   fHistNEvents->Fill(0); // count event
228   // Post the data already here
229   PostData(1,fOutput);
230   
231
232   TClonesArray *array3Prong = 0;
233   TClonesArray *arrayLikeSign =0;
234   if(!aod && AODEvent() && IsStandardAOD()) {
235     // In case there is an AOD handler writing a standard AOD, use the AOD 
236     // event in memory rather than the input (ESD) event.    
237     aod = dynamic_cast<AliAODEvent*> (AODEvent());
238     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
239     // have to taken from the AOD event hold by the AliAODExtension
240     AliAODHandler* aodHandler = (AliAODHandler*) 
241       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
242     if(aodHandler->GetExtensions()) {
243       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
244       AliAODEvent *aodFromExt = ext->GetAOD();
245       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
246       arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign3Prong");
247     }
248   } else {
249     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
250     arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign3Prong");
251   }
252
253   if(!array3Prong) {
254     printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
255     return;
256   }
257   if(!arrayLikeSign) {
258     printf("AliAnalysisTaskSEDplus::UserExec: LikeSign3Prong branch not found!\n");
259     // do in any case the OS analysis.
260   }
261
262  
263   TClonesArray *arrayMC=0;
264   AliAODMCHeader *mcHeader=0;
265
266   // AOD primary vertex
267   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
268   //    vtx1->Print();
269   
270   // load MC particles
271   if(fReadMC){
272     
273     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
274     if(!arrayMC) {
275       printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
276       return;
277     }
278     
279   // load MC header
280     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
281     if(!mcHeader) {
282       printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
283       return;
284     }
285   }
286   
287   Int_t n3Prong = array3Prong->GetEntriesFast();
288   if(fDebug>1) printf("Number of D+->Kpipi: %d\n",n3Prong);
289   
290   
291   Int_t nOS=0;
292   Int_t pdgDgDplustoKpipi[3]={321,211,211};
293   Double_t cutsDplus[12]={0.2,0.4,0.4,0.,0.,0.01,0.06,0.02,0.,0.85,0.,10000000000.};
294   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
295     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
296     
297     
298     Bool_t unsetvtx=kFALSE;
299     if(!d->GetOwnPrimaryVtx()){
300       d->SetOwnPrimaryVtx(vtx1);
301       unsetvtx=kTRUE;
302     }
303
304     if(d->SelectDplus(fVHF->GetDplusCuts())) {
305       Int_t iPtBin=0;
306       Double_t ptCand = d->Pt();
307       if(ptCand<2.){
308         iPtBin=0;
309         cutsDplus[7]=0.08;
310         cutsDplus[8]=0.5;
311         cutsDplus[9]=0.979;
312         cutsDplus[10]=0.0000055;
313       }
314       else if(ptCand>2. && ptCand<3){ 
315         iPtBin=1;
316         cutsDplus[7]=0.08;
317         cutsDplus[8]=0.5;
318         cutsDplus[9]=0.991;
319         cutsDplus[10]=0.000005;
320       }else if(ptCand>3. && ptCand<5){ 
321         iPtBin=2;
322         cutsDplus[7]=0.1;
323         cutsDplus[8]=0.5;
324         cutsDplus[9]=0.9955;
325         cutsDplus[10]=0.0000035;
326       }else{
327         iPtBin=3;
328         cutsDplus[7]=0.1;
329         cutsDplus[8]=0.5;
330         cutsDplus[9]=0.997;
331         cutsDplus[10]=0.000001;
332       }
333       Bool_t passTightCuts=d->SelectDplus(cutsDplus);
334       Int_t labDp=-1;
335       Float_t deltaPx=0.;
336       Float_t deltaPy=0.;
337       Float_t deltaPz=0.;
338       Float_t truePt=0.;
339       Float_t xDecay=0.;
340       Float_t yDecay=0.;
341       Float_t zDecay=0.;
342       Float_t pdgCode=-2;
343       if(fReadMC){
344         labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi);
345         if(labDp>=0){
346           AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
347           AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0));
348           deltaPx=partDp->Px()-d->Px();
349           deltaPy=partDp->Py()-d->Py();
350           deltaPz=partDp->Pz()-d->Pz();
351           truePt=partDp->Pt();
352           xDecay=dg0->Xv();       
353           yDecay=dg0->Yv();       
354           zDecay=dg0->Zv();
355           pdgCode=TMath::Abs(partDp->GetPdgCode());
356         }else{
357           pdgCode=-1;
358         }
359       }
360       Double_t invMass=d->InvMassDplus();
361
362       TString hisNameA(Form("hMassPt%d",iPtBin));
363       TString hisNameS(Form("hSigPt%d",iPtBin));
364       TString hisNameB(Form("hBkgPt%d",iPtBin));
365       TString hisNameATC(Form("hMassPt%dTC",iPtBin));
366       TString hisNameSTC(Form("hSigPt%dTC",iPtBin));
367       TString hisNameBTC(Form("hBkgPt%dTC",iPtBin));
368       
369       ((TH1F*)(fOutput->FindObject(hisNameA)))->Fill(invMass);
370       if(passTightCuts){
371         ((TH1F*)(fOutput->FindObject(hisNameATC)))->Fill(invMass);
372       }
373
374       Float_t tmp[22];
375       if(fFillNtuple){            
376         tmp[0]=pdgCode;
377         tmp[1]=deltaPx;
378         tmp[2]=deltaPy;
379         tmp[3]=deltaPz;
380         tmp[4]=truePt;
381         tmp[5]=xDecay;    
382         tmp[6]=yDecay;    
383         tmp[7]=zDecay;    
384         tmp[8]=d->PtProng(0);
385         tmp[9]=d->PtProng(1);
386         tmp[10]=d->PtProng(2);
387         tmp[11]=d->Pt();
388         tmp[12]=d->CosPointingAngle();
389         tmp[13]=d->DecayLength();
390         tmp[14]=d->Xv();
391         tmp[15]=d->Yv();
392         tmp[16]=d->Zv();
393         tmp[17]=d->InvMassDplus();
394         tmp[18]=d->GetSigmaVert();
395         tmp[19]=d->Getd0Prong(0);
396         tmp[20]=d->Getd0Prong(1);
397         tmp[21]=d->Getd0Prong(2);         
398         fNtupleDplus->Fill(tmp);
399         PostData(2,fNtupleDplus);
400       }
401
402       if(fReadMC){
403         if(labDp>=0) {
404           ((TH1F*)(fOutput->FindObject(hisNameS)))->Fill(invMass);
405           if(passTightCuts){
406             ((TH1F*)(fOutput->FindObject(hisNameSTC)))->Fill(invMass);
407           }
408           
409         }else{
410           ((TH1F*)(fOutput->FindObject(hisNameB)))->Fill(invMass);
411           if(passTightCuts){
412             ((TH1F*)(fOutput->FindObject(hisNameBTC)))->Fill(invMass);
413           }
414           fHistOSbkg->Fill(invMass);
415         }
416       }
417
418       fHistOS->Fill(invMass);
419       nOS++;
420     }
421     if(unsetvtx) d->UnsetOwnPrimaryVtx();
422   }
423  
424   //start LS analysis
425   if(fDoLS && arrayLikeSign) LSAnalysis(array3Prong,arrayLikeSign,aod,vtx1,nOS);
426    
427   PostData(1,fOutput);    
428   return;
429 }
430
431
432
433 //_________________________________________________________________
434 void AliAnalysisTaskSEDplus::LSAnalysis(TClonesArray *arrayOppositeSign,TClonesArray *arrayLikeSign,AliAODEvent *aod,AliAODVertex *vtx1, Int_t nDplusOS){
435
436 /*
437  * Fill the Like Sign histograms
438  */
439
440
441
442
443   //count pos/neg tracks
444   Int_t nPosTrks=0,nNegTrks=0;
445   //counter for particles passing single particle cuts
446   Int_t nspcplus=0;
447   Int_t nspcminus=0;
448
449   for(Int_t it=0;it<aod->GetNumberOfTracks();it++) {
450     AliAODTrack *track = aod->GetTrack(it);
451     if(track->Charge()>0){
452       nPosTrks++;
453       if(track->Pt()>=0.4){
454         nspcplus++;
455       }
456     }else{      
457       nNegTrks++;
458       if(track->Pt()>=0.4){
459         nspcminus++;
460       }
461     }
462   }
463
464   Int_t nOStriplets = arrayOppositeSign->GetEntriesFast();
465
466   Int_t nDplusLS=0;
467   Int_t nLikeSign = arrayLikeSign->GetEntriesFast();
468  
469   for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
470     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
471     Bool_t unsetvtx=kFALSE;
472     if(!d->GetOwnPrimaryVtx()) {
473       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
474       unsetvtx=kTRUE;
475     }
476     if(d->SelectDplus(fVHF->GetDplusCuts()))nDplusLS++;
477     if(unsetvtx) d->UnsetOwnPrimaryVtx();
478   }
479
480   Float_t wei2=0;
481   if(nLikeSign!=0)wei2 = (Float_t)nOStriplets/(Float_t)nLikeSign;
482   Float_t wei3=0;
483   if(nDplusLS!=0)wei3 = (Float_t)nDplusOS/(Float_t)nDplusLS;
484
485  // loop over  sign candidates
486   for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
487     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
488     Bool_t unsetvtx=kFALSE;
489     if(!d->GetOwnPrimaryVtx()) {
490       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
491       unsetvtx=kTRUE;
492     }
493  
494     if(d->SelectDplus(fVHF->GetDplusCuts())){
495       Int_t sign= d->GetCharge();
496       Float_t wei1=1;
497       Float_t wei4=1;
498       if(sign>0&&nPosTrks>2&&nspcplus>2){//wei* should be automatically protected, since to get a triplet there must be at least 3 good tracks in the event
499         wei1=3.*(Float_t)nNegTrks/((Float_t)nPosTrks-2.);
500         wei4=3.*(Float_t)nspcminus/((Float_t)nspcplus-2.);
501       }
502       if(sign<0&&nNegTrks>2&&nspcminus>2){
503         wei1=3.*(Float_t)nPosTrks/((Float_t)nNegTrks-2.);
504         wei4=3.*(Float_t)nspcplus/((Float_t)nspcminus-2.);
505       }
506       Double_t invMass=d->InvMassDplus();
507       fHistLS->Fill(invMass,wei1);
508       fHistLSnoweight->Fill(invMass);
509       fHistLStrip->Fill(invMass,wei2);
510       fHistLStripcut->Fill(invMass,wei3);
511       fHistLSsinglecut->Fill(invMass,wei4);
512      }
513     if(unsetvtx) d->UnsetOwnPrimaryVtx();
514   }
515   
516   //printf("------------ N. of positive tracks in Event ----- %d \n", nPosTrks);
517   //printf("------------ N. of negative tracks in Event ----- %d \n", nNegTrks);
518
519   //  printf("LS analysis...done\n");
520
521 }
522
523 //_________________________________________________________________
524
525
526 void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
527 {
528   // Terminate analysis
529   //
530   if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");
531   fOutput = dynamic_cast<TList*> (GetOutputData(1));
532   if (!fOutput) {     
533     printf("ERROR: fOutput not available\n");
534     return;
535   }
536   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
537   if(fFillNtuple){
538     fNtupleDplus = dynamic_cast<TNtuple*>(GetOutputData(2));
539   }
540   fHistOS =  dynamic_cast<TH1F*>(fOutput->FindObject("fHistOS"));
541   fHistOSbkg =  dynamic_cast<TH1F*>(fOutput->FindObject("fHistOSbkg"));
542   if(fDoLS){
543     fHistLS =  dynamic_cast<TH1F*>(fOutput->FindObject("fHistLS"));
544     fHistLStrip =  dynamic_cast<TH1F*>(fOutput->FindObject("fHistLStrip"));
545     fHistLStripcut =  dynamic_cast<TH1F*>(fOutput->FindObject("fHistLStripcut"));
546     fHistLSnoweight =  dynamic_cast<TH1F*>(fOutput->FindObject("fHistLSnoweight"));
547     fHistLSsinglecut =  dynamic_cast<TH1F*>(fOutput->FindObject("fHistLSsinglecut"));
548   }
549
550
551   return;
552 }