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     // Fix references to daughter tracks
 134     //AliAnalysisVertexingHF *fixer = new AliAnalysisVertexingHF();
 135     //fixer->FixReferences(aod);
 136     //delete fixer;
 137     //
 138     //AliRDHFCutsD0toKpi *mycuts=new AliRDHFCutsD0toKpi();
 139     //mycuts->SetFixRefs(kTRUE);
 140     //mycuts->IsEventSelected(aod);
 141
 142     // loop over D0->Kpi candidates
 143     Int_t nD0toKpi = arrayD0toKpi->GetEntriesFast();
 144     nTotD0toKpi += nD0toKpi;
 145     cout<<"Number of D0->Kpi: "<<nD0toKpi<<endl;
 146
 147     for (Int_t iD0toKpi = 0; iD0toKpi < nD0toKpi; iD0toKpi++) {
 148       AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->UncheckedAt(iD0toKpi);
 149
 150         // get daughter AOD tracks
 151         AliAODTrack *trk0 = (AliAODTrack*)d->GetDaughter(0);
 152         AliAODTrack *trk1 = (AliAODTrack*)d->GetDaughter(1);
 153
 154         if(trk0->GetStatus()) printf("ok %d\n",iD0toKpi);
 155     }
 156
 157
 158     // loop over D* candidates
 159     Int_t nDstar = arrayDstar->GetEntriesFast();
 160     nTotDstar += nDstar;
 161     cout<<"Number of D*->D0pi: "<<nDstar<<endl;
 162
 163   AliRDHFCutsDStartoKpipi *cuts = new AliRDHFCutsDStartoKpipi("CutsDStartoKpipi");
 164
 165
 166  AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
 167  //esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 168   //default
 169  esdTrackCuts->SetRequireTPCRefit(kTRUE);
 170   //esdTrackCuts->SetRequireITSRefit(kTRUE);
 171   //esdTrackCuts->SetMinNClustersTPC(70);
 172   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 173               AliESDtrackCuts::kAny);
 174  // default is kBoth, otherwise kAny
 175   esdTrackCuts->SetMinDCAToVertexXY(0.);
 176   esdTrackCuts->SetPtRange(0.3,1.e10);
 177
 178  // soft pion pre-selections
 179   AliESDtrackCuts* esdSoftPicuts=new AliESDtrackCuts();
 180   //esdSoftPicuts->SetRequireSigmaToVertex(kFALSE);
 181   //default
 182   //esdSoftPicuts->SetRequireTPCRefit(kFALSE);
 183   //esdSoftPicuts->SetRequireITSRefit(kFALSE);
 184   //esdSoftPicuts->SetMinNClustersITS(4); // default is 5
 185   //esdSoftPicuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 186   //                  AliESDtrackCuts::kAny); //test d0 asimmetry
 187   //esdSoftPicuts->SetPtRange(0.0,15.0);
 188
 189   cuts->AddTrackCuts(esdTrackCuts);
 190   cuts->AddTrackCutsSoftPi(esdSoftPicuts);
 191   Float_t cutsArrayDStartoKpipi[14]={0.15,0.07.,0.85,0.8,0.8,0.06.,0.06.,0.001,0.6,0.15, 0.03, 0.2, 5, 0.5}; // first 9 for D0 from D*, last 5 for D*
 192   cuts->SetCuts(14,cutsArrayDStartoKpipi);
 193
 194     for (Int_t iDstar = 0; iDstar < nDstar; iDstar++) {
 195       AliAODRecoCascadeHF *c = (AliAODRecoCascadeHF*)arrayDstar->UncheckedAt(iDstar);
 196
 197       if(cuts->IsSelected(c,AliRDHFCuts::kTracks)) printf("passed\n");
 198
 199     }
 200
 201
 202     // count 3prong candidates
 203     Int_t n3Prong = array3Prong->GetEntriesFast();
 204     nTot3Prong += n3Prong;
 205     cout<<"Number of Charm->3Prong: "<<n3Prong<<endl;
 206
 207     // loop over HF vertices
 208     Int_t nVtxsHF = arrayVerticesHF->GetEntriesFast();
 209     nTotHF += nVtxsHF;
 210     cout<<"Number of heavy-flavour vertices: "<<nVtxsHF<<endl;
 211     for (Int_t iVtx = 0; iVtx <0; iVtx++) {
 212       AliAODVertex *vtxHF = (AliAODVertex*)arrayVerticesHF->UncheckedAt(iVtx);
 213       // print info
 214       //cout << iVtx << ": vertex z position: " << vtxHF->GetZ() << endl;
 215     }
 216
 217
 218     // count cascade candidates
 219     if (arrayCascades){
 220       Int_t nCasc = arrayCascades->GetEntriesFast();
 221       nTotCasc+=nCasc;
 222       cout << "Number of Cascades: "<<nCasc<<endl;
 223     }
 224
 225   }
 226
 227   printf("\n Total HF vertices: %d\n",nTotHF);
 228   printf("\n Total D0->Kpi: %d\n",nTotD0toKpi);
 229   printf("\n Total D*->D0pi: %d\n",nTotDstar);
 230   printf("\n Total Charm->3Prong: %d\n",nTot3Prong);
 231   if (arrayCascades) printf("\n Total Cascades: %d\n",nTotCasc);
 232
 233   TCanvas *c1 = new TCanvas("c1","c1",0,0,800,700);
 234   c1->Divide(2,2);
 235   c1->cd(1);
 236   hCPtaVSd0d0->Draw("colz");
 237   c1->cd(2);
 238   hMass->SetFillColor(4);
 239   hMass->Draw();
 240   c1->cd(3);
 241   hSecVtxZ->SetFillColor(2);
 242   hSecVtxZ->Draw();
 243   c1->cd(4);
 244   hDeltaMassDstar->SetFillColor(3);
 245   hDeltaMassDstar->Draw();
 246
 247   return;
 248 }
 249 //------------------------------------------------------------------------
 250 void ReadAODVertexingHFsa(const char *aodHFFileName="AliAOD.VertexingHF.sa.root")
 251 {
 252   //
 253   // Example macro to read D0->Kpi candidates from a stand-alone
 254   // heavy-flavour AOD (i.e. without standard AOD) and apply cuts
 255   // Origin: A.Dainese
 256   //
 257   Bool_t useParFiles=kFALSE;
 258   gROOT->LoadMacro("$ALICE_ROOT/PWG3/vertexingHF/LoadLibraries.C");
 259   LoadLibraries(useParFiles);
 260
 261   // create a test histogram
 262   TH2F *hCPtaVSd0d0 = new TH2F("hCPtaVSd0d0","D^{0} correlation plot",1000,-50000,50000,1000,-1,1);
 263   hCPtaVSd0d0->SetXTitle("Product of impact parameters [#mu m^{2}]");
 264   hCPtaVSd0d0->SetYTitle("Cosine of pointing angle");
 265   TH1F *hMass = new TH1F("hMass","D^{0} mass plot",100,1.7,2);
 266   hMass->SetXTitle("Invariant mass [GeV]");
 267   hMass->SetYTitle("Entries");
 268   TH1F *hSecVtxZ = new TH1F("hSecVtxZ","D^{0} decay vertex z",1000,-10,10);
 269   hSecVtxZ->SetXTitle("z of decay vertex [cm]");
 270   hSecVtxZ->SetYTitle("Entries");
 271
 272   // open input file and get the TTree
 273   TFile inFile(aodHFFileName,"READ");
 274   if (!inFile.IsOpen()) return;
 275
 276   TTree *aodTree = (TTree*)inFile.Get("aodTree");
 277
 278   AliAODEvent *aod = new AliAODEvent();
 279
 280   aod->ReadFromTree(aodTree);
 281
 282   // load heavy flavour vertices
 283   TClonesArray *arrayVerticesHF =
 284     (TClonesArray*)aod->GetList()->FindObject("VerticesHF");
 285
 286   // load D0->Kpi candidates
 287   TClonesArray *arrayD0toKpi =
 288     (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
 289
 290
 291   Double_t cutsD0[9]=
 292   // cutsD0[0] = inv. mass half width [GeV]
 293   // cutsD0[1] = dca [cm]
 294   // cutsD0[2] = cosThetaStar
 295   // cutsD0[3] = pTK [GeV/c]
 296   // cutsD0[4] = pTPi [GeV/c]
 297   // cutsD0[5] = d0K [cm]   upper limit!
 298   // cutsD0[6] = d0Pi [cm]  upper limit!
 299   // cutsD0[7] = d0d0 [cm^2]
 300   // cutsD0[8] = cosThetaPoint
 301                      {1000.,
 302                       100000.,
 303                       1.1,
 304                       0.,
 305                       0.,
 306                       100000.,
 307                       100000.,
 308                       100000000.,
 309                       -1.1};
 310
 311   Int_t nTotHF=0,nTotD0toKpi=0;
 312   AliAODVertex *vtx1=0;
 313
 314   // loop over events
 315   Int_t nEvents = aodTree->GetEntries();
 316   for (Int_t nEv = 0; nEv < nEvents; nEv++) {
 317     cout<<"\n------------ Event: "<<nEv<<" ------------------"<<endl;
 318
 319     // read event
 320     aodTree->GetEvent(nEv);
 321     //aodTree->BranchRef();
 322
 323
 324     // loop over D0->Kpi candidates
 325     Int_t nD0toKpi = arrayD0toKpi->GetEntriesFast();
 326     nTotD0toKpi += nD0toKpi;
 327     cout<<"Number of D0->Kpi: "<<nD0toKpi<<endl;
 328
 329     for (Int_t iD0toKpi = 0; iD0toKpi < nD0toKpi; iD0toKpi++) {
 330       AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->UncheckedAt(iD0toKpi);
 331       Int_t okD0=0,okD0bar=0;
 332       if(d->SelectD0(cutsD0,okD0,okD0bar)) {
 333         //cout<<1e8*d->Prodd0d0()<<endl;
 334         hMass->Fill(d->InvMassD0(),0.5);
 335         hMass->Fill(d->InvMassD0bar(),0.5);
 336         hCPtaVSd0d0->Fill(1e8*d->Prodd0d0(),d->CosPointingAngle());
 337         hSecVtxZ->Fill(d->GetSecVtxZ());
 338         //cout<<d->GetSecVtxZ() <<endl;
 339
 340       }
 341     }
 342
 343   }
 344
 345   printf("\n Total D0->Kpi: %d\n",nTotD0toKpi);
 346
 347   TCanvas *c = new TCanvas("c","c",0,0,1000,1000);
 348   c->Divide(2,2);
 349   c->cd(1);
 350   hCPtaVSd0d0->Draw("colz");
 351   c->cd(2);
 352   hMass->SetFillColor(4);
 353   hMass->Draw();
 354   c->cd(3);
 355   hSecVtxZ->SetFillColor(2);
 356   hSecVtxZ->Draw();
 357
 358   return;
 359 }