0b711409b0993d4a6cb9930cb7e7d9df2faceacc
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowTasks / AliAnalysisTaskFlowK0Candidates.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 // AliAnalysisTaskFlowK0Candidates:
18 // Analysis task to select K0 candidates for flow analysis.
19 // Uses one AliESDtrackCuts object for both daughters and
20 // QA histograms to monitor the reconstruction.
21 // Author: Carlos Perez (cperez@cern.ch)
22 //////////////////////////////////////////////////////
23
24 #include "TChain.h"
25 #include "TList.h"
26 #include "TH1D.h"
27 #include "TVector3.h"
28
29 #include "AliAnalysisTaskSE.h"
30 #include "AliAnalysisManager.h"
31
32 #include "AliESDEvent.h"
33 #include "AliESDInputHandler.h"
34 #include "AliESDtrack.h"
35 #include "AliESDtrackCuts.h"
36
37 #include "AliAODEvent.h"
38 #include "AliAODInputHandler.h"
39 #include "AliAODTrack.h"
40
41 #include "AliCFManager.h"
42
43 #include "TMath.h"
44 #include "TObjArray.h"
45 #include "AliFlowCandidateTrack.h"
46 #include "AliFlowEventCuts.h"
47
48 #include "AliAnalysisTaskFlowK0Candidates.h"
49
50
51 ClassImp(AliAnalysisTaskFlowK0Candidates)
52
53 //_____________________________________________________________________________
54 AliAnalysisTaskFlowK0Candidates::AliAnalysisTaskFlowK0Candidates() :
55   AliAnalysisTaskSE(),
56   fCutsEvent(NULL),
57   fCuts(NULL),
58   fQAList(NULL),
59   fMassMin(0),
60   fMassMax(0),
61   fDLcut(0),
62   tEvent(NULL),
63   tMulti(NULL)
64 {
65   for(int i=0; i!=4; ++i) {
66     tMass[i] = NULL;
67     tDCA[i]  = NULL;
68     tDL[i]   = NULL;
69     tCTP[i]  = NULL;
70     td0d0[i] = NULL;
71     tPhi[i]  = NULL;
72     tEta[i]  = NULL;
73     tPt[i]   = NULL;
74     tAPhi[i] = NULL;
75     tAEta[i] = NULL;
76     tAPt[i]  = NULL;
77     tBPhi[i] = NULL;
78     tBEta[i] = NULL;
79     tBPt[i]  = NULL;
80   }
81 }
82
83 //_____________________________________________________________________________
84 AliAnalysisTaskFlowK0Candidates::AliAnalysisTaskFlowK0Candidates(const char *name, AliFlowEventCuts *cutsEvent, AliESDtrackCuts *cuts, Double_t MassMin, Double_t MassMax) :
85   AliAnalysisTaskSE(name),
86   fCutsEvent(cutsEvent),
87   fCuts(cuts),
88   fQAList(NULL),
89   fMassMin(MassMin),
90   fMassMax(MassMax),
91   fDLcut(0),
92   tEvent(NULL),
93   tMulti(NULL)
94 {
95   for(int i=0; i!=4; ++i) {
96     tMass[i] = NULL;
97     tDCA[i]  = NULL;
98     tDL[i]   = NULL;
99     tCTP[i]  = NULL;
100     td0d0[i] = NULL;
101     tPhi[i]  = NULL;
102     tEta[i]  = NULL;
103     tPt[i]   = NULL;
104     tAPhi[i] = NULL;
105     tAEta[i] = NULL;
106     tAPt[i]  = NULL;
107     tBPhi[i] = NULL;
108     tBEta[i] = NULL;
109     tBPt[i]  = NULL;
110   }
111   DefineInput(  0, TChain::Class() );
112   DefineOutput( 1, TObjArray::Class() );
113   DefineOutput( 2, TList::Class() );
114 }
115
116 //_____________________________________________________________________________
117 AliAnalysisTaskFlowK0Candidates::~AliAnalysisTaskFlowK0Candidates()
118 {
119   if(fQAList) delete fQAList;
120 }
121
122 //_____________________________________________________________________________
123 void AliAnalysisTaskFlowK0Candidates::UserCreateOutputObjects()
124 {
125   fQAList = new TList();
126   fQAList->SetOwner();
127   AddQAEvents();
128   AddQACandidates();
129   PostData( 2, fQAList);
130 }
131 //_____________________________________________________________________________
132 void AliAnalysisTaskFlowK0Candidates::AddQAEvents()
133 {
134   TList *tQAEvents = new TList();
135   tQAEvents->SetName("Events");
136   tQAEvents->SetOwner();
137   tEvent = new TH1D("Event", "Number of Events",    2, 0, 2);
138   tQAEvents->Add(tEvent);
139   tMulti = new TH1D("Multiplicity", "Multiplicity", 180, 0, 10000);
140   tQAEvents->Add(tMulti);
141   fQAList->Add(tQAEvents);
142 }  
143 //_____________________________________________________________________________
144 void AliAnalysisTaskFlowK0Candidates::AddQACandidates()
145 {
146   TList *tQACandidates[4];
147   TList *tQADaughters[4];
148   for(int i=0; i!=4; ++i) {
149     tQACandidates[i] = new TList();
150     tQACandidates[i]->SetOwner();
151     tQACandidates[i]->SetName(Form("Candidates%d",i));
152     tMass[i] = new TH1D( Form("Mass%i",i), "Mass;M_{#pi#pi} [GeV];Counts per MeV", 180, 0.41, 0.59); tQACandidates[i]->Add( tMass[i] );
153     tDCA[i]  = new TH1D( Form("DCA%i" ,i), "DCA;[cm];Counts per 10 um",            180, 0.00, 0.18); tQACandidates[i]->Add( tDCA[i]  );
154     tDL[i]   = new TH1D( Form("DL%i"  ,i), "Decay Length;[cm];Counts per 0.1 mm",  180, 0.00, 1.80); tQACandidates[i]->Add( tDL[i]   );
155     tCTP[i]  = new TH1D( Form("CTP%i" ,i), "Cos#theta_{p}",                        180,-1.10, 1.10); tQACandidates[i]->Add( tCTP[i]  );
156     td0d0[i] = new TH1D( Form("d0d0%i",i), "d_{0}xd_{0};[cm^{2}];Cnts 0.01 mm^{2}",180,-0.009,0.009);tQACandidates[i]->Add( td0d0[i] );
157     tPhi[i]  = new TH1D( Form("Phi%i" ,i), "Phi;[rad];Counts per degree",     180,0,TMath::TwoPi()); tQACandidates[i]->Add( tPhi[i]  );
158     tEta[i]  = new TH1D( Form("Eta%i" ,i), "Eta;;Counts per 0.04",                 180,-3.60, 3.60); tQACandidates[i]->Add( tEta[i]  );
159     tPt[i]   = new TH1D( Form("Pt%i"  ,i), "Pt;[GeV];Counts per 0.1 GeV",          180, 0.00,18.00); tQACandidates[i]->Add( tPt[i]   );
160     tQADaughters[i] = new TList();
161     tQADaughters[i]->SetOwner();
162     tQADaughters[i]->SetName(Form("Daughters%d",i));
163     tAPhi[i] = new TH1D( Form("PhiBef%i",i), "Phi prePropagation;[rad];Cnts per degree",180,0,TMath::TwoPi()); tQADaughters[i]->Add( tAPhi[i] );
164     tAEta[i] = new TH1D( Form("EtaBef%i",i), "Eta prePropagation;;Counts per 0.04"     ,180,-3.6,3.6); tQADaughters[i]->Add( tAEta[i] );
165     tAPt[i]  = new TH1D( Form("PtBef%i" ,i), "Pt prePropagation;[GeV];Counts per 0.1 GeV",180, 0, 18); tQADaughters[i]->Add( tAPt[i]  );
166     tBPhi[i] = new TH1D( Form("PhiAft%i",i), "Phi posPropagation;[rad];Cnts per degree",180,0,TMath::TwoPi()); tQADaughters[i]->Add( tBPhi[i] );
167     tBEta[i] = new TH1D( Form("EtaAft%i",i), "Eta posPropagation;;Counts per 0.04"     ,180,-3.6,3.6); tQADaughters[i]->Add( tBEta[i] );
168     tBPt[i]  = new TH1D( Form("PtAft%i" ,i), "Pt posPropagation;[GeV];Counts per 0.1 GeV",180, 0, 18); tQADaughters[i]->Add( tBPt[i]  );
169     tQACandidates[i]->Add(tQADaughters[i]);
170     fQAList->Add(tQACandidates[i]);
171   }
172 }
173 //_____________________________________________________________________________
174 void AliAnalysisTaskFlowK0Candidates::NotifyRun()
175 {
176 //  AliESDEvent *fESD = dynamic_cast<AliESDEvent*>(InputEvent());
177 //  if(!fESD) return;
178 }
179
180 //_____________________________________________________________________________
181 void AliAnalysisTaskFlowK0Candidates::UserExec(Option_t *)
182 {
183   AliESDEvent *fESD = dynamic_cast<AliESDEvent*>(InputEvent());
184   AliAODEvent *fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
185   if( (!fESD)&&(!fAOD) )
186     return;
187 //  printf("\nEvent found");
188   tEvent->Fill( 0 );
189   if(fCutsEvent)
190     if( !(fCutsEvent->IsSelected(InputEvent())) )
191       return;
192   tEvent->Fill( 1 );
193   if(fESD)
194     ReadFromESD(fESD);
195   else if(fAOD)
196     ReadFromAOD(fAOD);
197 }
198 //_____________________________________________________________________________
199 void AliAnalysisTaskFlowK0Candidates::ReadFromESD(AliESDEvent *fESD)
200 {
201   TObjArray *POIselection =  new TObjArray();
202   POIselection->SetOwner();
203   int nTracks = fESD->GetNumberOfTracks();
204   tMulti->Fill( nTracks );
205   for(int i=0; i!=nTracks; ++i) {
206     AliESDtrack *ioT = fESD->GetTrack(i);
207     if(fCuts)
208       if( !(fCuts->IsSelected(ioT)) )
209         continue;
210     // printf("\n first particle OK...");
211     for(int j=i+1; j!=nTracks; ++j) {
212       AliESDtrack *joT = fESD->GetTrack(j);
213       if( (ioT->Charge()*joT->Charge()) > 0 )
214         continue;
215       if(fCuts)
216         if( !(fCuts->IsSelected(joT)) )
217           continue;
218       // printf("\n second also...");
219       AliESDtrack *iT = new AliESDtrack(*ioT);
220       AliESDtrack *jT = new AliESDtrack(*joT);
221       // getting distance of closest approach
222       double DCA = iT->PropagateToDCA(jT,fESD->GetMagneticField());
223       // printf(" propagated...");
224       // getting decay length
225       TVector3 vp, vv, vl;
226       vp = TVector3( (iT->Xv()+jT->Xv())/2, (iT->Yv()+jT->Yv())/2, (iT->Zv()+jT->Zv())/2 ); // candidate position
227       vv = TVector3( fESD->GetPrimaryVertex()->GetX(), fESD->GetPrimaryVertex()->GetY(), fESD->GetPrimaryVertex()->GetZ() ); // vertex position
228       vl = vp - vv; // flight line
229       double DL;
230       DL  = vl.Mag();
231       // getting cos pointing angle
232       double CTP;
233       TVector3 vi, vj, vs;
234       vi = TVector3( iT->Px(), iT->Py(), iT->Pz() );
235       vj = TVector3( jT->Px(), jT->Py(), jT->Pz() );
236       vs = vi + vj;
237       CTP = TMath::Cos( vl.Angle(vs) ); // Marta says: "Carlos, it is faster if you calculate this by hand!"
238       double D0D0;
239       D0D0 = iT->GetD(vv.X(),vv.Y(),fESD->GetMagneticField())*jT->GetD(vv.X(),vv.Y(),fESD->GetMagneticField());
240       // getting invariant mass
241       double sum12 = iT->P()*iT->P()+jT->P()*jT->P();
242       double pro12 = iT->P()*iT->P()*jT->P()*jT->P();
243       double InvMass = 0;
244       InvMass += TMath::Power(0.13957018,2);
245       InvMass += sqrt( TMath::Power(0.13957018,4) + TMath::Power(0.13957018,2)*(sum12) + pro12 );
246       InvMass -= ( iT->Px()*jT->Px()+iT->Py()*jT->Py()+iT->Pz()*jT->Pz() );
247       InvMass *= 2;
248       InvMass = sqrt(InvMass);
249       // filtering
250       int iLevel = 3;
251       if( CTP<0.8 )
252         iLevel = 2;
253       if( DL<fDLcut ) // 0.5
254         iLevel = 1;
255       if( DCA>0.05 )
256         iLevel = 0;
257       // printf(" candidate at level %d...",iLevel);
258       for(int h=0; h!=iLevel+1; ++h) {
259         // printf(" %d",h);
260         // candidate
261         tDCA[h]->Fill( DCA );
262         tDL[h]->Fill( DL );
263         tCTP[h]->Fill( CTP );
264         td0d0[h]->Fill( D0D0 );
265         tMass[h]->Fill( InvMass );
266         tPhi[h]->Fill( vs.Phi()+TMath::Pi() );
267         tEta[h]->Fill( vs.Eta() );
268         tPt[h]->Fill( vs.Pt() );
269         // daughters
270         tAPhi[h]->Fill( iT->Phi() );
271         tAEta[h]->Fill( iT->Eta() );
272         tAPt[h]->Fill(  iT->Pt()  );
273         tAPhi[h]->Fill( jT->Phi() );
274         tAEta[h]->Fill( jT->Eta() );
275         tAPt[h]->Fill(  jT->Pt()  );
276         tBPhi[h]->Fill( ioT->Phi());
277         tBEta[h]->Fill( ioT->Eta());
278         tBPt[h]->Fill(  ioT->Pt() );
279         tBPhi[h]->Fill( joT->Phi());
280         tBEta[h]->Fill( joT->Eta());
281         tBPt[h]->Fill(  joT->Pt() );
282       }
283       // building candidate
284       if( (InvMass>fMassMin) && (InvMass<fMassMax) && (iLevel==3) ) {
285         AliFlowCandidateTrack *sTrack = new AliFlowCandidateTrack();
286         sTrack->SetMass( InvMass );
287         sTrack->SetPt( vs.Pt() );
288         sTrack->SetPhi( vs.Phi()+TMath::Pi() );
289         sTrack->SetEta( vs.Eta() );
290         sTrack->SetCharge( 0 );
291         sTrack->AddDaughter( iT->GetID() );
292         sTrack->AddDaughter( jT->GetID() );
293         POIselection->AddLast( sTrack );
294       }
295       delete iT;
296       delete jT;
297     } // endfor j
298   } // endfor i
299   PostData( 1, POIselection );
300   PostData( 2, fQAList);
301 }
302
303 //_____________________________________________________________________________
304 void AliAnalysisTaskFlowK0Candidates::ReadFromAOD(AliAODEvent *fAOD)
305 {
306   TObjArray *POIselection =  new TObjArray();
307   POIselection->SetOwner();
308   int nTracks = fAOD->GetNumberOfTracks();
309   tMulti->Fill( nTracks );
310   for(int i=0; i!=fAOD->GetNumberOfV0s(); ++i) {
311     AliAODv0 *iV0 = (AliAODv0*) fAOD->GetV0(i);
312     if ( iV0->GetOnFlyStatus() ) continue;
313     if ( (iV0->ChargeProng(0)*iV0->ChargeProng(1))>0 ) continue;
314     // getting distance of closest approach
315     double DCA = iV0->DcaV0Daughters();
316     // getting decay length
317     double DL;
318     double vv[3];
319     vv[0] = fAOD->GetPrimaryVertex()->GetX();
320     vv[1] = fAOD->GetPrimaryVertex()->GetY();
321     vv[2] = fAOD->GetPrimaryVertex()->GetZ();
322     DL = iV0->DecayLengthV0(vv);
323     // cos pointing angle
324     double CTP;
325     CTP = iV0->CosPointingAngle(vv);
326     double D0D0;
327     D0D0 = iV0->Prodd0d0();
328     // getting invariant mass
329     double InvMass = iV0->MassK0Short();
330     // filtering
331     int iLevel = 3;
332     if( CTP<0.8 )
333       iLevel = 2;
334     if( DL<fDLcut ) // 0.5
335       iLevel = 1;
336     if( DCA>0.05 )
337       iLevel = 0;
338     for(int h=0; h!=iLevel+1; ++h) {
339       // candidate
340       tDCA[h]->Fill( DCA );
341       tDL[h]->Fill( DL );
342       tCTP[h]->Fill( CTP );
343       td0d0[h]->Fill( D0D0 );
344       tMass[h]->Fill( InvMass );
345       tPhi[h]->Fill( iV0->Phi() );
346       tEta[h]->Fill( iV0->Eta() );
347       tPt[h]->Fill( iV0->Pt() );
348       // daughters
349       tAPhi[h]->Fill( ( (AliAODTrack*) iV0->GetDaughter(0) )->Phi() );
350       tAEta[h]->Fill( ( (AliAODTrack*) iV0->GetDaughter(0) )->Eta() );
351       tAPt[h]->Fill(  ( (AliAODTrack*) iV0->GetDaughter(0) )->Pt()  );
352       tAPhi[h]->Fill( ( (AliAODTrack*) iV0->GetDaughter(1) )->Phi() );
353       tAEta[h]->Fill( ( (AliAODTrack*) iV0->GetDaughter(1) )->Eta() );
354       tAPt[h]->Fill(  ( (AliAODTrack*) iV0->GetDaughter(1) )->Pt()  );
355       tBPhi[h]->Fill( iV0->PhiProng(0) );
356       tBEta[h]->Fill( iV0->EtaProng(0) );
357       tBPt[h]->Fill(  iV0->PtProng(0)  );
358       tBPhi[h]->Fill( iV0->PhiProng(1) );
359       tBEta[h]->Fill( iV0->EtaProng(1) );
360       tBPt[h]->Fill(  iV0->PtProng(1)  );
361     }
362     AliFlowCandidateTrack *sTrack = new AliFlowCandidateTrack();
363     sTrack->SetMass( InvMass );
364     sTrack->SetPt( iV0->Pt() );
365     sTrack->SetPhi( iV0->Phi() );
366     sTrack->SetEta( iV0->Eta() );
367     sTrack->SetCharge( 0 );
368     sTrack->AddDaughter( iV0->GetPosID() );
369     sTrack->AddDaughter( iV0->GetNegID() );
370     POIselection->AddLast( sTrack );
371   }
372   PostData( 1, POIselection );
373   PostData( 2, fQAList);
374 }
375
376 //_____________________________________________________________________________
377 void AliAnalysisTaskFlowK0Candidates::Terminate(Option_t *)
378 {
379 }
380