PWGLF/RESONANCES/macros/utils/ResonanceModel.C

   1 /* pt generation limits */
   2 const Double_t ptMin = 0.;
   3 const Double_t ptMax = 10.;
   4 /* y generation limits */
   5 const Double_t yMin = -0.1;
   6 const Double_t yMax = 0.1;
   7 /* phi generation limits */
   8 const Double_t phiMin = -TMath::Pi();
   9 const Double_t phiMax = TMath::Pi();
  10
  11 const Double_t pionMass = 139.57018e-3;
  12 const Double_t kaonMass = 493.677e-3;
  13 const Double_t protonMass = 938.272046e-3;
  14
  15 PhiModel(Int_t ngen = 1.e6)
  16 {
  17
  18     Double_t m = TDatabasePDG::Instance()->GetParticle("phi")->Mass();
  19     Double_t w = TDatabasePDG::Instance()->GetParticle("phi")->Width();
  20     Double_t dm[2] = {kaonMass, kaonMass};
  21     Int_t di[2] = {3, -3};
  22
  23     ResonanceModel(m, w, 2, dm, di, "PhiModel.root", ngen);
  24 }
  25
  26 KStarModel(Int_t ngen = 1.e6, Bool_t antiparticle = kFALSE)
  27 {
  28
  29     Double_t m = TDatabasePDG::Instance()->GetParticle("K*0")->Mass();
  30     Double_t w = TDatabasePDG::Instance()->GetParticle("K*0")->Width();
  31     Double_t dm[2] = {kaonMass, pionMass};
  32     Int_t di[2] = {3, -2};
  33     if (antiparticle) {di[0] = -di[0]; di[1] = -di[1];}
  34
  35     if (antiparticle) ResonanceModel(m, w, 2, dm, di, "KStarBarModel.root", ngen);
  36     else ResonanceModel(m, w, 2, dm, di, "KStarModel.root", ngen);
  37 }
  38
  39 Lambda1520Model(Int_t ngen = 1.e6, Bool_t antiparticle = kFALSE)
  40 {
  41
  42     Double_t m = 1519.54e-3; /* GeV */
  43     Double_t w = 15.73e-3; /* GeV */
  44     Double_t dm[2] = {protonMass, kaonMass};
  45     Int_t di[2] = {4, -3};
  46     if (antiparticle) {di[0] = -di[0]; di[1] = -di[1];}
  47
  48     if (antiparticle) ResonanceModel(m, w, 2, dm, di, "Lambda1520BarModel.root", ngen);
  49     else ResonanceModel(m, w, 2, dm, di, "Lambda1520Model.root", ngen);
  50 }
  51
  52 ResonanceModel(Double_t resonanceMass, Double_t resonanceWidth, Int_t nDaughters, Double_t *daughterMass, Int_t *daughterID, const Char_t *outname, Int_t ngen = 1.e6)
  53 {
  54
  55     /* functions */
  56     TF1 *fBW = new TF1("fBW", "[0] * TMath::BreitWigner(x, [1], [2])");
  57     fBW->SetParameter(0, ngen / 100.);
  58     fBW->SetParameter(1, resonanceMass);
  59     fBW->SetParameter(2, resonanceWidth);
  60
  61     /* histograms */
  62     Int_t massBins = 50;
  63     Int_t ptBins = 50;
  64     Double_t massMin = resonanceMass - 10. * resonanceWidth;
  65     Double_t massMax = resonanceMass + 10. * resonanceWidth;
  66
  67     TH1F *hGenMass = new TH1F("hGenMass", ";m (GeV/#it{c}^{2});", massBins, massMin, massMax);
  68     hGenMass->SetMarkerColor(2);
  69     hGenMass->SetMarkerStyle(25);
  70     hGenMass->Sumw2();
  71
  72     TH1F *hRecMass = new TH1F("hRecMass", ";m (GeV/#it{c}^{2});", massBins, massMin, massMax);
  73     hRecMass->SetMarkerColor(4);
  74     hRecMass->SetMarkerStyle(20);
  75     hRecMass->Sumw2();
  76
  77     TH1F *hGenPt = new TH1F("hGenPt", ";#it{p}_{T} (GeV/#it{c});", ptBins, ptMin, ptMax);
  78     hGenPt->SetMarkerColor(2);
  79     hGenPt->SetMarkerStyle(25);
  80     hGenPt->Sumw2();
  81
  82     TH1F *hRecPt = new TH1F("hRecPt", ";#it{p}_{T} (GeV/#it{c});", ptBins, ptMin, ptMax);
  83     hRecPt->SetMarkerColor(4);
  84     hRecPt->SetMarkerStyle(20);
  85     hRecPt->Sumw2();
  86
  87     /* init random generator */
  88     TRandom *gRandom = new TRandom();
  89     gRandom->SetSeed(123456789);
  90
  91     /* resonance four-momentum */
  92     TLorentzVector resonance;
  93     TLorentzVector recoP;
  94     TGenPhaseSpace decay;
  95
  96     /* benchmark */
  97     TBenchmark bm;
  98     bm.Start("time");
  99
 100     /* loop over generation */
 101     Int_t igen = 0;
 102     while (igen < ngen) {
 103
 104         /* generate mass */
 105         Double_t mass = gRandom->BreitWigner(resonanceMass, resonanceWidth);
 106         /* generate pt */
 107         Double_t pt = gRandom->Uniform(ptMin, ptMax);
 108         /* generate y */
 109         Double_t y = gRandom->Uniform(yMin, yMax);
 110         /* get eta */
 111         Double_t eta = y2eta(pt, resonanceMass, y);
 112         /* generate phi */
 113         Double_t phi = gRandom->Uniform(phiMin, phiMax);
 114
 115         /* init resonance four-momentum */
 116         resonance.SetPtEtaPhiM(pt, eta, phi, mass);
 117         /* decay particle, if allowed */
 118         if (!decay.SetDecay(resonance, nDaughters, daughterMass)) continue;
 119
 120         /* decay allowed */
 121         igen++;
 122         decay.Generate();
 123
 124         /* fill gen plots */
 125         hGenMass->Fill(mass);
 126         hGenPt->Fill(pt);
 127
 128         /* reconstructed invariant mass */
 129         recoP.SetPtEtaPhiM(0., 0., 0., 0.);
 130         TLorentzVector *daughterP;
 131         Bool_t daughterLost = kFALSE;
 132         for (Int_t idaugh = 0; idaugh < nDaughters; idaugh++) {
 133             daughterP = decay.GetDecay(idaugh);
 134             recoP += *daughterP;
 135             Double_t eff = GetEfficiency(daughterP, daughterID[idaugh]);
 136             if (gRandom->Uniform(0., 1.) > eff) daughterLost = kTRUE;
 137         }
 138
 139         /* fill reco plots */
 140         if (!daughterLost) {
 141             hRecMass->Fill(recoP.M());
 142             hRecPt->Fill(pt);
 143         }
 144     }
 145
 146     /* benchmark */
 147     bm.Stop("time");
 148     bm.Print("time");
 149
 150     gStyle->SetHistMinimumZero();
 151
 152     TCanvas *cGenRec = new TCanvas("cGenRec", "cRecGen", 600, 600);
 153     cGenRec->Divide(1, 2);
 154     cGenRec->cd(1);
 155     hGenMass->DrawCopy();
 156     hRecMass->DrawCopy("same");
 157     cGenRec->cd(2);
 158     hGenPt->DrawCopy();
 159     hRecPt->DrawCopy("same");
 160
 161     TCanvas *cEff = new TCanvas("cEff", "cEff", 600, 600);
 162     cEff->Divide(1, 2);
 163     cEff->cd(1);
 164     TH1 *hEffMass = (TH1 *)hRecMass->Clone("hEffMass");
 165     hEffMass->Divide(hEffMass, hGenMass, 1., 1., "B");
 166     hEffMass->DrawCopy();
 167     cEff->cd(2);
 168     TH1 *hEffPt = (TH1 *)hRecPt->Clone("hEffPt");
 169     hEffPt->Divide(hEffPt, hGenPt, 1., 1., "B");
 170     hEffPt->DrawCopy();
 171
 172     TFile *fout = TFile::Open(outname, "RECREATE");
 173     hGenMass->Write();
 174     hRecMass->Write();
 175     hGenPt->Write();
 176     hRecPt->Write();
 177     hEffPt->Write();
 178     fout->Close();
 179
 180 }
 181
 182 /*****************************************************************/
 183
 184 TFile *fTrackEff = NULL;
 185 TFile *fMatchEff = NULL;
 186
 187 Double_t
 188 GetEfficiency(TLorentzVector *lv, Int_t ipart)
 189 {
 190
 191     Char_t *pname[5] = {"", "", "pion", "kaon", "proton"};
 192     Char_t *cname[2] = {"positive", "negative"};
 193
 194     Int_t icharge = ipart > 0 ? 0 : 1;
 195     ipart = TMath::Abs(ipart);
 196
 197     Double_t pt = lv->Pt();
 198     Double_t eta = lv->Eta();
 199     Double_t phi = lv->Phi();
 200
 201     if (!fTrackEff) fTrackEff = TFile::Open("TOF_trackingEfficiency.root");
 202     if (!fMatchEff) fMatchEff = TFile::Open("TOF_matchingEfficiency.root");
 203     TH1 *hTrackEff = fTrackEff->Get(Form("hTrackingEff_MB_%s_%s", pname[ipart], cname[icharge]));
 204     TH1 *hMatchEff = fMatchEff->Get(Form("hMatchEff_MB_%s_%s", pname[ipart], cname[icharge]));
 205
 206     Double_t trackEff = hTrackEff->Interpolate(pt);
 207     Double_t matchEff = hMatchEff->Interpolate(pt);
 208
 209     Double_t trackEff_GF = TrackingEff_geantflukaCorrection(ipart, icharge)->Eval(pt);
 210     Double_t matchEff_GF = TOFmatchMC_geantflukaCorrection(ipart, icharge)->Eval(pt);
 211
 212     Double_t eff = trackEff;
 213     if (ipart == 3 && pt > 0.6) eff = trackEff * matchEff;
 214     if (ipart == 4 && pt > 1.1) eff = trackEff * matchEff;
 215
 216     Double_t pidEff = 0.954;
 217     eff *= pidEff;
 218
 219 //    return trackEff;
 220     return eff;
 221 }
 222
 223 /*****************************************************************/
 224
 225 Double_t
 226 y2eta(Double_t pt, Double_t mass, Double_t y){
 227     Double_t mt = TMath::Sqrt(mass * mass + pt * pt);
 228     return TMath::ASinH(mt / pt * TMath::SinH(y));
 229 }
 230 Double_t
 231 eta2y(Double_t pt, Double_t mass, Double_t eta){
 232     Double_t mt = TMath::Sqrt(mass * mass + pt * pt);
 233     return TMath::ASinH(pt / mt * TMath::SinH(eta));
 234 }
 235
 236 /*****************************************************************/
 237
 238 TF1 *fTrackEff_GF[5][2] = {NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL};
 239 TF1 *fMatchEff_GF[5][2] = {NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL};
 240
 241 TF1 *
 242 TrackingEff_geantflukaCorrection(Int_t ipart, Int_t icharge)
 243 {
 244
 245     if (fTrackEff_GF[ipart][icharge]) return fTrackEff_GF[ipart][icharge];
 246
 247     Char_t *pname[3] = {"pion", "kaon", "proton"};
 248     Char_t *cname[2] = {"positive", "negative"};
 249
 250     TF1 *f;
 251     if (ipart == 3 && icharge == 1)
 252         f = new TF1(Form("fGeantFluka_%s_%s", pname[ipart - 2], cname[icharge]), "TrackingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.);
 253     else if (ipart == 4 && icharge == 1)
 254         f = new TF1(Form("fGeantFluka_%s_%s", pname[ipart - 2], cname[icharge]), "TrackingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.);
 255     else
 256         f = new TF1(Form("fGeantFluka_%s_%s", pname[ipart - 2], cname[icharge]), "TrackingPtGeantFlukaCorrectionNull(x)", 0., 5.);
 257
 258     fTrackEff_GF[ipart][icharge] = f;
 259     return f;
 260 }
 261
 262 Double_t
 263 TrackingPtGeantFlukaCorrectionNull(Double_t pTmc)
 264 {
 265     return 1.;
 266 }
 267
 268 Double_t
 269 TrackingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
 270 {
 271     return (1 - 0.129758 *TMath::Exp(-pTmc*0.679612));
 272 }
 273
 274 Double_t
 275 TrackingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
 276 {
 277     return TMath::Min((0.972865 + 0.0117093*pTmc), 1.);
 278 }
 279
 280 /*****************************************************************/
 281
 282 TF1 *
 283 TOFmatchMC_geantflukaCorrection(Int_t ipart, Int_t icharge)
 284 {
 285
 286     if (fMatchEff_GF[ipart][icharge]) return fMatchEff_GF[ipart][icharge];
 287
 288     Char_t *pname[3] = {"pion", "kaon", "proton"};
 289     Char_t *cname[2] = {"positive", "negative"};
 290
 291     TF1 *f;
 292     if (ipart == 3 && icharge == 1)
 293         f = new TF1(Form("fGeantFluka_%s_%s", pname[ipart - 2], cname[icharge]), "MatchingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.);
 294     else if (ipart == 4 && icharge == 1)
 295         f = new TF1(Form("fGeantFluka_%s_%s", pname[ipart - 2], cname[icharge]), "MatchingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.);
 296     else
 297         f = new TF1(Form("fGeantFluka_%s_%s", pname[ipart - 2], cname[icharge]), "MatchingPtGeantFlukaCorrectionNull(x)", 0., 5.);
 298
 299     fMatchEff_GF[ipart][icharge] = f;
 300     return f;
 301 }
 302
 303 Double_t
 304 MatchingPtGeantFlukaCorrectionNull(Double_t pTmc)
 305 {
 306     return 1.;
 307 }
 308
 309 Double_t
 310 MatchingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
 311 {
 312     Float_t ptTPCoutP =pTmc*(1-6.81059e-01*TMath::Exp(-pTmc*4.20094));
 313     return (TMath::Power(1 - 0.129758*TMath::Exp(-ptTPCoutP*0.679612),0.07162/0.03471));
 314 }
 315
 316 Double_t
 317 MatchingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
 318 {
 319     Float_t ptTPCoutK=pTmc*(1- 3.37297e-03/pTmc/pTmc - 3.26544e-03/pTmc);
 320     return TMath::Min((TMath::Power(0.972865 + 0.0117093*ptTPCoutK,0.07162/0.03471)), 1.);
 321 }