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       Bool_t unsetvtx=kFALSE;
 124       if(!d->GetOwnPrimaryVtx()) {
 125         d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 126         unsetvtx=kTRUE;
 127       }
 128       Int_t okD0=0,okD0bar=0;
 129       if(d->SelectD0(cutsD0,okD0,okD0bar)) {
 130         //cout<<1e8*d->Prodd0d0()<<endl;
 131         hMass->Fill(d->InvMassD0(),0.5);
 132         hMass->Fill(d->InvMassD0bar(),0.5);
 133         hCPtaVSd0d0->Fill(1e8*d->Prodd0d0(),d->CosPointingAngle());
 134         hSecVtxZ->Fill(d->GetSecVtxZ());
 135         //cout<<d->GetSecVtxX() <<endl;
 136
 137         // get daughter AOD tracks
 138         AliAODTrack *trk0 = (AliAODTrack*)d->GetDaughter(0);
 139         AliAODTrack *trk1 = (AliAODTrack*)d->GetDaughter(1);
 140         if(!trk0 || !trk1) {
 141           trk0=aod->GetTrack(trkIDtoEntry[d->GetProngID(0)]);
 142           trk1=aod->GetTrack(trkIDtoEntry[d->GetProngID(1)]);
 143           cout<<"references to standard AOD not available"<<endl;
 144         }
 145         cout<<"pt of positive track: "<<trk0->Pt()<<endl;
 146
 147         // make a AliNeutralTrackParam from the D0
 148         // and calculate impact parameters to primary vertex
 149         AliNeutralTrackParam trackD0(d);
 150         cout<<"pt of D0: "<<d->Pt()<<" (AliAODRecoDecay); "<<trackD0.Pt()<<" (track)"<<endl;
 151         //trackD0.Print();
 152         Double_t dz[2],covdz[3];
 153         trackD0.PropagateToDCA(vtx1,aod->GetMagneticField(),1000.,dz,covdz);
 154         cout<<"D0 impact parameter rphi: "<<dz[0]<<" +- "<<TMath::Sqrt(covdz[0])<<endl;
 155       }
 156       if(unsetvtx) d->UnsetOwnPrimaryVtx();
 157     }
 158
 159
 160     // loop over D* candidates
 161     Int_t nDstar = arrayDstar->GetEntriesFast();
 162     nTotDstar += nDstar;
 163     cout<<"Number of D*->D0pi: "<<nDstar<<endl;
 164
 165     for (Int_t iDstar = 0; iDstar < nDstar; iDstar++) {
 166       AliAODRecoCascadeHF *c = (AliAODRecoCascadeHF*)arrayDstar->UncheckedAt(iDstar);
 167       Bool_t unsetvtx=kFALSE;
 168       if(!c->GetOwnPrimaryVtx()) {
 169         c->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 170         c->Get2Prong()->SetOwnPrimaryVtx(vtx1);
 171         unsetvtx=kTRUE;
 172       }
 173       if(c->SelectDstar(cutsDstar,cutsD0)) {
 174         hDeltaMassDstar->Fill(c->DeltaInvMass());
 175         // get daughters
 176         AliAODTrack *trk = (AliAODTrack*)c->GetBachelor();
 177         AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)c->Get2Prong();
 178         cout<<"pt of soft pi: "<<trk->Pt()<<endl;
 179         cout<<"pt of D0: "<<d->Pt()<<endl;
 180       }
 181
 182       if(unsetvtx) {
 183         c->UnsetOwnPrimaryVtx();
 184         c->Get2Prong()->UnsetOwnPrimaryVtx();
 185       }
 186     }
 187
 188
 189     // count 3prong candidates
 190     Int_t n3Prong = array3Prong->GetEntriesFast();
 191     nTot3Prong += n3Prong;
 192     cout<<"Number of Charm->3Prong: "<<n3Prong<<endl;
 193
 194     // loop over HF vertices
 195     Int_t nVtxsHF = arrayVerticesHF->GetEntriesFast();
 196     nTotHF += nVtxsHF;
 197     cout<<"Number of heavy-flavour vertices: "<<nVtxsHF<<endl;
 198     for (Int_t iVtx = 0; iVtx < nVtxsHF; iVtx++) {
 199       AliAODVertex *vtxHF = (AliAODVertex*)arrayVerticesHF->UncheckedAt(iVtx);
 200       // print info
 201       //cout << iVtx << ": vertex z position: " << vtxHF->GetZ() << endl;
 202     }
 203
 204   }
 205
 206   printf("\n Total HF vertices: %d\n",nTotHF);
 207   printf("\n Total D0->Kpi: %d\n",nTotD0toKpi);
 208   printf("\n Total D*->D0pi: %d\n",nTotDstar);
 209   printf("\n Total Charm->3Prong: %d\n",nTot3Prong);
 210
 211   TCanvas *c1 = new TCanvas("c1","c1",0,0,800,700);
 212   c1->Divide(2,2);
 213   c1->cd(1);
 214   hCPtaVSd0d0->Draw("colz");
 215   c1->cd(2);
 216   hMass->SetFillColor(4);
 217   hMass->Draw();
 218   c1->cd(3);
 219   hSecVtxZ->SetFillColor(2);
 220   hSecVtxZ->Draw();
 221   c1->cd(4);
 222   hDeltaMassDstar->SetFillColor(3);
 223   hDeltaMassDstar->Draw();
 224
 225   return;
 226 }
 227 //------------------------------------------------------------------------
 228 void ReadAODVertexingHFsa(const char *aodHFFileName="AliAOD.VertexingHF.sa.root")
 229 {
 230   //
 231   // Example macro to read D0->Kpi candidates from a stand-alone
 232   // heavy-flavour AOD (i.e. without standard AOD) and apply cuts
 233   // Origin: A.Dainese
 234   //
 235   Bool_t useParFiles=kFALSE;
 236   gROOT->LoadMacro("$ALICE_ROOT/PWG3/vertexingHF/LoadLibraries.C");
 237   LoadLibraries(useParFiles);
 238
 239   // create a test histogram
 240   TH2F *hCPtaVSd0d0 = new TH2F("hCPtaVSd0d0","D^{0} correlation plot",1000,-50000,50000,1000,-1,1);
 241   hCPtaVSd0d0->SetXTitle("Product of impact parameters [#mu m^{2}]");
 242   hCPtaVSd0d0->SetYTitle("Cosine of pointing angle");
 243   TH1F *hMass = new TH1F("hMass","D^{0} mass plot",100,1.7,2);
 244   hMass->SetXTitle("Invariant mass [GeV]");
 245   hMass->SetYTitle("Entries");
 246   TH1F *hSecVtxZ = new TH1F("hSecVtxZ","D^{0} decay vertex z",1000,-10,10);
 247   hSecVtxZ->SetXTitle("z of decay vertex [cm]");
 248   hSecVtxZ->SetYTitle("Entries");
 249
 250   // open input file and get the TTree
 251   TFile inFile(aodHFFileName,"READ");
 252   if (!inFile.IsOpen()) return;
 253
 254   TTree *aodTree = (TTree*)inFile.Get("aodTree");
 255
 256   AliAODEvent *aod = new AliAODEvent();
 257
 258   aod->ReadFromTree(aodTree);
 259
 260   // load heavy flavour vertices
 261   TClonesArray *arrayVerticesHF =
 262     (TClonesArray*)aod->GetList()->FindObject("VerticesHF");
 263
 264   // load D0->Kpi candidates
 265   TClonesArray *arrayD0toKpi =
 266     (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
 267
 268
 269   Double_t cutsD0[9]=
 270   // cutsD0[0] = inv. mass half width [GeV]
 271   // cutsD0[1] = dca [cm]
 272   // cutsD0[2] = cosThetaStar
 273   // cutsD0[3] = pTK [GeV/c]
 274   // cutsD0[4] = pTPi [GeV/c]
 275   // cutsD0[5] = d0K [cm]   upper limit!
 276   // cutsD0[6] = d0Pi [cm]  upper limit!
 277   // cutsD0[7] = d0d0 [cm^2]
 278   // cutsD0[8] = cosThetaPoint
 279                      {1000.,
 280                       100000.,
 281                       1.1,
 282                       0.,
 283                       0.,
 284                       100000.,
 285                       100000.,
 286                       100000000.,
 287                       -1.1};
 288
 289   Int_t nTotHF=0,nTotD0toKpi=0;
 290   AliAODVertex *vtx1=0;
 291
 292   // loop over events
 293   Int_t nEvents = aodTree->GetEntries();
 294   for (Int_t nEv = 0; nEv < nEvents; nEv++) {
 295     cout<<"\n------------ Event: "<<nEv<<" ------------------"<<endl;
 296
 297     // read event
 298     aodTree->GetEvent(nEv);
 299     //aodTree->BranchRef();
 300
 301
 302     // loop over D0->Kpi candidates
 303     Int_t nD0toKpi = arrayD0toKpi->GetEntriesFast();
 304     nTotD0toKpi += nD0toKpi;
 305     cout<<"Number of D0->Kpi: "<<nD0toKpi<<endl;
 306
 307     for (Int_t iD0toKpi = 0; iD0toKpi < nD0toKpi; iD0toKpi++) {
 308       AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->UncheckedAt(iD0toKpi);
 309       Int_t okD0=0,okD0bar=0;
 310       if(d->SelectD0(cutsD0,okD0,okD0bar)) {
 311         //cout<<1e8*d->Prodd0d0()<<endl;
 312         hMass->Fill(d->InvMassD0(),0.5);
 313         hMass->Fill(d->InvMassD0bar(),0.5);
 314         hCPtaVSd0d0->Fill(1e8*d->Prodd0d0(),d->CosPointingAngle());
 315         hSecVtxZ->Fill(d->GetSecVtxZ());
 316         //cout<<d->GetSecVtxZ() <<endl;
 317
 318       }
 319     }
 320
 321   }
 322
 323   printf("\n Total D0->Kpi: %d\n",nTotD0toKpi);
 324
 325   TCanvas *c = new TCanvas("c","c",0,0,1000,1000);
 326   c->Divide(2,2);
 327   c->cd(1);
 328   hCPtaVSd0d0->Draw("colz");
 329   c->cd(2);
 330   hMass->SetFillColor(4);
 331   hMass->Draw();
 332   c->cd(3);
 333   hSecVtxZ->SetFillColor(2);
 334   hSecVtxZ->Draw();
 335
 336   return;
 337 }