B2 analysis code
[u/mrichter/AliRoot.git] / 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 #include <Riostream.h>
20 #include <TSystem.h>
21 #include <TROOT.h>
22 #include <TFileMerger.h>
23 #include <TString.h>
24 #include <TFile.h>
25 #include <TGraphErrors.h>
26
27 #include "AliLnB2.h"
28 #include "B2.h"
29 #include "Config.h"
30
31 Double_t GetCd(Double_t z)
32 {
33 //
34 // parameterization of <Cd> as a function of multiplicity
35 // from ALICE Rlong and Rside measurements
36 //
37         return 0.046133 + 0.0484458*z;
38 }
39
40 Int_t B2Mult(const TString& pSpectra     = "~/alice/output/Proton-lhc10d-Mult-Spectra.root",
41              const TString& ptag         = "lhc10d",
42              const TString& dSpectra     = "~/alice/output/Deuteron-lhc10d-Mult-Spectra.root",
43              const TString& dtag         = "lhc10d",
44              const TString& outputMultPt = "~/alice/output/B2-Mult-Pt.root",
45              const TString& outputPtMult = "~/alice/output/B2-Pt-Mult.root",
46              const TString& otag         = "lhc10d")
47 {
48 //
49 // B2 as a function of multiplicity
50 //
51         using namespace B2mult;
52         
53         const Int_t kNpart = 2;
54         
55         const Int_t kNMinPt = 0;
56         const Int_t kNMaxPt = 6;
57         
58         const TString kPrefix[] = { "", "Anti"};
59         const TString kSuffix[] = { "", "bar" };
60         Int_t kCharge[]         = {1, -1};
61         
62         // B2 as a function of pt for each multiplicity class
63         
64         for(Int_t i=0; i<kNmult; ++i)
65         {
66                 TFileMerger m;
67                 
68                 for(Int_t j=0; j<kNpart; ++j)
69                 {
70                         TString b2file = kPrefix[j] + "B2.root";
71                         
72                         AliLnB2 b2(pSpectra, ptag + "-" + kMultClass[i], dSpectra, dtag + "-" + kMultClass[i], b2file, otag + "-" + kMultClass[i], 2, kCharge[j]);
73                         
74                         b2.SetCd(GetCd(kKNOmult[i]));
75                         
76                         b2.Run();
77                         
78                         m.AddFile(b2file.Data(),0);
79                 }
80                 
81                 // merge B2 and B2bar
82                 
83                 TString outputfile = otag + "-" + kMultClass[i] + "-B2.root";
84                 
85                 m.OutputFile(outputfile.Data());
86                 m.Merge();
87                 
88                 gSystem->Exec("rm -f B2.root AntiB2.root");
89         }
90         
91         // merge multiplicity classes
92         
93         TFileMerger m;
94         
95         for(Int_t i=0; i<kNmult; ++i)
96         {
97                 TString b2 = otag + "-" + kMultClass[i] + "-B2.root";
98                 m.AddFile(b2.Data(),0);
99         }
100         
101         m.OutputFile(outputMultPt.Data());
102         m.Merge();
103         
104         // delete tmp files
105         
106         for(Int_t i=0; i<kNmult; ++i)
107         {
108                 gSystem->Exec(Form("rm -f %s-%s-B2.root",otag.Data(),kMultClass[i].Data()));
109         }
110         
111         // B2 as a function of multiplicity for each pt
112         
113         TFile* finput = new TFile(outputMultPt.Data());
114         if (finput->IsZombie()) exit(1);
115         
116         TGraphErrors* grB2pt[kNpart][kNmult];
117         TGraphErrors* grR3pt[kNpart][kNmult];
118         
119         for(Int_t i=0; i<kNpart; ++i)
120         {
121                 for(Int_t j=0; j<kNmult; ++j)
122                 {
123                         grB2pt[i][j] = (TGraphErrors*)FindObj(finput, otag + "-" + kMultClass[j], Form("B2%s_Pt", kSuffix[i].Data()));
124                         grR3pt[i][j] = (TGraphErrors*)FindObj(finput, otag + "-" + kMultClass[j], Form("R3%s_Pt", kSuffix[i].Data()));
125                 }
126         }
127         
128         TFile* foutput = new TFile(outputPtMult.Data(),"recreate");
129         
130         Double_t* pt = grB2pt[0][0]->GetX();
131         TString ptLabel[kNMaxPt];
132         
133         if(kNMaxPt > grB2pt[0][0]->GetN())
134         {
135                 std::cerr << "max pt too big" << std::endl;
136                 exit(1);
137         }
138         
139         for(Int_t i=kNMinPt; i<kNMaxPt; ++i)
140         {
141                 ptLabel[i] = Form("pt%.02fA",pt[i]);
142                 foutput->mkdir(ptLabel[i].Data());
143         }
144         
145         for(Int_t i=0; i<kNpart; ++i)
146         {
147                 for(Int_t j=kNMinPt; j<kNMaxPt; ++j)
148                 {
149                         Double_t B2[kNmult];
150                         Double_t B2StatErr[kNmult];
151                         
152                         Double_t R3[kNmult];
153                         Double_t R3StatErr[kNmult];
154                         
155                         for(Int_t k=0; k<kNmult; ++k)
156                         {
157                                 Double_t x, y, ey;
158                                 grB2pt[i][k]->GetPoint(j,x,y);
159                                 ey = grB2pt[i][k]->GetErrorY(j);
160                                 
161                                 B2[k] = y;
162                                 B2StatErr[k] = ey;
163                                 
164                                 grR3pt[i][k]->GetPoint(j,x,y);
165                                 ey = grR3pt[i][k]->GetErrorY(j);
166                                 
167                                 R3[k] = y;
168                                 R3StatErr[k] = ey;
169                         }
170                         
171                         TGraphErrors* grB2Mult = new TGraphErrors(kNmult, kKNOmult, B2, kKNOmultErr, B2StatErr);
172                         grB2Mult->SetName(Form("B2%s_Zmult", kSuffix[i].Data()));
173                         
174                         TGraphErrors* grR3Mult = new TGraphErrors(kNmult, kKNOmult, R3, kKNOmultErr, R3StatErr);
175                         grR3Mult->SetName(Form("R3%s_Zmult", kSuffix[i].Data()));
176                         
177                         foutput->cd(ptLabel[j].Data());
178                         
179                         grB2Mult->Write();
180                         grR3Mult->Write();
181                         
182                         delete grB2Mult;
183                         delete grR3Mult;
184                 }
185         }
186         
187         delete foutput;
188         delete finput;
189         
190         // particle ratios
191         
192         gROOT->ProcessLine(Form(".x RatioMult.C+g(\"%s\",\"%s\",\"%s\",\"%s\")",pSpectra.Data(), ptag.Data(),dSpectra.Data(), dtag.Data()));
193         
194         // draw B2 as a function of pt
195         
196         for(Int_t i=0; i<kNpart; ++i)
197         {
198                 gROOT->ProcessLine(Form(".x DrawDir.C+g(\"%s\",\"B2%s_Pt\",\"\",0,2, 1.e-3, 7.e-2,\"p_{T}/A (GeV/c)\",\"B_{2} (GeV^{2}/c^{3})\", 0,\"c%d.B2pt\",\"B2%spt\")", outputMultPt.Data(), kSuffix[i].Data(), i, kSuffix[i].Data()));
199                 
200                 gROOT->ProcessLine(Form(".x DrawDir.C+g(\"%s\",\"R3%s_Pt\",\"\",0,2, 0, 1.7,\"p_{T}/A (GeV/c)\",\"R_{side}^{2} R_{long} (fm^{3})\", 0,\"c%d.R3pt\",\"R3%spt\")", outputMultPt.Data(), kSuffix[i].Data(), i, kSuffix[i].Data()));
201         }
202         
203         // draw B2 as a function of z
204         
205         for(Int_t i=0; i<kNpart; ++i)
206         {
207                 gROOT->ProcessLine(Form(".x DrawDir.C+g(\"%s\",\"B2%s_Zmult\",\"\",0,7, 1.e-3, 7.e-2,\"z\",\"B_{2} (GeV^{2}/c^{3})\", 0,\"c%d.B2z\",\"B2%sZ\")", outputPtMult.Data(), kSuffix[i].Data(), i, kSuffix[i].Data()));
208                 
209                 gROOT->ProcessLine(Form(".x DrawDir.C+g(\"%s\",\"R3%s_Zmult\",\"\",0,8, 0, 4,\"z\",\"R_{side}^{2} R_{long} (fm^{3})\", 0,\"c%d.R3z\",\"R3%sZ\")", outputPtMult.Data(), kSuffix[i].Data(), i, kSuffix[i].Data()));
210         }
211         
212         return 0;
213 }