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