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