PWG3/vertexingHF/ReadAODVertexingHF.C

   1 void ReadAODVertexingHF(const char *aodFileName="AliAOD.root",
   2                         const char *aodHFFileName="AliAOD.VertexingHF.root")
   3 {
   4   //
   5   // Example macro to read D0->Kpi candidates from AOD (having the
   6   // standard AOD + a friend heavy-flavour AOD) and apply cuts
   7   // Origin: A.Dainese
   8   //
   9
  10   Bool_t useParFiles=kFALSE;
  11   gROOT->LoadMacro("$ALICE_ROOT/PWG3/vertexingHF/LoadLibraries.C");
  12   LoadLibraries(useParFiles);
  13
  14   // create a test histogram
  15   TH2F *hCPtaVSd0d0 = new TH2F("hCPtaVSd0d0","D^{0} correlation plot",1000,-50000,50000,1000,-1,1);
  16   hCPtaVSd0d0->SetXTitle("Product of impact parameters [#mu m^{2}]");
  17   hCPtaVSd0d0->SetYTitle("Cosine of pointing angle");
  18   TH1F *hMass = new TH1F("hMass","D^{0} mass plot",100,1.7,2);
  19   hMass->SetXTitle("Invariant mass [GeV]");
  20   hMass->SetYTitle("Entries");
  21   TH1F *hSecVtxZ = new TH1F("hSecVtxZ","D^{0} decay vertex z",1000,-10,10);
  22   hSecVtxZ->SetXTitle("z of decay vertex [cm]");
  23   hSecVtxZ->SetYTitle("Entries");
  24   TH1F *hDeltaMassDstar = new TH1F("hDeltaMassDstar","D* delta mass plot",100,0,0.3);
  25   hDeltaMassDstar->SetXTitle("M(Kpipi)-M(Kpi) [GeV]");
  26   hDeltaMassDstar->SetYTitle("Entries");
  27
  28   // open input file and get the TTree
  29   TFile inFile(aodFileName,"READ");
  30   if (!inFile.IsOpen()) return;
  31
  32   TTree *aodTree = (TTree*)inFile.Get("aodTree");
  33   aodTree->AddFriend("aodTree",aodHFFileName);
  34
  35   AliAODEvent *aod = new AliAODEvent();
  36
  37   aod->ReadFromTree(aodTree);
  38
  39   // load heavy flavour vertices
  40   TClonesArray *arrayVerticesHF =
  41     (TClonesArray*)aod->GetList()->FindObject("VerticesHF");
  42
  43   // load D0->Kpi candidates
  44   TClonesArray *arrayD0toKpi =
  45     (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
  46
  47   // load 3prong candidates
  48   TClonesArray *array3Prong =
  49     (TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
  50
  51   // load D* candidates
  52   TClonesArray *arrayDstar =
  53     (TClonesArray*)aod->GetList()->FindObject("Dstar");
  54
  55
  56   Double_t cutsD0[9]=
  57   // cutsD0[0] = inv. mass half width [GeV]
  58   // cutsD0[1] = dca [cm]
  59   // cutsD0[2] = cosThetaStar
  60   // cutsD0[3] = pTK [GeV/c]
  61   // cutsD0[4] = pTPi [GeV/c]
  62   // cutsD0[5] = d0K [cm]   upper limit!
  63   // cutsD0[6] = d0Pi [cm]  upper limit!
  64   // cutsD0[7] = d0d0 [cm^2]
  65   // cutsD0[8] = cosThetaPoint
  66                      {1000.,
  67                       100000.,
  68                       1.1,
  69                       0.,
  70                       0.,
  71                       100000.,
  72                       100000.,
  73                       100000000.,
  74                       -1.1};
  75   Double_t cutsDstar[5]=
  76     // (to be passed to AliAODRecoCascadeHF::SelectDstar())
  77     // 0 = inv. mass half width of D* [GeV]
  78     // 1 = half width of (M_Kpipi-M_Kpi) [GeV]
  79     // 2 = PtMin of pi_s [GeV/c]
  80     // 3 = PtMax of pi_s [GeV/c]
  81     // 4 = theta, angle between the trace of pi_s and D0 decay plane [rad]
  82                      {999999.,
  83                       999999.,
  84                       0.1,
  85                       1.0,
  86                       0.5};
  87
  88
  89   Int_t nTotHF=0,nTotD0toKpi=0,nTotDstar=0,nTot3Prong=0;
  90   AliAODVertex *vtx1=0;
  91
  92   // loop over events
  93   Int_t nEvents = aodTree->GetEntries();
  94   for (Int_t nEv = 0; nEv < nEvents; nEv++) {
  95     cout<<"\n------------ Event: "<<nEv<<" ------------------"<<endl;
  96
  97     // read event
  98     aodTree->GetEvent(nEv);
  99     //aodTree->BranchRef();
 100
 101     //print event info
 102     aod->GetHeader()->Print();
 103
 104     // primary vertex
 105     vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 106     vtx1->Print();
 107
 108     // make trkIDtoEntry register (temporary)
 109     Int_t trkIDtoEntry[100000];
 110     for(Int_t it=0;it<aod->GetNumberOfTracks();it++) {
 111       AliAODTrack *track = aod->GetTrack(it);
 112       trkIDtoEntry[track->GetID()]=it;
 113     }
 114
 115
 116     // loop over D0->Kpi candidates
 117     Int_t nD0toKpi = arrayD0toKpi->GetEntriesFast();
 118     nTotD0toKpi += nD0toKpi;
 119     cout<<"Number of D0->Kpi: "<<nD0toKpi<<endl;
 120
 121     for (Int_t iD0toKpi = 0; iD0toKpi < nD0toKpi; iD0toKpi++) {
 122       AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->UncheckedAt(iD0toKpi);
 123
 124       Bool_t unsetvtx=kFALSE;
 125       if(!d->GetOwnPrimaryVtx()) {
 126         d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 127         unsetvtx=kTRUE;
 128       }
 129       Int_t okD0=0,okD0bar=0;
 130       if(d->SelectD0(cutsD0,okD0,okD0bar)) {
 131         //cout<<1e8*d->Prodd0d0()<<endl;
 132         hMass->Fill(d->InvMassD0(),0.5);
 133         hMass->Fill(d->InvMassD0bar(),0.5);
 134         hCPtaVSd0d0->Fill(1e8*d->Prodd0d0(),d->CosPointingAngle());
 135         hSecVtxZ->Fill(d->GetSecVtxZ());
 136         //cout<<d->GetSecVtxX() <<endl;
 137
 138         // get daughter AOD tracks
 139         AliAODTrack *trk0 = (AliAODTrack*)d->GetDaughter(0);
 140         AliAODTrack *trk1 = (AliAODTrack*)d->GetDaughter(1);
 141         if(!trk0 || !trk1) {
 142           trk0=aod->GetTrack(trkIDtoEntry[d->GetProngID(0)]);
 143           trk1=aod->GetTrack(trkIDtoEntry[d->GetProngID(1)]);
 144           cout<<"references to standard AOD not available"<<endl;
 145         }
 146         cout<<"pt of positive track: "<<trk0->Pt()<<endl;
 147
 148         // make a AliNeutralTrackParam from the D0
 149         // and calculate impact parameters to primary vertex
 150         AliNeutralTrackParam trackD0(d);
 151         cout<<"pt of D0: "<<d->Pt()<<" (AliAODRecoDecay); "<<trackD0.Pt()<<" (track)"<<endl;
 152         //trackD0.Print();
 153         Double_t dz[2],covdz[3];
 154         trackD0.PropagateToDCA(vtx1,aod->GetMagneticField(),1000.,dz,covdz);
 155         cout<<"D0 impact parameter rphi: "<<dz[0]<<" +- "<<TMath::Sqrt(covdz[0])<<endl;
 156       }
 157       if(unsetvtx) d->UnsetOwnPrimaryVtx();
 158     }
 159
 160
 161     // loop over D* candidates
 162     Int_t nDstar = arrayDstar->GetEntriesFast();
 163     nTotDstar += nDstar;
 164     cout<<"Number of D*->D0pi: "<<nDstar<<endl;
 165
 166     for (Int_t iDstar = 0; iDstar < nDstar; iDstar++) {
 167       AliAODRecoCascadeHF *c = (AliAODRecoCascadeHF*)arrayDstar->UncheckedAt(iDstar);
 168       Bool_t unsetvtx=kFALSE;
 169       if(!c->GetOwnPrimaryVtx()) {
 170         c->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 171         c->Get2Prong()->SetOwnPrimaryVtx(vtx1);
 172         unsetvtx=kTRUE;
 173       }
 174       if(c->SelectDstar(cutsDstar,cutsD0)) {
 175         hDeltaMassDstar->Fill(c->DeltaInvMass());
 176         // get daughters
 177         AliAODTrack *trk = (AliAODTrack*)c->GetBachelor();
 178         AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)c->Get2Prong();
 179         cout<<"pt of soft pi: "<<trk->Pt()<<endl;
 180         cout<<"pt of D0: "<<d->Pt()<<endl;
 181       }
 182
 183       if(unsetvtx) {
 184         c->UnsetOwnPrimaryVtx();
 185         c->Get2Prong()->UnsetOwnPrimaryVtx();
 186       }
 187     }
 188
 189
 190     // count 3prong candidates
 191     Int_t n3Prong = array3Prong->GetEntriesFast();
 192     nTot3Prong += n3Prong;
 193     cout<<"Number of Charm->3Prong: "<<n3Prong<<endl;
 194
 195     // loop over HF vertices
 196     Int_t nVtxsHF = arrayVerticesHF->GetEntriesFast();
 197     nTotHF += nVtxsHF;
 198     cout<<"Number of heavy-flavour vertices: "<<nVtxsHF<<endl;
 199     for (Int_t iVtx = 0; iVtx < nVtxsHF; iVtx++) {
 200       AliAODVertex *vtxHF = (AliAODVertex*)arrayVerticesHF->UncheckedAt(iVtx);
 201       // print info
 202       //cout << iVtx << ": vertex z position: " << vtxHF->GetZ() << endl;
 203     }
 204
 205   }
 206
 207   printf("\n Total HF vertices: %d\n",nTotHF);
 208   printf("\n Total D0->Kpi: %d\n",nTotD0toKpi);
 209   printf("\n Total D*->D0pi: %d\n",nTotDstar);
 210   printf("\n Total Charm->3Prong: %d\n",nTot3Prong);
 211
 212   TCanvas *c1 = new TCanvas("c1","c1",0,0,800,700);
 213   c1->Divide(2,2);
 214   c1->cd(1);
 215   hCPtaVSd0d0->Draw("colz");
 216   c1->cd(2);
 217   hMass->SetFillColor(4);
 218   hMass->Draw();
 219   c1->cd(3);
 220   hSecVtxZ->SetFillColor(2);
 221   hSecVtxZ->Draw();
 222   c1->cd(4);
 223   hDeltaMassDstar->SetFillColor(3);
 224   hDeltaMassDstar->Draw();
 225
 226   return;
 227 }
 228 //------------------------------------------------------------------------
 229 void ReadAODVertexingHFsa(const char *aodHFFileName="AliAOD.VertexingHF.sa.root")
 230 {
 231   //
 232   // Example macro to read D0->Kpi candidates from a stand-alone
 233   // heavy-flavour AOD (i.e. without standard AOD) and apply cuts
 234   // Origin: A.Dainese
 235   //
 236   Bool_t useParFiles=kFALSE;
 237   gROOT->LoadMacro("$ALICE_ROOT/PWG3/vertexingHF/LoadLibraries.C");
 238   LoadLibraries(useParFiles);
 239
 240   // create a test histogram
 241   TH2F *hCPtaVSd0d0 = new TH2F("hCPtaVSd0d0","D^{0} correlation plot",1000,-50000,50000,1000,-1,1);
 242   hCPtaVSd0d0->SetXTitle("Product of impact parameters [#mu m^{2}]");
 243   hCPtaVSd0d0->SetYTitle("Cosine of pointing angle");
 244   TH1F *hMass = new TH1F("hMass","D^{0} mass plot",100,1.7,2);
 245   hMass->SetXTitle("Invariant mass [GeV]");
 246   hMass->SetYTitle("Entries");
 247   TH1F *hSecVtxZ = new TH1F("hSecVtxZ","D^{0} decay vertex z",1000,-10,10);
 248   hSecVtxZ->SetXTitle("z of decay vertex [cm]");
 249   hSecVtxZ->SetYTitle("Entries");
 250
 251   // open input file and get the TTree
 252   TFile inFile(aodHFFileName,"READ");
 253   if (!inFile.IsOpen()) return;
 254
 255   TTree *aodTree = (TTree*)inFile.Get("aodTree");
 256
 257   AliAODEvent *aod = new AliAODEvent();
 258
 259   aod->ReadFromTree(aodTree);
 260
 261   // load heavy flavour vertices
 262   TClonesArray *arrayVerticesHF =
 263     (TClonesArray*)aod->GetList()->FindObject("VerticesHF");
 264
 265   // load D0->Kpi candidates
 266   TClonesArray *arrayD0toKpi =
 267     (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
 268
 269
 270   Double_t cutsD0[9]=
 271   // cutsD0[0] = inv. mass half width [GeV]
 272   // cutsD0[1] = dca [cm]
 273   // cutsD0[2] = cosThetaStar
 274   // cutsD0[3] = pTK [GeV/c]
 275   // cutsD0[4] = pTPi [GeV/c]
 276   // cutsD0[5] = d0K [cm]   upper limit!
 277   // cutsD0[6] = d0Pi [cm]  upper limit!
 278   // cutsD0[7] = d0d0 [cm^2]
 279   // cutsD0[8] = cosThetaPoint
 280                      {1000.,
 281                       100000.,
 282                       1.1,
 283                       0.,
 284                       0.,
 285                       100000.,
 286                       100000.,
 287                       100000000.,
 288                       -1.1};
 289
 290   Int_t nTotHF=0,nTotD0toKpi=0;
 291   AliAODVertex *vtx1=0;
 292
 293   // loop over events
 294   Int_t nEvents = aodTree->GetEntries();
 295   for (Int_t nEv = 0; nEv < nEvents; nEv++) {
 296     cout<<"\n------------ Event: "<<nEv<<" ------------------"<<endl;
 297
 298     // read event
 299     aodTree->GetEvent(nEv);
 300     //aodTree->BranchRef();
 301
 302
 303     // loop over D0->Kpi candidates
 304     Int_t nD0toKpi = arrayD0toKpi->GetEntriesFast();
 305     nTotD0toKpi += nD0toKpi;
 306     cout<<"Number of D0->Kpi: "<<nD0toKpi<<endl;
 307
 308     for (Int_t iD0toKpi = 0; iD0toKpi < nD0toKpi; iD0toKpi++) {
 309       AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->UncheckedAt(iD0toKpi);
 310       Int_t okD0=0,okD0bar=0;
 311       if(d->SelectD0(cutsD0,okD0,okD0bar)) {
 312         //cout<<1e8*d->Prodd0d0()<<endl;
 313         hMass->Fill(d->InvMassD0(),0.5);
 314         hMass->Fill(d->InvMassD0bar(),0.5);
 315         hCPtaVSd0d0->Fill(1e8*d->Prodd0d0(),d->CosPointingAngle());
 316         hSecVtxZ->Fill(d->GetSecVtxZ());
 317         //cout<<d->GetSecVtxZ() <<endl;
 318
 319       }
 320     }
 321
 322   }
 323
 324   printf("\n Total D0->Kpi: %d\n",nTotD0toKpi);
 325
 326   TCanvas *c = new TCanvas("c","c",0,0,1000,1000);
 327   c->Divide(2,2);
 328   c->cd(1);
 329   hCPtaVSd0d0->Draw("colz");
 330   c->cd(2);
 331   hMass->SetFillColor(4);
 332   hMass->Draw();
 333   c->cd(3);
 334   hSecVtxZ->SetFillColor(2);
 335   hSecVtxZ->Draw();
 336
 337   return;
 338 }