ade99e1c0012db2014562a31967fd2ee62ce7093
[u/mrichter/AliRoot.git] / EVE / alice-macros / aod_HF.C
1 // $Id$
2 // Main author: Davide Caffarri 2009
3
4 /**************************************************************************
5  * Copyright(c) 1998-2008, ALICE Experiment at CERN, all rights reserved. *
6  * See http://aliceinfo.cern.ch/Offline/AliRoot/License.html for          *
7  * full copyright notice.                                                 *
8  **************************************************************************/
9 #if !defined(__CINT__) || defined(__MAKECINT__)
10 #include <TClonesArray.h>
11 #include <TROOT.h>
12 #include <TSystem.h>
13 #include <TEveVSDStructs.h>
14 #include <TEveTrackPropagator.h>
15 #include <TEvePointSet.h>
16 #include <TEveManager.h>
17 #include <TEveUtil.h>
18
19 #include <AliExternalTrackParam.h>
20 #include <AliVVertex.h>
21 #include <AliAODVertex.h>
22 #include <AliAODEvent.h>
23 #include <AliAODTrack.h>
24 #include <AliAODMCParticle.h>
25 #include <AliESDtrack.h>
26 #include <AliESDEvent.h>
27 #include <PWG3/vertexingHF/AliAODRecoDecayHF.h>
28 #include <PWG3/vertexingHF/AliAODRecoDecayHF2Prong.h>
29 #include <PWG3/vertexingHF/macros/LoadLibraries.C>
30 #include <AliEveHF.h>
31 #include <AliEveEventManager.h>
32 #else
33 class AliAODRecoDecayHF;
34 #endif
35
36 void aod_hf_init_rectrack(TEveRecTrack& rt, AliExternalTrackParam* tp)
37 {
38   Double_t      pbuf[3], vbuf[3];
39   rt.fSign = tp->GetSign();
40   tp->GetXYZ(vbuf);     rt.fV.Set(vbuf);
41   tp->GetPxPyPz(pbuf);  rt.fP.Set(pbuf);
42   // Double_t ep = at->GetP(), mc = at->GetMass();
43   rt.fBeta = 1; // ep/TMath::Sqrt(ep*ep + mc*mc);
44 }
45
46
47 AliEveHF* aod_make_HF(TEveTrackPropagator* rnrStyle, AliAODVertex* primAODVtx,
48                       AliESDtrack* neg, AliESDtrack* pos, AliAODRecoDecayHF* rd, Int_t i)
49 {
50   TEveRecTrack  rcPos;
51   TEveRecTrack  rcNeg;
52   //TEveRecV0     rcV0;
53
54   /*
55     Double_t p[3];
56     p[0]=rd->PxProng(0);  p[1]=rd->PyProng(0);  p[2]=rd->PzProng(0);
57     rcV0.fPPos.Set(p);
58     p[0]=rd->PxProng(1);  p[1]=rd->PyProng(1);  p[2]=rd->PzProng(1);
59     rcV0.fPNeg.Set(p);
60
61     p[0]=rd->Px();  p[1]=rd->Py();  p[2]=rd->Pz();
62   */
63
64   Double_t primVtx[3]={primAODVtx->GetX(), primAODVtx->GetY(), primAODVtx->GetZ()};
65
66
67   Double_t v[3] = {rd->Xv(),rd->Yv(),rd->Zv()};
68   printf("vertex %f %f %f\n",v[0],v[1],v[2]);
69
70   aod_hf_init_rectrack(rcNeg, neg);
71   //rcNeg.fIndex = v0->GetNindex();
72   aod_hf_init_rectrack(rcPos, pos);
73   //rcPos.fIndex = v0->GetPindex();
74
75   AliEveHF* myHF = new AliEveHF(&rcNeg, &rcPos, primVtx ,rd, rnrStyle);
76   myHF->SetElementName(Form("D0->Kpi %d", i));
77   myHF->SetElementTitle(Form("CosPointingAngle %f", rd->CosPointingAngle()));
78   myHF->SetAODIndex(i);
79   return myHF;
80 }
81
82 AliEveHFList* aod_HF()
83 {
84   Bool_t useParFiles=kFALSE;
85   
86   TEveUtil::LoadMacro("$ALICE_ROOT/PWG3/vertexingHF/macros/LoadLibraries.C+");
87   LoadLibraries(useParFiles);
88   
89   AliAODEvent* aod = AliEveEventManager::AssertAOD();
90   AliESDEvent* esd = AliEveEventManager::AssertESD();
91
92   /*
93     gSystem->Load("libANALYSIS");
94     gSystem->Load("libANALYSISalice");
95     gSystem->Load("libCORRFW");
96     gSystem->Load("libPWG3base");
97     gSystem->Load("libPWG3vertexingHF");
98   */
99
100   // load MC particles
101   TClonesArray *mcArray =
102     (TClonesArray*) aod->FindListObject(AliAODMCParticle::StdBranchName());
103   if (!mcArray) {
104     printf("MC particles branch not found!\n");
105     return 0;
106   }
107
108   AliAODVertex* primVtx_aod = (AliAODVertex*) aod->GetPrimaryVertex();
109   // AliESDVertex *primVtx_esd = (AliESDVertex*) esd->GetPrimaryVertex();
110
111   AliEveHFList* cont = new AliEveHFList("AOD HF vertices");
112   cont->SetMainColor(2);
113   TEveTrackPropagator* rnrStyle = cont->GetPropagator();
114   rnrStyle->SetMagField( 0.1*aod->GetMagneticField() );
115
116   gEve->AddElement(cont);
117
118   TEvePointSet* pointsD0toKpi = new TEvePointSet("D0->Kpi vertex locations");
119
120   // load D0->Kpi candidates
121   TClonesArray *arrayD0toKpi = (TClonesArray*) aod->FindListObject("D0toKpi");
122
123   // load 3prong candidates
124   // TClonesArray *array3Prong =
125   // (TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
126
127   Int_t countD0 = 0;
128   for (Int_t iD0toKpi=0; iD0toKpi<arrayD0toKpi->GetEntriesFast(); iD0toKpi++)
129   {
130     AliAODRecoDecayHF2Prong *rd = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->UncheckedAt(iD0toKpi);
131     Bool_t unsetvtx=kFALSE;
132     if (!rd->GetOwnPrimaryVtx()) {
133       rd->SetOwnPrimaryVtx(primVtx_aod);
134       unsetvtx=kTRUE;
135     }
136     // REAL D0 particle. If you want to draw only real D0 un-comment these lines
137     //Int_t labD0 = rd->MatchToMC(421,mcArray);
138     //if(labD0<0) continue;
139
140     AliAODTrack *negAODTr = dynamic_cast<AliAODTrack *>(rd->GetDaughter(0));
141     AliAODTrack *posAODTr = dynamic_cast<AliAODTrack *>(rd->GetDaughter(1));
142
143     AliVVertex  *secv = rd->GetSecondaryVtx();
144
145     AliESDtrack *negTr = new AliESDtrack(negAODTr);
146     AliESDtrack *posTr = new AliESDtrack(posAODTr);
147
148     negTr->PropagateToDCA((AliAODVertex*)secv,aod->GetMagneticField(),100.);
149     posTr->PropagateToDCA((AliAODVertex*)secv,aod->GetMagneticField(),100.);
150
151     AliEveHF* myD0 = aod_make_HF(rnrStyle,primVtx_aod,negTr,posTr,rd,iD0toKpi);
152     if (myD0) {
153       gEve->AddElement(myD0,cont);
154       countD0++;
155     }
156
157     pointsD0toKpi->SetNextPoint(rd->Xv(),rd->Yv(),rd->Zv());
158     pointsD0toKpi->SetPointId(rd);
159
160     if(unsetvtx) {
161       rd->UnsetOwnPrimaryVtx();
162     }
163   }
164
165   //cont->SetTitle("test");
166
167   cont->MakeHFs();
168   gEve->Redraw3D();
169
170   pointsD0toKpi->SetTitle(Form("N=%d", pointsD0toKpi->Size()));
171   pointsD0toKpi->SetMarkerStyle(4);
172   pointsD0toKpi->SetMarkerSize(1.5);
173   pointsD0toKpi->SetMarkerColor(kViolet);
174
175   gEve->AddElement(pointsD0toKpi);
176   gEve->Redraw3D();
177
178   return cont;
179 }