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   gSystem->Load("libANALYSIS.so");
  10   gSystem->Load("libANALYSISalice.so");
  11   gSystem->Load("libAOD.so");
  12   gSystem->Load("libPWG3base.so");
  13   gSystem->Load("libPWG3vertexingHF.so");
  14
  15   // create a test histogram
  16   TH2F *hCPtaVSd0d0 = new TH2F("hCPtaVSd0d0","D^{0} correlation plot",1000,-50000,50000,1000,-1,1);
  17   hCPtaVSd0d0->SetXTitle("Product of impact parameters [#mu m^{2}]");
  18   hCPtaVSd0d0->SetYTitle("Cosine of pointing angle");
  19   TH1F *hMass = new TH1F("hMass","D^{0} mass plot",100,1.7,2);
  20   hMass->SetXTitle("Invariant mass [GeV]");
  21   hMass->SetYTitle("Entries");
  22   TH1F *hSecVtxZ = new TH1F("hSecVtxZ","D^{0} decay vertex z",1000,-10,10);
  23   hSecVtxZ->SetXTitle("z of decay vertex [cm]");
  24   hSecVtxZ->SetYTitle("Entries");
  25
  26   // open input file and get the TTree
  27   TFile inFile(aodFileName,"READ");
  28   if (!inFile.IsOpen()) return;
  29
  30   TTree *aodTree = (TTree*)inFile.Get("aodTree");
  31   aodTree->AddFriend("aodTree",aodHFFileName);
  32
  33   AliAODEvent *aod = new AliAODEvent();
  34
  35   aod->ReadFromTree(aodTree);
  36
  37   // load heavy flavour vertices
  38   TClonesArray *arrayVerticesHF =
  39     (TClonesArray*)aod->GetList()->FindObject("VerticesHF");
  40
  41   // load D0->Kpi candidates
  42   TClonesArray *arrayD0toKpi =
  43     (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
  44
  45   // load 3prong candidates
  46   TClonesArray *array3Prong =
  47     (TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
  48
  49
  50   Double_t cutsD0[9]=
  51   // cutsD0[0] = inv. mass half width [GeV]
  52   // cutsD0[1] = dca [cm]
  53   // cutsD0[2] = cosThetaStar
  54   // cutsD0[3] = pTK [GeV/c]
  55   // cutsD0[4] = pTPi [GeV/c]
  56   // cutsD0[5] = d0K [cm]   upper limit!
  57   // cutsD0[6] = d0Pi [cm]  upper limit!
  58   // cutsD0[7] = d0d0 [cm^2]
  59   // cutsD0[8] = cosThetaPoint
  60                      {1000.,
  61                       100000.,
  62                       1.1,
  63                       0.,
  64                       0.,
  65                       100000.,
  66                       100000.,
  67                       100000000.,
  68                       -1.1};
  69
  70   Int_t nTotHF=0,nTotD0toKpi=0,nTot3Prong=0;
  71   AliAODVertex *vtx1=0;
  72
  73   // loop over events
  74   Int_t nEvents = aodTree->GetEntries();
  75   for (Int_t nEv = 0; nEv < nEvents; nEv++) {
  76     cout<<"\n------------ Event: "<<nEv<<" ------------------"<<endl;
  77
  78     // read event
  79     aodTree->GetEvent(nEv);
  80     //aodTree->BranchRef();
  81
  82     //print event info
  83     aod->GetHeader()->Print();
  84
  85     // primary vertex
  86     vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
  87     vtx1->Print();
  88
  89     // make trkIDtoEntry register (temporary)
  90     Int_t trkIDtoEntry[100000];
  91     for(Int_t it=0;it<aod->GetNumberOfTracks();it++) {
  92       AliAODTrack *track = aod->GetTrack(it);
  93       trkIDtoEntry[track->GetID()]=it;
  94     }
  95
  96     // loop over D0->Kpi candidates
  97     Int_t nD0toKpi = arrayD0toKpi->GetEntriesFast();
  98     nTotD0toKpi += nD0toKpi;
  99     cout<<"Number of D0->Kpi: "<<nD0toKpi<<endl;
 100
 101     for (Int_t iD0toKpi = 0; iD0toKpi < nD0toKpi; iD0toKpi++) {
 102       AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->UncheckedAt(iD0toKpi);
 103       Bool_t unsetvtx=kFALSE;
 104       if(!d->GetOwnPrimaryVtx()) {
 105         d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 106         unsetvtx=kTRUE;
 107       }
 108       Int_t okD0=0,okD0bar=0;
 109       if(d->SelectD0(cutsD0,okD0,okD0bar)) {
 110         //cout<<1e8*d->Prodd0d0()<<endl;
 111         hMass->Fill(d->InvMassD0(),0.5);
 112         hMass->Fill(d->InvMassD0bar(),0.5);
 113         hCPtaVSd0d0->Fill(1e8*d->Prodd0d0(),d->CosPointingAngle());
 114         hSecVtxZ->Fill(d->GetSecVtxZ());
 115         //cout<<d->GetSecVtxZ() <<endl;
 116
 117         // get daughter AOD tracks
 118         AliAODTrack *trk0=aod->GetTrack(trkIDtoEntry[d->GetProngID(0)]);
 119         AliAODTrack *trk1=aod->GetTrack(trkIDtoEntry[d->GetProngID(1)]);
 120         cout<<"pt of positive track: "<<trk0->Pt()<<endl;
 121         // this will be replaced by
 122         //AliAODTrack *trk0 = (AliAODTrack*)(d->GetSecondaryVtx()->GetDaughter(0));
 123       }
 124       if(unsetvtx) d->UnsetOwnPrimaryVtx();
 125     }
 126
 127     // count 3prong candidates
 128     Int_t n3Prong = array3Prong->GetEntriesFast();
 129     nTot3Prong += n3Prong;
 130     cout<<"Number of Charm->3Prong: "<<n3Prong<<endl;
 131
 132     // loop over HF vertices
 133     Int_t nVtxsHF = arrayVerticesHF->GetEntriesFast();
 134     nTotHF += nVtxsHF;
 135     cout<<"Number of heavy-flavour vertices: "<<nVtxsHF<<endl;
 136     for (Int_t iVtx = 0; iVtx < nVtxsHF; iVtx++) {
 137       AliAODVertex *vtxHF = (AliAODVertex*)arrayVerticesHF->UncheckedAt(iVtx);
 138       // print info
 139       //cout << iVtx << ": vertex z position: " << vtxHF->GetZ() << endl;
 140     }
 141
 142   }
 143
 144   printf("\n Total HF vertices: %d\n",nTotHF);
 145   printf("\n Total D0->Kpi: %d\n",nTotD0toKpi);
 146   printf("\n Total Charm->3Prong: %d\n",nTot3Prong);
 147
 148   TCanvas *c = new TCanvas("c","c",0,0,1000,1000);
 149   c->Divide(2,2);
 150   c->cd(1);
 151   hCPtaVSd0d0->Draw("colz");
 152   c->cd(2);
 153   hMass->SetFillColor(4);
 154   hMass->Draw();
 155   c->cd(3);
 156   hSecVtxZ->SetFillColor(2);
 157   hSecVtxZ->Draw();
 158
 159   return;
 160 }
 161 //------------------------------------------------------------------------
 162 void ReadAODVertexingHFsa(const char *aodHFFileName="AliAOD.VertexingHF.sa.root")
 163 {
 164   //
 165   // Example macro to read D0->Kpi candidates from a stand-alone
 166   // heavy-flavour AOD (i.e. without standard AOD) and apply cuts
 167   // Origin: A.Dainese
 168   //
 169   gSystem->Load("libANALYSIS.so");
 170   gSystem->Load("libANALYSISalice.so");
 171   gSystem->Load("libAOD.so");
 172   gSystem->Load("libPWG3base.so");
 173
 174   // create a test histogram
 175   TH2F *hCPtaVSd0d0 = new TH2F("hCPtaVSd0d0","D^{0} correlation plot",1000,-50000,50000,1000,-1,1);
 176   hCPtaVSd0d0->SetXTitle("Product of impact parameters [#mu m^{2}]");
 177   hCPtaVSd0d0->SetYTitle("Cosine of pointing angle");
 178   TH1F *hMass = new TH1F("hMass","D^{0} mass plot",100,1.7,2);
 179   hMass->SetXTitle("Invariant mass [GeV]");
 180   hMass->SetYTitle("Entries");
 181   TH1F *hSecVtxZ = new TH1F("hSecVtxZ","D^{0} decay vertex z",1000,-10,10);
 182   hSecVtxZ->SetXTitle("z of decay vertex [cm]");
 183   hSecVtxZ->SetYTitle("Entries");
 184
 185   // open input file and get the TTree
 186   TFile inFile(aodHFFileName,"READ");
 187   if (!inFile.IsOpen()) return;
 188
 189   TTree *aodTree = (TTree*)inFile.Get("aodTree");
 190
 191   AliAODEvent *aod = new AliAODEvent();
 192
 193   aod->ReadFromTree(aodTree);
 194
 195   // load heavy flavour vertices
 196   TClonesArray *arrayVerticesHF =
 197     (TClonesArray*)aod->GetList()->FindObject("VerticesHF");
 198
 199   // load D0->Kpi candidates
 200   TClonesArray *arrayD0toKpi =
 201     (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
 202
 203
 204   Double_t cutsD0[9]=
 205   // cutsD0[0] = inv. mass half width [GeV]
 206   // cutsD0[1] = dca [cm]
 207   // cutsD0[2] = cosThetaStar
 208   // cutsD0[3] = pTK [GeV/c]
 209   // cutsD0[4] = pTPi [GeV/c]
 210   // cutsD0[5] = d0K [cm]   upper limit!
 211   // cutsD0[6] = d0Pi [cm]  upper limit!
 212   // cutsD0[7] = d0d0 [cm^2]
 213   // cutsD0[8] = cosThetaPoint
 214                      {1000.,
 215                       100000.,
 216                       1.1,
 217                       0.,
 218                       0.,
 219                       100000.,
 220                       100000.,
 221                       100000000.,
 222                       -1.1};
 223
 224   Int_t nTotHF=0,nTotD0toKpi=0;
 225   AliAODVertex *vtx1=0;
 226
 227   // loop over events
 228   Int_t nEvents = aodTree->GetEntries();
 229   for (Int_t nEv = 0; nEv < nEvents; nEv++) {
 230     cout<<"\n------------ Event: "<<nEv<<" ------------------"<<endl;
 231
 232     // read event
 233     aodTree->GetEvent(nEv);
 234     //aodTree->BranchRef();
 235
 236
 237     // loop over D0->Kpi candidates
 238     Int_t nD0toKpi = arrayD0toKpi->GetEntriesFast();
 239     nTotD0toKpi += nD0toKpi;
 240     cout<<"Number of D0->Kpi: "<<nD0toKpi<<endl;
 241
 242     for (Int_t iD0toKpi = 0; iD0toKpi < nD0toKpi; iD0toKpi++) {
 243       AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->UncheckedAt(iD0toKpi);
 244       Int_t okD0=0,okD0bar=0;
 245       if(d->SelectD0(cutsD0,okD0,okD0bar)) {
 246         //cout<<1e8*d->Prodd0d0()<<endl;
 247         hMass->Fill(d->InvMassD0(),0.5);
 248         hMass->Fill(d->InvMassD0bar(),0.5);
 249         hCPtaVSd0d0->Fill(1e8*d->Prodd0d0(),d->CosPointingAngle());
 250         hSecVtxZ->Fill(d->GetSecVtxZ());
 251         //cout<<d->GetSecVtxZ() <<endl;
 252
 253       }
 254     }
 255
 256   }
 257
 258   printf("\n Total D0->Kpi: %d\n",nTotD0toKpi);
 259
 260   TCanvas *c = new TCanvas("c","c",0,0,1000,1000);
 261   c->Divide(2,2);
 262   c->cd(1);
 263   hCPtaVSd0d0->Draw("colz");
 264   c->cd(2);
 265   hMass->SetFillColor(4);
 266   hMass->Draw();
 267   c->cd(3);
 268   hSecVtxZ->SetFillColor(2);
 269   hSecVtxZ->Draw();
 270
 271   return;
 272 }