PWGLF/SPECTRA/Nuclei/B2/macros/B2Mult.C

   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 // B2 as a function of multiplicity
  17 // author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
  18
  19 #if !defined(__CINT__) || defined(__MAKECINT__)
  20 #include <Riostream.h>
  21 #include <TSystem.h>
  22 #include <TROOT.h>
  23 #include <TFileMerger.h>
  24 #include <TString.h>
  25 #include <TFile.h>
  26 #include <TGraphErrors.h>
  27 #include "AliLnBA.h"
  28 #endif
  29
  30 #include "B2.h"
  31 #include "Config.h"
  32
  33 Double_t GetCd(Double_t z)
  34 {
  35 //
  36 // parameterization of <Cd> as a function of multiplicity
  37 // from ALICE Rlong and Rside measurements
  38 //
  39         return 0.046133 + 0.0484458*z;
  40 }
  41
  42 extern void RatioMult(  const TString&, const TString&, const TString&, const TString&, const Int_t, const TString*, const Double_t*, const Double_t*, const TString&, const Bool_t, const TString& , const TString& );
  43
  44 Int_t B2Mult(  const TString& pSpectra     = "~/alice/output/Proton-lhc10d-Mult-Spectra.root"
  45              , const TString& ptag         = "lhc10d"
  46              , const TString& dSpectra     = "~/alice/output/Deuteron-lhc10d-Mult-Spectra.root"
  47              , const TString& dtag         = "lhc10d"
  48              , const TString& outputMultPt = "~/alice/output/B2-lhc10d-MultPt.root"
  49              , const TString& outputPtMult = "~/alice/output/B2-lhc10d-PtMult.root"
  50              , const TString& otag         = "lhc10d"
  51              , const Bool_t qTsallis       = 0
  52              , const TString& tsallisTag   = "Tsallis"
  53              , const TString& oratio       = "~/alice/output/Particle-Ratios-lhc10d-Tsallis.root")
  54 {
  55 //
  56 // B2 as a function of multiplicity
  57 //
  58         using namespace B2mult;
  59         using namespace std;
  60
  61         const Int_t kNpart = 2;
  62
  63         const Int_t kNMinPt = 0;
  64         const Int_t kNMaxPt = 6;
  65
  66         const TString kNucleus[kNpart] = { "Deuteron", "AntiDeuteron" };
  67
  68         // B2 as a function of pt for each multiplicity class
  69
  70         for(Int_t i=0; i<kNmult; ++i)
  71         {
  72                 TFileMerger m;
  73
  74                 const Int_t kZ[kNpart] = { 1, -1 };
  75                 const TString kB2File[kNpart] = {"B2.root", "AntiB2.root" };
  76
  77                 for(Int_t j=0; j<kNpart; ++j)
  78                 {
  79                         cout << kMultTag[i] << endl;
  80
  81                         AliLnBA b2(pSpectra, ptag + kMultTag[i], dSpectra, dtag + kMultTag[i], kB2File[j], otag + kMultTag[i], 2, kZ[j]);
  82
  83                         b2.SetCd(GetCd(kKNOmult[i]));
  84
  85                         b2.Run();
  86
  87                         m.AddFile(kB2File[j].Data(),0);
  88                 }
  89
  90                 // merge B2 and B2bar
  91
  92                 TString outputfile = otag + kMultTag[i] + "-B2.root";
  93
  94                 m.OutputFile(outputfile.Data());
  95                 m.Merge();
  96
  97                 gSystem->Exec(Form("rm -f %s %s",kB2File[0].Data(),kB2File[1].Data()));
  98         }
  99
 100         // merge multiplicity classes
 101
 102         TFileMerger m;
 103
 104         for(Int_t i=0; i<kNmult; ++i)
 105         {
 106                 TString b2 = otag + kMultTag[i] + "-B2.root";
 107                 m.AddFile(b2.Data(),0);
 108         }
 109
 110         m.OutputFile(outputMultPt.Data());
 111         m.Merge();
 112
 113         // delete tmp files
 114
 115         for(Int_t i=0; i<kNmult; ++i)
 116         {
 117                 gSystem->Exec(Form("rm -f %s%s-B2.root",otag.Data(),kMultTag[i].Data()));
 118         }
 119
 120         // B2 as a function of multiplicity for each pt
 121
 122         TFile* finput = new TFile(outputMultPt.Data());
 123         if (finput->IsZombie()) exit(1);
 124
 125         TGraphErrors* grB2pt[kNpart][kNmult];
 126         TGraphErrors* grR3pt[kNpart][kNmult];
 127         TGraphErrors* grSysB2pt[kNpart][kNmult];
 128         TGraphErrors* grSysR3pt[kNpart][kNmult];
 129
 130         for(Int_t i=0; i<kNpart; ++i)
 131         {
 132                 for(Int_t j=0; j<kNmult; ++j)
 133                 {
 134                         grB2pt[i][j] = FindObj<TGraphErrors>(finput, otag + kMultTag[j], kNucleus[i] + "_B2_Pt");
 135                         grR3pt[i][j] = FindObj<TGraphErrors>(finput, otag + kMultTag[j], kNucleus[i] + "_R3_Pt");
 136
 137                         grSysB2pt[i][j] = FindObj<TGraphErrors>(finput, otag + kMultTag[j], kNucleus[i] + "_SystErr_B2_Pt");
 138                         grSysR3pt[i][j] = FindObj<TGraphErrors>(finput, otag + kMultTag[j], kNucleus[i] + "_SystErr_R3_Pt");
 139                 }
 140         }
 141
 142         TFile* foutput = new TFile(outputPtMult.Data(),"recreate");
 143
 144         Double_t* pt = grB2pt[0][0]->GetX();
 145         TString ptLabel[kNMaxPt];
 146
 147         if(kNMaxPt > grB2pt[0][0]->GetN())
 148         {
 149                 std::cerr << "max pt too big" << std::endl;
 150                 exit(1);
 151         }
 152
 153         for(Int_t i=kNMinPt; i<kNMaxPt; ++i)
 154         {
 155                 ptLabel[i] = Form("pT %.02fA",pt[i]);
 156                 foutput->mkdir(ptLabel[i].Data());
 157         }
 158
 159         for(Int_t i=0; i<kNpart; ++i)
 160         {
 161                 for(Int_t j=kNMinPt; j<kNMaxPt; ++j)
 162                 {
 163                         Double_t B2[kNmult];
 164                         Double_t B2StatErr[kNmult];
 165                         Double_t B2SystErr[kNmult];
 166
 167                         Double_t R3[kNmult];
 168                         Double_t R3StatErr[kNmult];
 169                         Double_t R3SystErr[kNmult];
 170
 171                         for(Int_t k=0; k<kNmult; ++k)
 172                         {
 173                                 Double_t x, y, staterr, systerr;
 174                                 grB2pt[i][k]->GetPoint(j,x,y);
 175                                 staterr = grB2pt[i][k]->GetErrorY(j);
 176                                 systerr = grSysB2pt[i][k]->GetErrorY(j);
 177
 178                                 B2[k] = y;
 179                                 B2StatErr[k] = staterr;
 180                                 B2SystErr[k] = systerr;
 181
 182                                 grR3pt[i][k]->GetPoint(j,x,y);
 183                                 staterr = grR3pt[i][k]->GetErrorY(j);
 184                                 systerr = grSysR3pt[i][k]->GetErrorY(j);
 185
 186                                 R3[k] = y;
 187                                 R3StatErr[k] = staterr;
 188                                 R3SystErr[k] = systerr;
 189                         }
 190
 191                         TGraphErrors* grB2Mult = new TGraphErrors(kNmult, kKNOmult, B2, kKNOmultErr, B2StatErr);
 192                         grB2Mult->SetName(Form("%s_B2_Zmult", kNucleus[i].Data()));
 193
 194                         TGraphErrors* grR3Mult = new TGraphErrors(kNmult, kKNOmult, R3, kKNOmultErr, R3StatErr);
 195                         grR3Mult->SetName(Form("%s_R3_Zmult", kNucleus[i].Data()));
 196
 197                         Double_t zMultSystErr[kNmult];
 198                         for(Int_t k=0; k<kNmult; ++k) zMultSystErr[k]=0.07;
 199
 200                         TGraphErrors* grSysB2Mult = new TGraphErrors(kNmult, kKNOmult, B2, zMultSystErr, B2SystErr);
 201                         grSysB2Mult->SetName(Form("%s_SystErr_B2_Zmult", kNucleus[i].Data()));
 202
 203                         TGraphErrors* grSysR3Mult = new TGraphErrors(kNmult, kKNOmult, R3, zMultSystErr, R3SystErr);
 204                         grSysR3Mult->SetName(Form("%s_SystErr_R3_Zmult", kNucleus[i].Data()));
 205
 206                         foutput->cd(ptLabel[j].Data());
 207
 208                         grB2Mult->Write();
 209                         grR3Mult->Write();
 210                         grSysB2Mult->Write();
 211                         grSysR3Mult->Write();
 212
 213                         delete grB2Mult;
 214                         delete grR3Mult;
 215                         delete grSysB2Mult;
 216                         delete grSysR3Mult;
 217                 }
 218         }
 219
 220         delete foutput;
 221         delete finput;
 222
 223         //
 224         // Particle ratios
 225         // -----------------------------------
 226
 227         RatioMult(pSpectra, dSpectra, ptag, dtag, kNmult, kMultTag, kKNOmult, kKNOmultErr, kKNOmultName, qTsallis, oratio, tsallisTag);
 228
 229
 230         // draw B2 as a function of pt
 231
 232         for(Int_t i=0; i<kNpart; ++i)
 233         {
 234                 gROOT->ProcessLine(Form(".x DrawDir.C+g(\"%s\",\"%s_B2_Pt\",\"\",0,2, 1.e-3, 7.e-2,\"p_{T}/A (GeV/c)\",\"B_{2} (GeV^{2}/c^{3})\", 0,\"c%d.B2pt\",\"%s_B2_Pt\")", outputMultPt.Data(), kNucleus[i].Data(), i, kNucleus[i].Data()));
 235
 236                 gROOT->ProcessLine(Form(".x DrawDir.C+g(\"%s\",\"%s_R3_Pt\",\"\",0,2, 0, 1.7,\"p_{T}/A (GeV/c)\",\"R_{side}^{2} R_{long} (fm^{3})\", 0,\"c%d.R3pt\",\"%s_R3_Pt\")", outputMultPt.Data(), kNucleus[i].Data(), i, kNucleus[i].Data()));
 237         }
 238
 239         // draw B2 as a function of z
 240
 241         for(Int_t i=0; i<kNpart; ++i)
 242         {
 243                 gROOT->ProcessLine(Form(".x DrawDir.C+g(\"%s\",\"%s_B2_Zmult\",\"\",0,5, 3.e-3, 6.e-2,\"z\",\"B_{2} (GeV^{2}/c^{3})\", 0,\"c%d.B2z\",\"%s_B2_Zmult\")", outputPtMult.Data(), kNucleus[i].Data(), i, kNucleus[i].Data()));
 244
 245                 gROOT->ProcessLine(Form(".x DrawDir.C+g(\"%s\",\"%s_R3_Zmult\",\"\",0,5, 0, 4,\"z\",\"R_{side}^{2} R_{long} (fm^{3})\", 0,\"c%d.R3z\",\"%s_R3_Zmult\")", outputPtMult.Data(), kNucleus[i].Data(), i, kNucleus[i].Data()));
 246         }
 247
 248         return 0;
 249 }