Fix
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / charmFlow / AliAnalysisTaskFlowD2H.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 // FlowD2H main task:
18 // >> Make flowEvent with RPcuts and POIcuts given in constructor and passes it
19 //    to the daughter tasks.
20 // >> The POIcuts are polymorphic based on the AliRDHFCuts class allowing the 
21 //    use of all charmed candidates reconstructed in the central barrel.
22 // Author: Carlos Perez (cperez@cern.ch)
23 //==============================================================================
24
25 /* $Id$ */
26
27 #include "TChain.h"
28 #include "TList.h"
29 #include "TH2D.h"
30 #include "TProfile.h"
31
32 #include "AliAnalysisTaskSE.h"
33 #include "AliAnalysisManager.h"
34
35 #include "AliCentrality.h"
36
37 #include "AliAODEvent.h"
38 #include "AliAODInputHandler.h"
39 #include "AliAODTrack.h"
40
41 #include "TMath.h"
42 #include "TObjArray.h"
43
44 #include "AliFlowCandidateTrack.h"
45 #include "AliFlowTrackCuts.h"
46 #include "AliFlowEvent.h"
47 #include "AliFlowCommonConstants.h"
48
49 #include "AliRDHFCuts.h"
50 #include "AliRDHFCutsD0toKpi.h"
51 #include "AliRDHFCutsDStartoKpipi.h"
52 #include "AliRDHFCutsDplustoKpipi.h"
53
54 #include "AliAODRecoDecayHF2Prong.h"
55 #include "AliAODRecoDecayHF3Prong.h"
56 #include "AliAODRecoCascadeHF.h"
57
58 #include "AliAnalysisTaskFlowD2H.h"
59
60 ClassImp(AliAnalysisTaskFlowD2H)
61
62 //_____________________________________________________________________________
63 AliAnalysisTaskFlowD2H::AliAnalysisTaskFlowD2H() :
64   AliAnalysisTaskSE(), fCutsRP(NULL), fCutsPOI(NULL), fSource(0),
65   fDebugV2(kFALSE), fHList(NULL), fAnaCuts(NULL)
66 {
67 // Default constructor
68   for(int i=0; i!=2; ++i)
69     fEvent[i] = NULL;
70   for(int i=0; i!=2; ++i)
71     fMass[i] = NULL;
72   for(int i=0; i!=2; ++i)
73     fPOIEta[i] = 0;
74   for(int i=0; i!=4; ++i)
75     fFlowEta[i] = 0;
76   for(int x=0; x!=2; ++x)
77     fFlowPts[x] = 0;
78   for(int x=0; x!=2; ++x)
79     for(int m=0; m!=5; ++m)
80       fFlowBands[x][m] = 0;
81 }
82
83 //_____________________________________________________________________________
84 AliAnalysisTaskFlowD2H::AliAnalysisTaskFlowD2H(const char *name,
85     AliFlowTrackCuts *cutsRPs, AliRDHFCuts *cutsPOIs, Int_t specie) :
86   AliAnalysisTaskSE(name), fCutsRP(cutsRPs), 
87   fCutsPOI(cutsPOIs), fSource(specie),
88   fDebugV2(kFALSE), fHList(NULL), fAnaCuts(NULL)
89 {
90 // Standard constructor
91   for(int i=0; i!=2; ++i)
92     fEvent[i] = NULL;
93   for(int i=0; i!=2; ++i)
94     fMass[i] = NULL;
95   for(int i=0; i!=2; ++i)
96     fPOIEta[i] = 0;
97   for(int i=0; i!=4; ++i)
98     fFlowEta[i] = 0;
99   for(int x=0; x!=2; ++x)
100     fFlowPts[x] = 0;
101   for(int x=0; x!=2; ++x)
102     for(int m=0; m!=5; ++m)
103       fFlowBands[x][m] = 0;
104
105   DefineInput( 0,TChain::Class());
106   DefineOutput(1,TList::Class());
107   DefineOutput(2,AliFlowEventSimple::Class()); // first band
108   DefineOutput(3,AliFlowEventSimple::Class()); // second band
109   DefineOutput(4,AliFlowEventSimple::Class()); // third band
110   DefineOutput(5,AliFlowEventSimple::Class()); // fourth band
111   DefineOutput(6,AliFlowEventSimple::Class()); // fifth band
112 }
113
114 //_____________________________________________________________________________
115 AliAnalysisTaskFlowD2H::~AliAnalysisTaskFlowD2H(){
116   // delete objects
117   if(fHList)
118     delete fHList;
119 }
120
121 //_____________________________________________________________________________
122 void AliAnalysisTaskFlowD2H::UserCreateOutputObjects(){
123   // Define output objects + parameters
124   if(fDebugV2)
125     printf("DEBUGGER: Creating output\n");
126   fHList = new TList();
127   fHList->SetOwner();
128   AddHistograms();
129   PostData(1,fHList);
130
131   AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
132   cc->SetNbinsMult(10000);
133   cc->SetMultMin(0);
134   cc->SetMultMax(10000);
135
136   cc->SetNbinsPt(fFlowPts[1]-fFlowPts[0]);
137   cc->SetPtMin(fFlowPts[0]);
138   cc->SetPtMax(fFlowPts[1]);
139
140   cc->SetNbinsPhi(180);
141   cc->SetPhiMin(0.0);
142   cc->SetPhiMax(TMath::TwoPi());
143
144   cc->SetNbinsEta(200);
145   cc->SetEtaMin(-5.0);
146   cc->SetEtaMax(+5.0);
147
148   cc->SetNbinsQ(500);
149   cc->SetQMin(0.0);
150   cc->SetQMax(3.0);
151 }
152 //_____________________________________________________________________________
153 void AliAnalysisTaskFlowD2H::AddHistograms()
154 {
155   // Create histograms
156   TList *tEvents = new TList();
157     tEvents->SetName("Events");
158     tEvents->SetOwner();
159     for(int i=0; i!=2; ++i) {
160       fEvent[i] = new TH2D(Form("Event%d",i),"Events;V0M;Arb",10,0,100,10,0,12000);
161       tEvents->Add(fEvent[i]);
162     }
163     fAnaCuts=new TProfile("Cuts","Analysis Cuts", 10,0,10);
164       fAnaCuts->Fill(0.5,fPOIEta[0],1);fAnaCuts->GetXaxis()->SetBinLabel(1,"ETAm");
165       fAnaCuts->Fill(1.5,fPOIEta[1],1);fAnaCuts->GetXaxis()->SetBinLabel(2,"ETAM");
166     tEvents->Add(fAnaCuts);
167     fHList->Add(tEvents);
168   TList *tCandidates;
169     tCandidates = new TList();
170     tCandidates->SetOwner();
171     tCandidates->SetName(Form("Candidates%d",fSource));
172     Double_t dBinMin=0, dBinMax=3;
173     Int_t nBins=3;
174     switch (fSource) {
175       case (AliRDHFCuts::kD0toKpiCuts): // 360/72 (bw=5MeV)
176         dBinMin = 1.695; dBinMax = 2.055; nBins=72; break;
177       case (AliRDHFCuts::kDstarCuts): // 36/72 (bw=0.5MeV)
178         dBinMin = 0.137; dBinMax = 0.173; nBins=72; break;
179       case (AliRDHFCuts::kDplusCuts): // 360/72 (bw=5MeV)
180         dBinMin = 1.695; dBinMax = 2.055; nBins=72; break;
181     }
182     for(int i=0; i!=2; ++i) {
183       fMass[i] = new TH2D( Form("Mass%d",i),
184                            Form("Mass%d;Mass [GeV];Pt [GeV]",i),
185                            nBins, dBinMin, dBinMax,
186                            fFlowPts[1]-fFlowPts[0], fFlowPts[0], fFlowPts[1] );
187       tCandidates->Add(fMass[i]);
188     }
189     fHList->Add(tCandidates);
190   if (fDebugV2) printf("DEBUGGER: Created histos for DMeson %d\n", fSource );
191 }
192 //_____________________________________________________________________________
193 void AliAnalysisTaskFlowD2H::NotifyRun()
194 {
195   // Do nothing
196 }
197 //_____________________________________________________________________________
198 void AliAnalysisTaskFlowD2H::UserExec(Option_t *)
199 {
200   // Do analysis + fIll histograms
201   AliAODEvent *fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
202
203   if(!fAOD)
204     return;
205   Int_t iMulti = fAOD->GetNTracks();  /// TO DO
206   fEvent[0]->Fill( fCutsPOI->GetCentrality(fAOD), iMulti );
207   if(!fCutsPOI->IsEventSelected(fAOD)) return;
208   if(fCutsPOI->IsEventSelectedInCentrality(fAOD)>0) return;
209   fEvent[1]->Fill( fCutsPOI->GetCentrality(fAOD), iMulti );
210
211   if (fDebugV2) printf("Event selected\n");
212   fCutsRP->SetEvent(fAOD,MCEvent());
213   AliFlowTrackCuts* dummy = new AliFlowTrackCuts("null_cuts");
214   dummy->SetParamType(AliFlowTrackCuts::kGlobal);
215   dummy->SetPtRange(+1,-1); // select nothing QUICK
216   dummy->SetEtaRange(+1,-1); // select nothing VZERO
217   dummy->SetEvent(fAOD,MCEvent());
218   AliFlowEvent *flowEvent[5];
219   for(int r=0; r!=5; ++r) {
220     flowEvent[r] = new AliFlowEvent(fCutsRP,dummy);
221     flowEvent[r]->SetReferenceMultiplicity( iMulti );
222     flowEvent[r]->DefineDeadZone(0,0,0,0);
223     flowEvent[r]->TagSubeventsInEta( fFlowEta[0], fFlowEta[1],
224                                      fFlowEta[2], fFlowEta[3] );
225   }
226   delete dummy;
227   if (fDebugV2) printf(" ᶫFlowEvent has %d RPs\n", flowEvent[0]->NumberOfTracks() );
228
229   switch (fSource) {
230     case (AliRDHFCuts::kD0toKpiCuts):
231       FillD0toKpi(fAOD,flowEvent); break;
232     case (AliRDHFCuts::kDstarCuts):
233       FillDStartoKpipi(fAOD,flowEvent); break;
234     case (AliRDHFCuts::kDplusCuts):
235       FillDplustoKpipi(fAOD,flowEvent); break;
236   }
237
238   for(int m=0; m!=5; ++m)
239     PostData(2+m,flowEvent[m]);
240   PostData(1,fHList);
241
242 }
243 //______________________________________________________________________________
244 void AliAnalysisTaskFlowD2H::FillD0toKpi(AliAODEvent *theAOD, 
245                                          AliFlowEvent *theMB[5] )
246 {
247   // Fill D0->Kpi histos
248   TList *listHF = (TList*) theAOD->GetList();
249   if(!listHF) return;
250   TClonesArray *listDzero = (TClonesArray*) listHF->FindObject("D0toKpi");
251   if(!listDzero) return;
252   int nEntries = listDzero->GetEntriesFast();
253   if( !nEntries ) return;
254   AliRDHFCutsD0toKpi *fCutsD0toKpi = (AliRDHFCutsD0toKpi*) fCutsPOI;
255   if (fDebugV2) printf("  ᶫ%d candidates found. Looping...\n", nEntries);
256   const Int_t ndght=2;
257   Int_t nIDs[ndght];
258   for( int iEntry=0; iEntry!=nEntries; ++iEntry ) {
259     AliAODRecoDecayHF2Prong *d0cand = (AliAODRecoDecayHF2Prong*) listDzero->UncheckedAt( iEntry );
260     if( !d0cand ) continue;
261     if( !d0cand->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts) ) continue;
262     if( !fCutsD0toKpi->IsInFiducialAcceptance(d0cand->Pt(),d0cand->Y(421)) )continue;
263     int topCut  = fCutsD0toKpi->IsSelected( d0cand, AliRDHFCuts::kAll, theAOD );
264     int nLevel=topCut>0?1:0;
265     if( (d0cand->Eta()<fPOIEta[0])||(d0cand->Eta()>fPOIEta[1]) )
266       nLevel=0;
267     Double_t massD0=topCut>1?d0cand->InvMassD0bar():d0cand->InvMassD0();
268     for(int h=0; h!=nLevel+1; ++h)
269       fMass[h]->Fill(massD0,d0cand->Pt());
270     if( (d0cand->Pt()<fFlowPts[0]) || (d0cand->Pt()>fFlowPts[1]) ) continue;
271     AliAODTrack* iT;
272     for(Int_t i=0; i!=ndght; ++i) {
273       iT = (AliAODTrack*)d0cand->GetDaughter(i);
274       nIDs[i] = iT->GetID();
275     }
276     // Candidates Insertion (done in filling method: faster)
277     if(nLevel)
278       for(Int_t r=0; r!=5; ++r)
279         if( (massD0>=fFlowBands[0][r]) && (massD0<fFlowBands[1][r]) ) {
280           AliFlowCandidateTrack *sTrack = (AliFlowCandidateTrack*)
281                                  MakeTrack(massD0, d0cand->Pt(), d0cand->Phi(),
282                                            d0cand->Eta(), ndght, nIDs);
283           if(fDebugV2) printf("   ᶫInjecting D0 candidate on band %d \n", r);
284           for(Int_t iDau=0; iDau!=ndght; ++iDau)
285             for(Int_t iRPs=0; iRPs!=theMB[r]->NumberOfTracks(); ++iRPs) {
286               AliFlowTrack *iRP = (AliFlowTrack*) (theMB[r]->GetTrack(iRPs));
287               if(!iRP->InRPSelection()) continue;
288               if( TMath::Abs(sTrack->GetIDDaughter(iDau)) == TMath::Abs(iRP->GetID()) ) {
289                 sTrack->SetDaughter(iDau,iRP);
290                 iRP->SetForRPSelection(kFALSE);
291                 if(fDebugV2) printf("    ᶫdaughter%d with fID %d was removed from this RP set\n", iDau, sTrack->GetIDDaughter(iDau));
292               }
293             }
294           theMB[r]->AddTrack(sTrack);
295         }
296     // END of Candidates Insertion
297   }
298 }
299 //______________________________________________________________________________
300 void AliAnalysisTaskFlowD2H::FillDStartoKpipi(const AliAODEvent *theAOD, 
301                                               AliFlowEvent *theMB[5] )
302 {
303   // Fills D* histos
304   TList *listHF = (TList*) theAOD->GetList();
305   if(!listHF) return;
306   TClonesArray *listDstar = (TClonesArray*) listHF->FindObject("Dstar");
307   if(!listDstar) return;
308   int nEntries = listDstar->GetEntriesFast();
309   if( !nEntries ) return;
310   AliRDHFCutsDStartoKpipi *fCutsDStartoKpipi = (AliRDHFCutsDStartoKpipi*) fCutsPOI;
311   if (fDebugV2) printf("  ᶫ%d candidates found. Looping...\n", nEntries);
312   const Int_t ndght=3;
313   Int_t nIDs[ndght];
314   for( int iEntry=0; iEntry!=nEntries; ++iEntry ) {
315     AliAODRecoCascadeHF *dst = (AliAODRecoCascadeHF*) listDstar->UncheckedAt( iEntry );
316     if( !dst ) continue;
317     if( !dst->HasSelectionBit(AliRDHFCuts::kDstarCuts) ) continue;
318     if( !fCutsDStartoKpipi->IsInFiducialAcceptance(dst->Pt(),dst->YDstar()) )continue;
319     int topCut = fCutsDStartoKpipi->IsSelected( dst, AliRDHFCuts::kAll );
320     int nLevel=topCut>0?1:0;
321     if( (dst->Eta()<fPOIEta[0])||(dst->Eta()>fPOIEta[1]) )
322       nLevel=0;
323     Double_t massDS=dst->DeltaInvMass();
324     for(int h=0; h!=nLevel+1; ++h)
325       fMass[h]->Fill(massDS,dst->Pt());
326     if( (dst->Pt()<fFlowPts[0]) || (dst->Pt()>fFlowPts[1]) ) continue;
327     AliAODRecoDecayHF2Prong *d0cand = (AliAODRecoDecayHF2Prong*)dst->Get2Prong();
328     nIDs[0] = ((AliAODTrack*)d0cand->GetDaughter(0))->GetID();
329     nIDs[1] = ((AliAODTrack*)d0cand->GetDaughter(1))->GetID();
330     nIDs[2] = ((AliAODTrack*)dst->GetBachelor() )->GetID();
331     // Candidates Insertion (done in filling method: faster)
332     if(nLevel)
333       for(Int_t r=0; r!=5; ++r)
334         if( (massDS>=fFlowBands[0][r]) && (massDS<fFlowBands[1][r]) ) {
335           AliFlowCandidateTrack *sTrack = (AliFlowCandidateTrack*)
336                                  MakeTrack(massDS, dst->Pt(), dst->Phi(),
337                                            dst->Eta(), ndght, nIDs);
338           if(fDebugV2) printf("   ᶫInjecting DStar candidate on band %d \n", r);
339           for(Int_t iDau=0; iDau!=ndght; ++iDau)
340             for(Int_t iRPs=0; iRPs!=theMB[r]->NumberOfTracks(); ++iRPs) {
341               AliFlowTrack *iRP = (AliFlowTrack*) (theMB[r]->GetTrack(iRPs));
342               if(!iRP->InRPSelection()) continue;
343               if( TMath::Abs(sTrack->GetIDDaughter(iDau)) == TMath::Abs(iRP->GetID()) ) {
344                 sTrack->SetDaughter(iDau,iRP);
345                 iRP->SetForRPSelection(kFALSE);
346                 if(fDebugV2) printf("    ᶫdaughter%d with fID %d was removed from this RP set\n", iDau, sTrack->GetIDDaughter(iDau));
347               }
348             }
349           theMB[r]->AddTrack(sTrack);
350         }
351     // END of Candidates Insertion
352   }
353 }
354 //______________________________________________________________________________
355 void AliAnalysisTaskFlowD2H::FillDplustoKpipi(const AliAODEvent *theAOD, 
356                                               AliFlowEvent *theMB[5] )
357 {
358   // Fill D+ histos
359   TList *listHF = (TList*) theAOD->GetList();
360   if(!listHF) return;
361   TClonesArray *listDplus = (TClonesArray*) listHF->FindObject("Charm3Prong");
362   if(!listDplus) return;
363   int nEntries = listDplus->GetEntriesFast();
364   if( !nEntries ) return;
365   AliRDHFCutsDplustoKpipi *fCutsDStartoKpipi = (AliRDHFCutsDplustoKpipi*) fCutsPOI;
366   if (fDebugV2) printf("  ᶫ%d candidates found. Looping...\n", nEntries);
367   const Int_t ndght=3;
368   Int_t nIDs[ndght];
369   for( int iEntry=0; iEntry!=nEntries; ++iEntry ) {
370     AliAODRecoDecayHF3Prong *dplu = (AliAODRecoDecayHF3Prong*) listDplus->UncheckedAt( iEntry );
371     if( !dplu ) continue;
372     if( !dplu->HasSelectionBit(AliRDHFCuts::kDplusCuts) ) continue;
373     if( !fCutsDStartoKpipi->IsInFiducialAcceptance(dplu->Pt(),dplu->YDplus()) )continue;
374     int topCut = fCutsDStartoKpipi->IsSelected( dplu, AliRDHFCuts::kAll );
375     int nLevel=topCut>0?1:0;
376     if( (dplu->Eta()<fPOIEta[0])||(dplu->Eta()>fPOIEta[1]) )
377       nLevel=0;
378     Double_t massDp=dplu->InvMassDplus();
379     for(int h=0; h!=nLevel+1; ++h)
380       fMass[h]->Fill(massDp,dplu->Pt());
381     if( (dplu->Pt()<fFlowPts[0]) || (dplu->Pt()>fFlowPts[1]) ) continue;
382     nIDs[0] = ((AliAODTrack*)dplu->GetDaughter(0))->GetID();
383     nIDs[1] = ((AliAODTrack*)dplu->GetDaughter(1))->GetID();
384     nIDs[2] = ((AliAODTrack*)dplu->GetDaughter(2))->GetID();
385     // Candidates Insertion (done in filling method: faster)
386     if(nLevel)
387       for(Int_t r=0; r!=5; ++r)
388         if( (massDp>=fFlowBands[0][r]) && (massDp<fFlowBands[1][r]) ) {
389           AliFlowCandidateTrack *sTrack = (AliFlowCandidateTrack*)
390                                  MakeTrack(massDp, dplu->Pt(), dplu->Phi(),
391                                            dplu->Eta(), ndght, nIDs);
392           if(fDebugV2) printf("   ᶫInjecting DStar candidate on band %d \n", r);
393           for(Int_t iDau=0; iDau!=ndght; ++iDau)
394             for(Int_t iRPs=0; iRPs!=theMB[r]->NumberOfTracks(); ++iRPs) {
395               AliFlowTrack *iRP = (AliFlowTrack*) (theMB[r]->GetTrack(iRPs));
396               if(!iRP->InRPSelection()) continue;
397               if( TMath::Abs(sTrack->GetIDDaughter(iDau)) == TMath::Abs(iRP->GetID()) ) {
398                 sTrack->SetDaughter(iDau,iRP);
399                 iRP->SetForRPSelection(kFALSE);
400                 if(fDebugV2) printf("    ᶫdaughter%d with fID %d was removed from this RP set\n", iDau, sTrack->GetIDDaughter(iDau));
401               }
402             }
403           theMB[r]->AddTrack(sTrack);
404         }
405     // END of Candidates Insertion
406   }
407 }
408 //______________________________________________________________________________
409 AliFlowCandidateTrack* AliAnalysisTaskFlowD2H::MakeTrack( Double_t mass, 
410                           Double_t pt, Double_t phi, Double_t eta, 
411                           Int_t nDau, const Int_t *iID ) {
412   // create track for flow tasks
413   AliFlowCandidateTrack *sTrack = new AliFlowCandidateTrack();
414   sTrack->SetMass(mass);
415   sTrack->SetPt(pt);
416   sTrack->SetPhi(phi);
417   sTrack->SetEta(eta);
418   for(Int_t iDau=0; iDau!=nDau; ++iDau)
419     sTrack->AddDaughter(iID[iDau]);
420   sTrack->SetForPOISelection(kTRUE);
421   sTrack->SetForRPSelection(kFALSE);
422   return sTrack;
423 }
424 //_____________________________________________________________________________
425 void AliAnalysisTaskFlowD2H::Terminate(Option_t *)
426 {
427   // Close the analysis, empty for now
428 }
429