PWG3/vertexingHF/CompareBkgLikeSign_JPSI.C

   1 void CompareBkgLikeSign_JPSI(const char *aodFileName="AliAOD.root",
   2                              const char *aodHFFileName="AliAOD.VertexingHF.root")
   3 {
   4   //
   5   // Example macro to read D0->Kpi, JPSI->ee, Like Sign
   6   // candidates from AOD (having the standard AOD +
   7   // a friend heavy-flavour AOD) and apply cuts
   8   // C. Di Giglio
   9   //
  10   Bool_t useParFiles=kFALSE;
  11   gROOT->LoadMacro("$ALICE_ROOT/PWG3/vertexingHF/LoadLibraries.C");
  12   LoadLibraries(useParFiles);
  13
  14   //create histograms
  15   TH2F *hCPtaVSd0d0 = new TH2F("hCPtaVSd0d0","Like Sign pairs 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","Like Sign mass plot",100,0.,3.2);
  19   hMass->SetXTitle("Invariant mass [GeV]");
  20   hMass->SetYTitle("Entries");
  21   hMass->Sumw2();
  22   TH1F *hMassJPSI = new TH1F("hMassJPSI","JPSI mass plot",100,0.,3.2);
  23   hMassJPSI->SetXTitle("Invariant mass [GeV]");
  24   hMassJPSI->SetYTitle("Entries");
  25   hMassJPSI->Sumw2();
  26   TH1F *hSecVtxZ = new TH1F("hSecVtxZ","LikeSign decay vertex z",1000,-10,10);
  27   hSecVtxZ->SetXTitle("z of decay vertex [cm]");
  28   hSecVtxZ->SetYTitle("Entries");
  29   // Cosine theta star
  30   TH1F *hCtsJPSI = new TH1F("hCtsJPSI","JPSI cosine of decay angle",50,-1.,1.);
  31   hCtsJPSI->SetXTitle("cos #theta^{*}");
  32   hCtsJPSI->SetYTitle("Entries");
  33   hCtsJPSI->Sumw2();
  34   TH1F *hCtsLikeSign = new TH1F("hCtsLikeSign","LikeSign cosine of decay angle",50,-1.,1.);
  35   hCtsLikeSign->SetXTitle("cos #theta^{*}");
  36   hCtsLikeSign->SetYTitle("Entries");
  37   hCtsLikeSign->Sumw2();
  38   TH1F *hCtsLikeSignPos = new TH1F("hCtsLikeSignPos","LikeSign cosine of decay angle for ++ pairs",50,-1.,1.);
  39   hCtsLikeSignPos->SetXTitle("cos #theta^{*}");
  40   hCtsLikeSignPos->SetYTitle("Entries");
  41   hCtsLikeSignPos->Sumw2();
  42   TH1F *hCtsLikeSignNeg = new TH1F("hCtsLikeSignNeg","LikeSign cosine of decay angle for -- pairs",50,-1.,1.);
  43   hCtsLikeSignNeg->SetXTitle("cos #theta^{*}");
  44   hCtsLikeSignNeg->SetYTitle("Entries");
  45   hCtsLikeSignNeg->Sumw2();
  46   // Cosine of pointing angle
  47   TH1F *hCPtaJPSI = new TH1F("hCPtaJPSI","JPSI cosine of pointing angle distribution",100,-1.,1.);
  48   hCPtaJPSI->SetXTitle("cos #theta_{point}");
  49   hCPtaJPSI->SetYTitle("Entries");
  50   hCPtaJPSI->Sumw2();
  51   TH1F *hCPtaLikeSign = new TH1F("hCPtaLikeSign","LikeSign cosine of pointing angle distribution",100,-1.,1.);
  52   hCPtaLikeSign->SetXTitle("cos #theta_{point}");
  53   hCPtaLikeSign->SetYTitle("Entries");
  54   hCPtaLikeSign->Sumw2();
  55   // d0 x d0
  56   TH1F *hd0d0JPSI = new TH1F("hd0d0JPSI","JPSI Product of the impact parameters",100,-100000,100000);
  57   hd0d0JPSI->SetXTitle("d_{0}^{k} #times d_{0}^{#pi} [#mu m^{2}]");
  58   hd0d0JPSI->SetYTitle("Entries");
  59   hd0d0JPSI->Sumw2();
  60   TH1F *hd0d0LikeSign = new TH1F("hd0d0LikeSign","LikeSign Product of the impact parameters",100,-100000,100000);
  61   hd0d0LikeSign->SetXTitle("d_{0}^{+/-} #times d_{0}^{+/-} [#mu m^{2}]");
  62   hd0d0LikeSign->SetYTitle("Entries");
  63   hd0d0LikeSign->Sumw2();
  64   // DCA
  65   TH1F *hDCAJPSI = new TH1F("hDCAJPSI","JPSI distance of closest approach",100,0.,5.);
  66   hDCAJPSI->SetXTitle("dca [10^{2}#mu m]");
  67   hDCAJPSI->SetYTitle("Entries*10^{2}");
  68   hDCAJPSI->Sumw2();
  69   TH1F *hDCALikeSign = new TH1F("hDCALikeSign","LikeSign distance of closest approach",100,0.,5.);
  70   hDCALikeSign->SetXTitle("dca [10^{2}#mu m]");
  71   hDCALikeSign->SetYTitle("Entries*10^{2}");
  72   hDCALikeSign->Sumw2();
  73
  74   Int_t nTotHF=0,nTotLikeSign=0;
  75   //Double_t normalizationFactorAvg=0;
  76   Double_t weightLS=1.;
  77   Double_t nTotPosPairs=0;
  78   Double_t nTotNegPairs=0;
  79   //Double_t nTotPosPairsAvg=0;
  80   //Double_t nTotNegPairsAvg=0;//comment in case of one directory processing
  81   Int_t totEvents=0;
  82
  83   // open input file and get the TTree
  84   TFile inFile(aodFileName,"READ");
  85   if (!inFile.IsOpen()) return;
  86
  87   TTree *aodTree = (TTree*)inFile.Get("aodTree");
  88   aodTree->AddFriend("aodTree",aodHFFileName);
  89
  90   AliAODEvent *aod = new AliAODEvent();
  91
  92   aod->ReadFromTree(aodTree);
  93
  94   // load heavy flavour vertices
  95   TClonesArray *arrayVerticesHF =
  96     (TClonesArray*)aod->GetList()->FindObject("VerticesHF");
  97
  98   // load like sign candidates
  99   TClonesArray *arrayLikeSign =
 100     (TClonesArray*)aod->GetList()->FindObject("LikeSign2Prong");
 101
 102   // load B->JPSI->e+e- candidates
 103   TClonesArray *arrayBtoJpsiToEle =
 104     (TClonesArray*)aod->GetList()->FindObject("JPSItoEle");
 105
 106   Bool_t normalize = kTRUE;
 107
 108   Double_t cuts[9]=
 109   // cuts[0] = inv. mass half width [GeV]
 110   // cuts[1] = dca [cm]
 111   // cuts[2] = cosThetaStar (negative electron)
 112   // cuts[3] = pTP [GeV/c]
 113   // cuts[4] = pTN [GeV/c]
 114   // cuts[5] = d0P [cm]   upper limit!
 115   // cuts[6] = d0N [cm]  upper limit!
 116   // cuts[7] = d0d0 [cm^2]
 117   // cuts[8] = cosThetaPoint
 118                      {1000.,
 119                       100000.,
 120                       1.1,
 121                       0.,
 122                       0.,
 123                       100000.,
 124                       100000.,
 125                       100000000.,
 126                       -1.1};
 127   //Double_t cuts[9]= {1000., 100000., 0.8, 0., 0., 100000., 100000., 80000., 0.8};
 128
 129   Int_t nTotBtoJpsiToEle=0;
 130   Int_t nTotD0toKpi=0;
 131   AliAODVertex *vtx1=0;
 132
 133   // loop over events
 134   Int_t nEvents = aodTree->GetEntries();
 135   for (Int_t nEv = 0; nEv < nEvents; nEv++) {
 136     cout<<"\n------------ Event: "<<nEv<<" ------------------"<<endl;
 137
 138     // read event
 139     aodTree->GetEvent(nEv);
 140
 141     //print event info
 142     aod->GetHeader()->Print();
 143
 144     // primary vertex
 145     vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 146     vtx1->Print();
 147
 148     // make trkIDtoEntry register (temporary)
 149     Int_t trkIDtoEntry[100000];
 150     for(Int_t it=0;it<aod->GetNumberOfTracks();it++) {
 151       AliAODTrack *track = aod->GetTrack(it);
 152       trkIDtoEntry[track->GetID()]=it;
 153     }
 154
 155     // loop over Like sign candidates
 156     Int_t nLikeSign = arrayLikeSign->GetEntriesFast();
 157     nTotLikeSign += nLikeSign;
 158     cout<<"Number of LikeSign pairs: "<<nLikeSign<<endl;
 159     Int_t nPosPairs=0,nNegPairs=0;
 160
 161     for (Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
 162       AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
 163       Bool_t unsetvtx=kFALSE;
 164       if(!d->GetOwnPrimaryVtx()) {
 165         d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 166         unsetvtx=kTRUE;
 167       }
 168       Int_t okBtoJPSIls=0;
 169       if(d->SelectBtoJPSI(cuts,okBtoJPSIls)) {
 170         hMass->Fill(d->InvMassJPSIee());
 171         hCPtaVSd0d0->Fill(1e8*d->Prodd0d0(),d->CosPointingAngle());
 172         hSecVtxZ->Fill(d->GetSecVtxZ());
 173         hCPtaLikeSign->Fill(d->CosPointingAngle());
 174         hd0d0LikeSign->Fill(1e8*d->Prodd0d0());
 175         hCtsLikeSign->Fill(d->CosThetaStarJPSI());
 176         hDCALikeSign->Fill(100*d->GetDCA());
 177         AliAODTrack *trk0 = (AliAODTrack*)d->GetDaughter(0);
 178         if(!trk0) {
 179           trk0=aod->GetTrack(trkIDtoEntry[d->GetProngID(0)]);
 180           cout<<"references to standard AOD not available"<<endl;
 181         }
 182         if((trk0->Charge())==1) {
 183           nPosPairs++;
 184           hCtsLikeSignPos->Fill(d->CosThetaStarJPSI());
 185         } else {
 186           nNegPairs++;
 187           hCtsLikeSignNeg->Fill(d->CosThetaStarJPSI());
 188         }
 189       }
 190       if(unsetvtx) d->UnsetOwnPrimaryVtx();
 191     }
 192     cout<<"\n------------ N. of positive pairs in Event: "<<nEv<<" ----- "<< nPosPairs<<endl;
 193     cout<<"\n------------ N. of negative pairs in Event: "<<nEv<<" ----- "<< nNegPairs<<endl;
 194
 195     nTotPosPairs += nPosPairs;
 196     nTotNegPairs += nNegPairs;
 197
 198     // loop over JPSI candidates
 199     Int_t nBtoJpsiToEle = arrayBtoJpsiToEle->GetEntriesFast();
 200     nTotBtoJpsiToEle += nBtoJpsiToEle;
 201     cout<<"Number of B->J/Psi->e+e-: "<<nBtoJpsiToEle<<endl;
 202
 203     for (Int_t iBtoJpsiToEle = 0; iBtoJpsiToEle < nBtoJpsiToEle; iBtoJpsiToEle++) {
 204       AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)arrayBtoJpsiToEle->UncheckedAt(iBtoJpsiToEle);
 205       Bool_t unsetvtx=kFALSE;
 206       if(!d->GetOwnPrimaryVtx()) {
 207         d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 208         unsetvtx=kTRUE;
 209       }
 210       Int_t okBtoJPSI=0;
 211       if(d->SelectBtoJPSI(cuts,okBtoJPSI)) {
 212         hMassJPSI->Fill(d->InvMassJPSIee());
 213         hCtsJPSI->Fill(d->CosThetaStarJPSI());
 214         hd0d0JPSI->Fill(1e8*d->Prodd0d0());
 215         hCPtaJPSI->Fill(d->CosPointingAngle());
 216         hDCAJPSI->Fill(100*d->GetDCA());
 217       }
 218       if(unsetvtx) d->UnsetOwnPrimaryVtx();
 219     }
 220
 221     Int_t nJPSItoEle = arrayBtoJpsiToEle->GetEntriesFast();
 222     nTotBtoJpsiToEle += nJPSItoEle;
 223     cout<<"Number of JPSI->ee: "<<nJPSItoEle<<endl;
 224
 225   }
 226
 227   //nTotPosPairsAvg = nTotPosPairs/nEvents;
 228   //nTotNegPairsAvg = nTotNegPairs/nEvents;
 229   //normalizationFactorAvg = 2.*TMath::Sqrt(nTotPosPairsAvg*nTotNegPairsAvg);
 230   Double_t normalizationFactor = 2.*TMath::Sqrt(nTotPosPairs*nTotNegPairs);
 231
 232   //cout << "\n------------ Normalization factor averaged over events: -------- " << normalizationFactorAvg << endl;
 233
 234   if(normalize) weightLS = normalizationFactor;
 235
 236   cout <<"Applied weight to like-sign spectrum ------- " << weightLS << endl;
 237   printf("\n Total HF vertices: %d\n",nTotHF);
 238   printf("\n Total LikeSign pairs: %d\n",nTotLikeSign);
 239
 240   TCanvas *c1 = new TCanvas("c1","c1",0,0,1000,1000);
 241   c1->Divide(2,3);
 242   c1->cd(1);
 243   hCPtaLikeSign->SetMarkerStyle(20);
 244   hCPtaLikeSign->SetMarkerSize(0.7);
 245   hCPtaLikeSign->SetMarkerColor(4);
 246   hCPtaLikeSign->Scale((1/weightLS)*hCPtaLikeSign->GetEntries());
 247   hCPtaLikeSign->Draw();
 248   hCPtaJPSI->SetMarkerStyle(20);
 249   hCPtaJPSI->SetMarkerSize(0.7);
 250   hCPtaJPSI->SetMarkerColor(2);
 251   hCPtaJPSI->Draw("same");
 252   c1->cd(2)->SetLogy();
 253   hd0d0LikeSign->SetMarkerStyle(20);
 254   hd0d0LikeSign->SetMarkerSize(0.7);
 255   hd0d0LikeSign->SetMarkerColor(4);
 256   hd0d0LikeSign->Scale((1/weightLS)*hd0d0LikeSign->GetEntries());
 257   hd0d0LikeSign->Draw();
 258   hd0d0JPSI->SetMarkerStyle(20);
 259   hd0d0JPSI->SetMarkerSize(0.7);
 260   hd0d0JPSI->SetMarkerColor(2);
 261   hd0d0JPSI->Draw("same");
 262   c1->cd(3)->SetLogy();
 263   hCtsLikeSign->SetMarkerStyle(20);
 264   hCtsLikeSign->SetMarkerSize(0.7);
 265   hCtsLikeSign->SetMarkerColor(4); // all LS couples
 266   hCtsLikeSign->Scale((1/weightLS)*hCtsLikeSign->GetEntries());
 267   hCtsLikeSignPos->SetMarkerStyle(20);
 268   hCtsLikeSignPos->SetMarkerSize(0.6);
 269   hCtsLikeSignPos->SetMarkerColor(3); // ++
 270   hCtsLikeSignPos->Scale((1/weightLS)*hCtsLikeSignPos->GetEntries());
 271   hCtsLikeSignNeg->SetMarkerStyle(21);
 272   hCtsLikeSignNeg->SetMarkerSize(0.6);
 273   hCtsLikeSignNeg->SetMarkerColor(7); // --
 274   hCtsLikeSignNeg->Scale((1/weightLS)*hCtsLikeSignNeg->GetEntries());
 275   hCtsLikeSign->Draw();
 276   hCtsLikeSignPos->Draw("same");
 277   hCtsLikeSignNeg->Draw("same");
 278   hCtsJPSI->SetMarkerStyle(20);
 279   hCtsJPSI->SetMarkerSize(0.7);
 280   hCtsJPSI->SetMarkerColor(2); // JPSI is red
 281   hCtsJPSI->Draw("same");
 282   c1->cd(4);
 283   hDCALikeSign->SetMarkerStyle(20);
 284   hDCALikeSign->SetMarkerSize(0.7);
 285   hDCALikeSign->SetMarkerColor(4);
 286   hDCALikeSign->Scale((1/weightLS)*hDCALikeSign->GetEntries());
 287   hDCALikeSign->Draw();
 288   hDCAJPSI->SetMarkerStyle(20);
 289   hDCAJPSI->SetMarkerSize(0.7);
 290   hDCAJPSI->SetMarkerColor(2);
 291   hDCAJPSI->Draw("same");
 292   c1->cd(5)->SetLogy();
 293   hMass->SetMarkerStyle(20);
 294   hMass->SetMarkerSize(0.7);
 295   hMass->SetMarkerColor(4);
 296   hMass->Scale((1/weightLS)*hMass->GetEntries());
 297   hMass->Draw();
 298   hMassJPSI->SetMarkerStyle(20);
 299   hMassJPSI->SetMarkerSize(0.7);
 300   hMassJPSI->SetMarkerColor(2);
 301   hMassJPSI->Draw("same");
 302   c1->cd(6)->SetLogy();
 303   hMassJPSI->Draw();
 304
 305   return;
 306 }