PWGLF/SPECTRA/Nuclei/B2/AliLnBA.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 // coalescence parameter
  17 // author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
  18
  19 #include <TMath.h>
  20 #include <TFile.h>
  21 #include <TString.h>
  22 #include <TGraphErrors.h>
  23 #include <TF1.h>
  24
  25 #include "AliLnBA.h"
  26 #include "B2.h"
  27
  28 ClassImp(AliLnBA)
  29
  30 AliLnBA::AliLnBA(const TString& protonSpectra, const TString& protonTag, const TString& nucleusSpectra, const TString& nucleusTag, const TString& outputFilename, const TString& otag, Int_t a, Int_t z)
  31 : TObject()
  32 , fProtonSpectra(protonSpectra)
  33 , fProtonTag(protonTag)
  34 , fNucleusSpectra(nucleusSpectra)
  35 , fNucleusTag(nucleusTag)
  36 , fOutputFilename(outputFilename)
  37 , fOutputTag(otag)
  38 , fA(2)
  39 , fZ(1)
  40 , fNucleusName("Deuteron")
  41 , fCd(0.15)
  42 {
  43 //
  44 // constructor
  45 //
  46         this->SetNucleus(a,z);
  47 }
  48
  49 AliLnBA::~AliLnBA()
  50 {
  51 //
  52 // destructor
  53 //
  54 }
  55
  56 void AliLnBA::SetNucleus(Int_t a, Int_t z)
  57 {
  58 //
  59 // set nucleus mass and name
  60 //
  61         fA = a;
  62         fZ = z;
  63
  64         Int_t zz = TMath::Abs(z);
  65
  66         if(a==1 && zz==1)      fNucleusName = "Proton";
  67         else if(a==2 && zz==1) fNucleusName = "Deuteron";
  68         else if(a==3 && zz==1) fNucleusName = "Triton";
  69         else if(a==3 && zz==2) fNucleusName = "He3";
  70         else if(a==4 && zz==2) fNucleusName = "Alpha";
  71         else
  72         {
  73                 this->Warning("SetNucleus", "unknown nucleus A = %d, Z = %d", a, z);
  74                 fNucleusName = "Unknown";
  75         }
  76
  77         if(z<0) fNucleusName.Prepend("Anti");
  78 }
  79
  80 Int_t AliLnBA::Run()
  81 {
  82 //
  83 // coalescence parameter
  84 //
  85         TFile* finput1 = new TFile(fProtonSpectra.Data(), "read");
  86         if (finput1->IsZombie()) exit(1);
  87
  88         TFile* finputA = new TFile(fNucleusSpectra.Data(), "read");
  89         if (finputA->IsZombie()) exit(1);
  90
  91         // invariant differential yields
  92
  93         TString nucleonName = (fZ > 0) ? "Proton" : "AntiProton";
  94
  95         TGraphErrors* grPrtInvDYieldPt = FindObj<TGraphErrors>(finput1, fProtonTag, nucleonName + "_InvDiffYield_Pt");
  96         TGraphErrors* grNucInvDYieldPt = FindObj<TGraphErrors>(finputA, fNucleusTag, fNucleusName + "_InvDiffYield_Pt");
  97
  98         TFile* foutput = new TFile(fOutputFilename.Data(),"recreate");
  99
 100         foutput->mkdir(fOutputTag.Data());
 101         foutput->cd(fOutputTag.Data());
 102
 103         // coalescence parameter
 104
 105         TGraphErrors* grBAPt = this->GetBAPt(grPrtInvDYieldPt, grNucInvDYieldPt, fNucleusName + Form("_B%d_Pt",fA));
 106
 107         grBAPt->Write();
 108
 109         // homogeneity volume
 110         TGraphErrors* grR3Pt = 0;
 111         if(fA==2)
 112         {
 113                 grR3Pt = this->Rside2Rlong(grBAPt, fNucleusName + "_R3_Pt", fCd);
 114                 grR3Pt->Write();
 115         }
 116
 117         // propagate systematic errors if possible
 118
 119         TString grPrtName = fProtonTag +  "/" + nucleonName + "_SystErr_InvDiffYield_Pt;1";
 120         TString grNucName = fNucleusTag + "/" + fNucleusName + "_SystErr_InvDiffYield_Pt;1";
 121
 122         TGraphErrors* grSystErrPrtInvDYieldPt = dynamic_cast<TGraphErrors*>(finput1->Get(grPrtName.Data()));
 123         TGraphErrors* grSystErrNucInvDYieldPt = dynamic_cast<TGraphErrors*>(finputA->Get(grNucName.Data()));
 124
 125         if( (grSystErrPrtInvDYieldPt != 0) && (grSystErrNucInvDYieldPt != 0) )
 126         {
 127                 TGraphErrors* grSystErrBAPt = this->GetBAPt(grSystErrPrtInvDYieldPt, grSystErrNucInvDYieldPt, fNucleusName + Form("_SystErr_B%d_Pt",fA));
 128
 129                 if(fA==2)
 130                 {
 131                         TGraphErrors* grSystErrR3Pt = this->Rside2Rlong(grSystErrBAPt, fNucleusName + "_SystErr_R3_Pt", fCd);
 132                         grSystErrR3Pt->Write();
 133
 134                         delete grSystErrR3Pt;
 135                 }
 136
 137                 grSystErrBAPt->Write();
 138                 delete grSystErrBAPt;
 139         }
 140         else
 141         {
 142                 this->Warning("Run", "systematic errors are not propagated");
 143         }
 144
 145         delete grBAPt;
 146         delete grR3Pt;
 147
 148         delete foutput;
 149         delete finput1;
 150         delete finputA;
 151
 152         return 0;
 153 }
 154
 155 TGraphErrors* AliLnBA::GetBAPt(const TGraphErrors* grPrtInvDYieldPt, const TGraphErrors* grNucInvDYieldPt, const TString& name) const
 156 {
 157 //
 158 // coalescence parameter
 159 //
 160         TGraphErrors* grBAPt = new TGraphErrors();
 161         grBAPt->SetName(name.Data());
 162
 163         for(Int_t i=0, j=0; i < grNucInvDYieldPt->GetN(); ++i)
 164         {
 165                 Double_t ptNuc, yNuc;
 166
 167                 grNucInvDYieldPt->GetPoint(i, ptNuc, yNuc);
 168
 169                 if(ptNuc/fA < 0.4) continue; // acceptance
 170
 171                 Double_t yPrt = grPrtInvDYieldPt->Eval(ptNuc/fA); // interpolate
 172
 173                 if(yPrt == 0 || yNuc == 0 ) continue;
 174
 175                 Double_t bA = yNuc/TMath::Power(yPrt,fA);
 176
 177                 // error
 178                 Double_t ePrt = this->GetErrorY(grPrtInvDYieldPt, ptNuc/fA);
 179                 Double_t eNuc = grNucInvDYieldPt->GetErrorY(i);
 180
 181                 Double_t errPt = grNucInvDYieldPt->GetErrorX(i)/fA;
 182                 Double_t errBA = bA*TMath::Sqrt(TMath::Power(eNuc/yNuc,2) + TMath::Power(fA*ePrt/yPrt,2));
 183
 184                 grBAPt->SetPoint(j, ptNuc/fA, bA);
 185                 grBAPt->SetPointError(j++, errPt, errBA);
 186         }
 187
 188         return grBAPt;
 189 }
 190
 191 Double_t AliLnBA::GetErrorY(const TGraphErrors* gr, Double_t x0) const
 192 {
 193 //
 194 // estimate error of gr(x0) with the closest point to x0
 195 //
 196         const Double_t kEpsilon  = 1.e-6;
 197         const Double_t kUnknownError = 1.e+6;
 198
 199         for(Int_t i=0; i<gr->GetN(); ++i)
 200         {
 201                 Double_t x = gr->GetX()[i];
 202
 203                 if( TMath::Abs(x-x0) < kEpsilon) return gr->GetErrorY(i);
 204
 205                 if( ((i == 0) && (x > x0)) || ((i == (gr->GetN()-1)) && (x < x0)) )
 206                 {
 207                         this->Warning("GetErrorY", "%f is out of bounds",x);
 208                         return kUnknownError;
 209                 }
 210
 211                 if( x > x0 )
 212                 {
 213                         this->Warning("GetErrorY", "Interpolating error at %f",x0);
 214                         return (gr->GetErrorY(i)+gr->GetErrorY(i-1))/2;
 215                 }
 216         }
 217
 218         return 0;
 219 }
 220
 221 Double_t AliLnBA::Rside2Rlong(Double_t pt, Double_t B2, Double_t Cd) const
 222 {
 223 //
 224 // Rside^2*Rlong from B2 value
 225 // Phys. Rev. C, vol. 59, no. 3, pp. 1585--1602, 1999
 226 // (for pp ignore the exponential term)
 227 //
 228         Double_t hbarc = 0.1973269631; // GeV*fm
 229         Double_t m = 0.938272013;  // GeV/c^2
 230         Double_t pi = TMath::Pi();
 231
 232         if(B2==0) return 0;
 233
 234         Double_t mt = TMath::Sqrt(pt*pt + m*m);
 235         Double_t r3 = 3.*TMath::Power(pi,3./2.)*TMath::Power(hbarc,3.)*Cd/(2.*mt*B2);
 236
 237         return r3;
 238 }
 239
 240 TGraphErrors* AliLnBA::Rside2Rlong(const TGraphErrors* grB2, const TString& name, Double_t Cd) const
 241 {
 242 //
 243 // Rside^2*Rlong from B2 value
 244 // Phys. Rev. C, vol. 59, no. 3, pp. 1585--1602, 1999
 245 // (for pp ignore the exponential term)
 246 //
 247         TGraphErrors* grR3 = new TGraphErrors(*grB2);
 248         grR3->SetName(name.Data());
 249
 250         Double_t pt, b2;
 251
 252         for(Int_t i=0; i < grB2->GetN(); ++i)
 253         {
 254                 grB2->GetPoint(i, pt, b2);
 255
 256                 if(b2==0) continue;
 257                 Double_t r3 = Rside2Rlong(pt, b2, Cd);
 258
 259                 grR3->SetPoint(i, pt, r3);
 260
 261                 // only statistical error propagation of B2
 262                 Double_t errPt = grB2->GetErrorX(i);
 263                 Double_t errB2 = grB2->GetErrorY(i);
 264                 Double_t errR = r3*errB2/b2;
 265
 266                 grR3->SetPointError(i, errPt, errR);
 267         }
 268
 269         return grR3;
 270 }