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