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