PWG3/vertexingOld/AliD0toKpiPlots.C

   1 void AliD0toKpiPlots(const Char_t *inName="AliD0toKpi.root",
   2                      const Char_t *outName="D0histograms.root") {
   3   //--------------------------------------------------------------------------
   4   // This macro histograms many variables of D0->Kpi candidates
   5   //
   6   //     Andrea Dainese, andrea.dainese@lnl.infn.it
   7   //--------------------------------------------------------------------------
   8
   9   gSystem->Load("libANALYSIS");
  10   gSystem->Load("libAOD.so");
  11   gSystem->Load("libPWG3base.so");
  12   gSystem->Load("libPWG3vertexingOld.so");
  13
  14   // set of cuts
  15   Double_t D0Cuts[9] = {0.1,         // mass [GeV]
  16                         1000000.,          // dca [micron]
  17                         1.1,           // cosThetaStar
  18                         0.,           // pT K [GeV/c]
  19                         0.,           // pT Pi [GeV/c]
  20                         100000.,       // d0K upper [micron]
  21                         100000.,       // d0Pi upper [micron]
  22                         10000000000.,            // d0d0 [micron^2]
  23                         -1.1};          // cosThetaPointing
  24
  25   // number of events (for normalization)
  26   Bool_t normalize = kFALSE;
  27   Double_t  events = 1.;
  28
  29
  30
  31   // define histograms
  32   TH1F *hptK = new TH1F("hptK","\"K\" p_{t} distribution",50,0,10);
  33   hptK->SetXTitle("p_{t} [GeV]");
  34
  35   TH1F *hptPi = new TH1F("hptPi","\"#pi\" p_{t} distribution",50,0,10);
  36   hptPi->SetXTitle("p_{t} [GeV]");
  37
  38   TH1F *hDCA = new TH1F("hDCA","DCA",50,0,1000);
  39   hDCA->SetXTitle("dca [#mu m]");
  40
  41   TH1F *hptD0 = new TH1F("hptD0","D^{0} p_{t} distribution",40,0,40);
  42   hptD0->SetXTitle("p_{t} [GeV]");
  43
  44   TH1F *hyD0 = new TH1F("hyD0","D^{0} rapidity distribution",50,-2,2);
  45   hyD0->SetXTitle("y");
  46
  47   TH1F *hCPtaD0 = new TH1F("hCPtaD0","cosine of pointing angle distribution",100,-1,1);
  48   hCPtaD0->SetXTitle("cos #theta_{point}");
  49
  50   TH1F *hCPtaXY = new TH1F("hCPtaXY","cosine of pointing angle in (x,y) plane",100,-1,1);
  51   hCPtaXY->SetXTitle("cos #theta_{point}");
  52
  53   TH1F *hCts = new TH1F("hCts","cosine of decay angle",50,-1.2,1.2);
  54   hCts->SetXTitle("cos #theta^{*}");
  55
  56   TH2F *hCtsVsPtK = new TH2F("hCtsVsPtK","cosine of decay angle VS \"K\" p_{t}",50,0,5,50,-1,1);
  57   hCtsVsPtK->SetYTitle("cos #theta^{*}");
  58   hCtsVsPtK->SetXTitle("p_{t} [GeV]");
  59
  60   TH1F *hd0d0 = new TH1F("hd0d0","Product of the impact parameters",100,-100000,100000);
  61   hd0d0->SetXTitle("d_{0}^{K} #times d_{0}^{#pi} [#mu m^{2}]");
  62
  63   TH1F *hd0K = new TH1F("hd0K","Impact parameter of \"K\"",100,-5000,5000);
  64   hd0K->SetXTitle("d_{0}^{K} [#mu m]");
  65
  66   TH1F *hd0Pi = new TH1F("hd0Pi","Impact parameter of \"#pi\"",100,-5000,5000);
  67   hd0Pi->SetXTitle("d_{0}^{#pi} [#mu m]");
  68
  69   TH2F *hCPtaVsd0d0 = new TH2F("hCPtaVsd0d0","cos #theta_{point} vs d_{0}^{K} #times d_{0}^{#pi}",100,-100000,100000,100,-1,1);
  70   hCPtaVsd0d0->SetXTitle("d_{0}^{K} #times d_{0}^{#pi} [#mu m^{2}]");
  71   hCPtaVsd0d0->SetYTitle("cos #theta_{point}");
  72
  73   TH2F *hCPtaVsd0d0zoom = new TH2F("hCPtaVsd0d0zoom","cos #theta_{point} vs d_{0}^{K} #times d_{0}^{#pi}",100,-100000,0,100,.9,1);
  74   hCPtaVsd0d0zoom->SetXTitle("d_{0}^{K} #times d_{0}^{#pi} [#mu m^{2}]");
  75   hCPtaVsd0d0zoom->SetYTitle("cos #theta_{point}");
  76
  77   TH2F *hd0d0VSptD0 = new TH2F("hd0d0VSptD0","d_{0}^{K} #times d_{0}^{#pi} VS D^{0} p_{t}",50,0,25,100,-120000,120000);
  78   hd0d0VSptD0->SetYTitle("d_{0}^{K} #times d_{0}^{#pi} [#mu m^{2}]");
  79   hd0d0VSptD0->SetXTitle("D^{0} p_{t} [GeV]");
  80
  81   TH1F *hMass = new TH1F("hMass","Invariant mass distribution",50,1.765,1.965);
  82   hMass->SetXTitle("M[K,#pi] [GeV]");
  83
  84   TH2F *hArm = new TH2F("hArm","Armenteros plot",50,-2,2,50,0,1);
  85   hArm->SetXTitle("#alpha");
  86   hArm->SetYTitle("q_{t}");
  87
  88   // open input file and get tree
  89   TFile *inFile = TFile::Open(inName);
  90
  91   TTree *treeD0 = (TTree*)inFile->Get("TreeD0");
  92   AliD0toKpi *D = 0;
  93   treeD0->SetBranchAddress("D0toKpi",&D);
  94   Int_t entries = (Int_t)treeD0->GetEntries();
  95
  96   printf("+++\n+++ Number of D0 in tree:  %d\n+++\n",entries);
  97
  98   Double_t MD0,MD0bar,ctsD0,ctsD0bar,ctsPiD0,ctsPiD0bar;
  99   Double_t WgtD0,WgtD0bar;
 100   Double_t sampleABC=0.;
 101   Int_t okD0=0,okD0bar=0;
 102   Int_t nSel = 0;
 103   Int_t ptbin;
 104
 105   // loop on D0
 106   for(Int_t i=0; i<entries; i++) {
 107     if(i%10000==0) printf(" candidate %d of %d\n",i,entries);
 108
 109     // get event from tree
 110     treeD0->GetEvent(i);
 111     //--- select the PID strategy & compute weights
 112     //    D->ApplyPID("TOFparam_PbPb");
 113     //    D->ComputeWgts();
 114     // get weights for the three samples A+B+C
 115     //    D->GetWgts(WgtD0,WgtD0bar,"ABC");
 116     WgtD0 = 1.; WgtD0bar = 1.;
 117
 118     // normalize to 1 event
 119     if(normalize) { WgtD0 /= events; WgtD0bar /= events; }
 120
 121     // check if candidate passes selection (as D0 or D0bar)
 122     D->Select(D0Cuts,okD0,okD0bar);
 123
 124     // set weights to 0 if the candidate doesn't pass selection
 125     if(!okD0)    WgtD0=0.;
 126     if(!okD0bar) WgtD0bar=0.;
 127     if(okD0 || okD0bar) nSel++;
 128
 129     // count selected candidates
 130     sampleABC += WgtD0 + WgtD0bar;
 131
 132     // inv mass and cosThetaStar
 133     D->InvMass(MD0,MD0bar);
 134     D->CosThetaStar(ctsD0,ctsD0bar);
 135
 136     // fill histograms
 137     hptK->Fill(D->PtChild(1),WgtD0);
 138     hptK->Fill(D->PtChild(0),WgtD0bar);
 139     hptPi->Fill(D->PtChild(0),WgtD0);
 140     hptPi->Fill(D->PtChild(1),WgtD0bar);
 141     hd0K->Fill(D->Getd0Child(1),WgtD0);
 142     hd0K->Fill(D->Getd0Child(0),WgtD0bar);
 143     hd0Pi->Fill(D->Getd0Child(0),WgtD0);
 144     hd0Pi->Fill(D->Getd0Child(1),WgtD0bar);
 145     hMass->Fill(MD0,WgtD0);
 146     hMass->Fill(MD0bar,WgtD0bar);
 147     hCts->Fill(ctsD0,WgtD0);
 148     hCts->Fill(ctsD0bar,WgtD0bar);
 149     hCtsVsPtK->Fill(D->PtChild(1),ctsD0,WgtD0);
 150     hCtsVsPtK->Fill(D->PtChild(0),ctsD0bar,WgtD0bar);
 151     hDCA->Fill(D->GetDCA(),WgtD0+WgtD0bar);
 152     hptD0->Fill(D->Pt(),WgtD0+WgtD0bar);
 153     hyD0->Fill(D->Rapidity(),WgtD0+WgtD0bar);
 154     hd0d0->Fill(D->ProdImpParams(),WgtD0+WgtD0bar);
 155     hCPtaD0->Fill(D->CosPointing(),WgtD0+WgtD0bar);
 156     hCPtaXY->Fill(D->CosPointingXY(),WgtD0+WgtD0bar);
 157     hCPtaVsd0d0->Fill(D->ProdImpParams(),D->CosPointing(),WgtD0+WgtD0bar);
 158     hd0d0VSptD0->Fill(D->Pt(),D->ProdImpParams(),WgtD0+WgtD0bar);
 159     hCPtaVsd0d0zoom->Fill(D->ProdImpParams(),D->CosPointing(),WgtD0+WgtD0bar);
 160     hArm->Fill(D->Alpha(),D->Qt(),WgtD0+WgtD0bar);
 161
 162
 163   } // end loop on D0 candidates
 164
 165   inFile->Close();
 166
 167   printf("\n\n --- Total number of candidates passing selection: %d\n\n --- Sum of weights sample A+B+C: %f\n\n",nSel,sampleABC);
 168
 169   // draw histograms
 170   TCanvas *c1 = new TCanvas("c1","pt K & pi",0,0,700,700);
 171   c1->SetLogy();
 172   hptK->Draw();
 173   hptPi->Draw("same");
 174
 175   TCanvas *c2 = new TCanvas("c2","pt D0",0,0,700,700);
 176   c2->SetLogy();
 177   hptD0->Draw();
 178
 179   TCanvas *c3 = new TCanvas("c3","rapidity D0",0,0,700,700);
 180   hyD0->Draw();
 181
 182   TCanvas *c4 = new TCanvas("c4","pointing angle",0,0,700,700);
 183   hCPtaD0->Draw();
 184
 185   TCanvas *c5 = new TCanvas("c5","d0 x d0",0,0,700,700);
 186   c5->SetLogy();
 187   hd0d0->Draw();
 188
 189   TCanvas *c6 = new TCanvas("c6","pointing angle VS d0d0",0,0,700,700);
 190   c6->SetLogz();
 191   hCPtaVsd0d0->Draw("box");
 192
 193   TCanvas *c7 = new TCanvas("c7","mass",0,0,700,700);
 194   hMass->Draw();
 195
 196   TCanvas *c8 = new TCanvas("c8","armenteros",0,0,700,700);
 197   hArm->Draw("box");
 198
 199   TCanvas *c9 = new TCanvas("c9","decay angle",0,0,700,700);
 200   hCts->Draw();
 201
 202   TCanvas *c10 = new TCanvas("c10","dca",0,0,700,700);
 203   c10->SetLogy();
 204   hDCA->Draw();
 205
 206   TCanvas *c11 = new TCanvas("c11","d0 K & pi",0,0,700,700);
 207   c11->SetLogy();
 208   hd0K->Draw();
 209   hd0Pi->Draw("same");
 210
 211   // write all histograms to file
 212   TFile *outFile = new TFile(outName,"recreate");
 213   hMass->Write();
 214   hDCA->Write();
 215   hCts->Write();
 216   hCtsVsPtK->Write();
 217   hArm->Write();
 218   hCPtaVsd0d0->Write();
 219   hd0d0VSptD0->Write();
 220   hCPtaVsd0d0zoom->Write();
 221   hd0d0->Write();
 222   hCPtaD0->Write();
 223   hCPtaXY->Write();
 224   hptK->Write();
 225   hptPi->Write();
 226   hptD0->Write();
 227   hyD0->Write();
 228   hd0K->Write();
 229   hd0Pi->Write();
 230   outFile->Close();
 231
 232   return;
 233 }
 234
 235
 236
 237