]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/FORWARD/analysis/AliFMDDndeta.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
add maximum M02 band cut, retune fit param, define temporary m02 cut for eta and...
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis / AliFMDDndeta.cxx
CommitLineData
7c1a1f1d 1// Code to analyse dN/deta from the forward analysis
2// This can plot the results
3// Also works for MC data
4//
5// -- Author: Hans Hjersing Dalsgaard <canute@gmail.com>
c5058a61 6#include "AliFMDDndeta.h"
7#include "TFile.h"
8#include "AliLog.h"
9#include "TH1.h"
10#include "AliFMDAnaParameters.h"
7e2bf482 11#include "AliFMDAnaCalibSharingEfficiency.h"
c5058a61 12#include "TStyle.h"
13//#include "TObjArray.h"
14#include "TCanvas.h"
15#include "TLine.h"
2dcdd4eb 16#include "TGraphErrors.h"
541c19ed 17#include "TGraphAsymmErrors.h"
c5058a61 18#include "TPad.h"
19#include "iostream"
059c7c6b 20#include "TH3.h"
c5058a61 21#include "TMath.h"
059c7c6b 22#include "TProfile.h"
23#include "TProfile2D.h"
24#include "TProfile3D.h"
c5058a61 25#include "TLegend.h"
059c7c6b 26#include "TPad.h"
c5058a61 27#include "TLatex.h"
059c7c6b 28#include "TStyle.h"
29#include "TF1.h"
2dcdd4eb 30
31#define SMALLNUMBER 0.0001
32
c5058a61 33ClassImp(AliFMDDndeta)
34//_____________________________________________________________________
35
36AliFMDDndeta::AliFMDDndeta()
37: TObject(),
38 fList(0),
39 fMultList(),
40 fNbinsToCut(2),
2dcdd4eb 41 fVtxCut1(-10),
42 fVtxCut2(10),
c5058a61 43 fIsInit(kFALSE),
44 fIsGenerated(),
ae26bdd7 45 fPrimEvents(),
46 fEvents(),
059c7c6b 47 fPrimdNdeta(),
48 fDrawAll(kFALSE)
c5058a61 49{
541c19ed 50 //AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
059c7c6b 51 /* fDataObject = new TProfile3D("dataObject","dataObject",
52 pars->GetNetaBins(),-6,6,
53 20,0,2*TMath::Pi(),
54 pars->GetNvtxBins(),-0.5,pars->GetNvtxBins()-0.5);
55 fDataObject->SetXTitle("#eta");
56 fDataObject->SetYTitle("#varphi [radians]");
57 fDataObject->SetZTitle("v_{z} [cm]");*/
58
70d74659 59 fAnalysisNames[0] = "Hits";
60 fAnalysisNames[1] = "HitsTrVtx";
61 fAnalysisNames[2] = "dNdeta";
62 fAnalysisNames[3] = "dNdetaTrVtx";
059c7c6b 63 fAnalysisNames[4] = "dNdetaNSD";
2dcdd4eb 64
059c7c6b 65 for(Int_t i=0; i<5;i++)
70d74659 66 fMultList[i] = new TList();
c5058a61 67}
68//_____________________________________________________________________
69void AliFMDDndeta::SetNames(Analysis what) {
2dcdd4eb 70 // Set names of histograms from analysis
c5058a61 71
72 switch(what) {
73 case kHits :
74 fPrimEvents.Form("nMCEventsNoCuts"); //was nMCEvents
75 fEvents.Form("nEvents");
76 fPrimdNdeta.Form("hMultvsEtaNoCuts");
77 break;
7e2bf482 78 case kHitsTrVtx :
79 fPrimEvents.Form("nMCEvents");
80 fEvents.Form("nEvents");
81 fPrimdNdeta.Form("hMultvsEtaNoCuts");
82 break;
c5058a61 83 case kMult :
84 fPrimEvents.Form("nMCEventsNoCuts");
85 fEvents.Form("nEvents");
86 fPrimdNdeta.Form("hMultvsEtaNoCuts");
87 break;
88 case kMultTrVtx :
89 fPrimEvents.Form("nMCEvents");
90 fEvents.Form("nEvents");
91 fPrimdNdeta.Form("hMultvsEta");
92 break;
059c7c6b 93 case kMultNSD :
94 fPrimEvents.Form("nMCEventsNSDNoCuts");
95 fEvents.Form("nNSDEvents");
96 fPrimdNdeta.Form("hMultvsEtaNSDNoCuts");
97 break;
98
c5058a61 99 default:
100 break;
101 }
102}
103//_____________________________________________________________________
104const char* AliFMDDndeta::GetAnalysisName(Analysis what, UShort_t det, Char_t ring, Int_t vtxbin) {
2dcdd4eb 105 //Get names of histograms
ae26bdd7 106 const char* name = "";
c5058a61 107
108 switch(what) {
109 case kHits :
110 name = Form("hits_NoCuts_FMD%d%c_vtxbin%d_proj",det,ring,vtxbin); //NoCuts added
111 break;
7e2bf482 112 case kHitsTrVtx :
113 name = Form("hits_FMD%d%c_vtxbin%d_proj",det,ring,vtxbin);
114 break;
c5058a61 115 case kMult :
116 name = Form("dNdeta_FMD%d%c_vtxbin%d_proj",det,ring,vtxbin);
117 break;
118 case kMultTrVtx :
119 name = Form("dNdeta_FMD%d%c_TrVtx_vtxbin%d_proj",det,ring,vtxbin);
120 break;
059c7c6b 121 case kMultNSD :
122 name = Form("dNdetaNSD_FMD%d%c_vtxbin%d_proj",det,ring,vtxbin);
123 break;
124
c5058a61 125 default:
126 break;
127 }
128 //std::cout<<name.Data()<<std::endl;
129 return name;
130}
131//_____________________________________________________________________
132const char* AliFMDDndeta::GetPrimName(Analysis what, UShort_t det, Char_t ring, Int_t vtxbin) {
2dcdd4eb 133 //Get names of primaries
ae26bdd7 134 const char* name = "";
c5058a61 135
136 switch(what) {
137 case kHits :
138 name = Form("hMCHits_nocuts_FMD%d%c_vtxbin%d",det,ring,vtxbin); //nocuts added
139 break;
7e2bf482 140 case kHitsTrVtx :
141 name = Form("hMCHits_FMD%d%c_vtxbin%d",det,ring,vtxbin);
142 break;
c5058a61 143 case kMult :
144 name = Form("primmult_vtxbin%d",vtxbin);
145 break;
146 case kMultTrVtx :
147 name = Form("primmult_NoCuts_vtxbin%d",vtxbin);
148 break;
059c7c6b 149 case kMultNSD :
150 name = Form("primmult_NSD_vtxbin%d",vtxbin);
151 break;
c5058a61 152 default:
153 break;
154 }
155 return name;
156}
157//_____________________________________________________________________
2dcdd4eb 158void AliFMDDndeta::GenerateHits(Analysis what) {
159 //Generate the hits distributions
c5058a61 160 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
161 TH2F* hTmp = pars->GetBackgroundCorrection(1,'I',0);
162 Int_t nVertexBins = pars->GetNvtxBins();
163 Int_t nEtaBins = hTmp->GetNbinsX();
164 Float_t etaMin = hTmp->GetXaxis()->GetXmin();
165 Float_t etaMax = hTmp->GetXaxis()->GetXmax();
166
059c7c6b 167 for(Int_t i = 0; i<nVertexBins;i++) { //
df2a9c32 168 TH1F* hHits = new TH1F(Form("hMCHits_vtxbin%d_%s",i,fAnalysisNames[what].Data()),Form("hHits_vtxbin%d_%s",i,fAnalysisNames[what].Data()),nEtaBins,etaMin,etaMax);
c5058a61 169 hHits->Sumw2();
70d74659 170 fMultList[what]->Add(hHits);
c5058a61 171 }
2dcdd4eb 172
c5058a61 173 for(Int_t det = 1; det<=3; det++) {
174 Int_t maxRing = (det == 1 ? 0 : 1);
175 for(Int_t ring = 0;ring<=maxRing;ring++) {
176
177 Char_t ringChar = (ring == 0 ? 'I' : 'O');
178 for(Int_t v=0; v< nVertexBins; v++) {
2dcdd4eb 179 TH1F* hits = (TH1F*)fList->FindObject(GetPrimName(what,det,ringChar,v));
c5058a61 180
181 Int_t nNonZero = 0, nNonZeroInData = 0;
182
183 //removing edges
184 for(Int_t i =1;i<=hits->GetNbinsX();i++) {
185 if(hits->GetBinContent(i) !=0)
186 nNonZero++;
187 }
188
df2a9c32 189 TH1F* sumMultHist = (TH1F*)fMultList[what]->FindObject(Form("hMCHits_vtxbin%d_%s",v,fAnalysisNames[what].Data()));
c5058a61 190 for(Int_t i =1;i<=sumMultHist->GetNbinsX();i++) {
191 if(hits->GetBinContent(i) == 0 ) continue;
192
193 nNonZeroInData++;
194
195 if(nNonZeroInData <=fNbinsToCut || nNonZeroInData > (nNonZero - fNbinsToCut)) {
196 continue;
197 }
2dcdd4eb 198 if(sumMultHist->GetBinContent(i) < SMALLNUMBER && TMath::Abs(hits->GetBinContent(i)) > 0){
c5058a61 199 sumMultHist->SetBinContent(i,hits->GetBinContent(i));
200 sumMultHist->SetBinError(i,hits->GetBinError(i));
201 }
202 else {
203
204
205 Float_t sumofw = (1/TMath::Power(hits->GetBinError(i),2)) + (1/TMath::Power(sumMultHist->GetBinError(i),2));
206 Float_t wav = (hits->GetBinContent(i)*(1/TMath::Power(hits->GetBinError(i),2)) + sumMultHist->GetBinContent(i)*(1/TMath::Power(sumMultHist->GetBinError(i),2)))/sumofw;
207 Float_t error = 1/TMath::Sqrt(sumofw);
208 sumMultHist->SetBinContent(i, wav);
209 sumMultHist->SetBinError(i, error);
210 }
211 }
212 }
213 }
214 }
215}
216
217//_____________________________________________________________________
218void AliFMDDndeta::Init(const Char_t* filename) {
2dcdd4eb 219 //Initialize everything from file
c5058a61 220 TFile* fin = TFile::Open(filename);
221
222 if(!fin) {
2dcdd4eb 223 AliWarning("No file - exiting !");
c5058a61 224 return;
225 }
c5058a61 226 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
707ec2bd 227 //pars->Init();
c5058a61 228
507687cd 229 TList* list = (TList*)fin->Get(Form("%s/BackgroundCorrected",pars->GetDndetaAnalysisName()));
2dcdd4eb 230
231 if(!list) //an old file ? Perhaps...
232 list = (TList*)fin->Get("BackgroundCorrected");
70d74659 233
2dcdd4eb 234 Init(list);
707ec2bd 235
2dcdd4eb 236}
237//_____________________________________________________________________
238void AliFMDDndeta::Init(TList* list) {
239 //Initialize everything from TList
240
241 if(!list) {
242 AliWarning("No list - exiting !");
243 return;
244 }
70d74659 245
246 fList = (TList*)list->Clone("inputList");
707ec2bd 247
c5058a61 248 fIsGenerated[kHits] = kFALSE;
249 fIsGenerated[kMult] = kFALSE;
250 fIsGenerated[kMultTrVtx] = kFALSE;
2dcdd4eb 251 fIsGenerated[kHitsTrVtx] = kFALSE;
059c7c6b 252 fIsGenerated[kMultNSD] = kFALSE;
c5058a61 253
254 fIsInit = kTRUE;
255}
256//_____________________________________________________________________
257void AliFMDDndeta::GenerateMult(Analysis what) {
2dcdd4eb 258 //Generate dNdeta
c5058a61 259
260 if(!fIsInit) {
261 AliWarning("Not initialised - call Init to remedy");
262 return;
263 }
264
265 if(fIsGenerated[what]) {
266 AliWarning("Already generated - have a look at the results!");
267 return;
268 }
269 else fIsGenerated[what] = kTRUE;
270
271 SetNames(what);
272
2dcdd4eb 273 if(what == kHits || what == kHitsTrVtx)
274 GenerateHits(what);
c5058a61 275
276 TH1I* hMCEvents = (TH1I*)fList->FindObject(fPrimEvents.Data());
277
278 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
279
280 TH2F* hTmp = pars->GetBackgroundCorrection(1,'I',0);
281 Int_t nVertexBins = pars->GetNvtxBins();
282 Int_t nEtaBins = hTmp->GetNbinsX();
283 Float_t etaMin = hTmp->GetXaxis()->GetXmin();
284 Float_t etaMax = hTmp->GetXaxis()->GetXmax();
285
286 for(Int_t i = 0; i<nVertexBins;i++) {
df2a9c32 287 TH1F* hMult = new TH1F(Form("hMult_vtxbin%d_%s",i,fAnalysisNames[what].Data()),Form("hMult_vtxbin%d_%s",i,fAnalysisNames[what].Data()),nEtaBins,etaMin,etaMax);
c5058a61 288 hMult->Sumw2();
70d74659 289 fMultList[what]->Add(hMult);
c5058a61 290 }
291
292 for(Int_t det = 1; det<=3;det++) {
293 Int_t maxRing = (det == 1 ? 0 : 1);
294 for(Int_t iring = 0; iring<=maxRing; iring++) {
295 Char_t ringChar = (iring == 0 ? 'I' : 'O');
059c7c6b 296 //Int_t nsec = (iring == 0 ? 20 : 40);
df2a9c32 297 TH1F* hRingMult= new TH1F(Form("hRingMult_FMD%d%c_%s",det,ringChar,fAnalysisNames[what].Data()),Form("hRingMult_FMD%d%c_%s",det,ringChar,fAnalysisNames[what].Data()),nEtaBins,etaMin,etaMax);
70d74659 298 fMultList[what]->Add(hRingMult);
059c7c6b 299 // TProfile* phiprofile = new TProfile(Form("dNdphiFMD%d%c",det,ringChar), Form("dNdphiFMD%d%c;#Phi",det,ringChar), nsec , 0, 2*TMath::Pi());
300 //fMultList[what]->Add(phiprofile);
c5058a61 301 }
302 }
059c7c6b 303 TH1I* hEvents = (TH1I*)fList->FindObject(fEvents.Data());
c5058a61 304 TH1F* hPrimVtx = 0;
305
306
307
308 for(Int_t det = 1; det<=3; det++) {
309 Int_t maxRing = (det == 1 ? 0 : 1);
310 for(Int_t ring = 0;ring<=maxRing;ring++) {
311
312 Char_t ringChar = (ring == 0 ? 'I' : 'O');
313 for(Int_t v=0; v< nVertexBins; v++) {
314 if(det == 1) {
2dcdd4eb 315 if(what == kHits || what == kHitsTrVtx)
df2a9c32 316 hPrimVtx = (TH1F*)fMultList[what]->FindObject(Form("hMCHits_vtxbin%d_%s",v,fAnalysisNames[what].Data()));
c5058a61 317 else
318 hPrimVtx = (TH1F*)fList->FindObject(GetPrimName(what,det,ringChar,v));
319
320
321 Float_t nPrimevents = hMCEvents->GetBinContent(v+1);
322 Float_t xb = hPrimVtx->GetNbinsX();
323 Float_t xr = hPrimVtx->GetXaxis()->GetXmax() - hPrimVtx->GetXaxis()->GetXmin();
324 hPrimVtx->Scale(xb / xr );
325 if(nPrimevents > 0)
326 hPrimVtx->Scale(1/nPrimevents);
327
328 }
329
330 Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
331 Float_t vtxZ1 = (delta*v) - pars->GetVtxCutZ();
332 Float_t vtxZ2 = (delta*(v+1)) - pars->GetVtxCutZ();
333
2dcdd4eb 334 if(vtxZ1< fVtxCut1 || vtxZ2 >fVtxCut2)
c5058a61 335 continue;
336 Float_t nEvents = (Float_t)hEvents->GetBinContent(v+1);
337
338 //TH1F* multhistproj = (TH1F*)fList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d_proj",det,ringChar,v));
059c7c6b 339
340
341
342
343
344
345 TH1F* multhistproj = (TH1F*)fList->FindObject(GetAnalysisName(what,det,ringChar,v));
c5058a61 346 if(nEvents)
347 multhistproj->Scale(1/nEvents);
348 Float_t xnBins = multhistproj->GetNbinsX();
349 Float_t xrange = multhistproj->GetXaxis()->GetXmax() - multhistproj->GetXaxis()->GetXmin();
350 multhistproj->Scale(xnBins / xrange);
351
541c19ed 352 //TH2F* multhist = (TH2F*)fList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,v));
c5058a61 353
354
059c7c6b 355 //if(nEvents)
356 // multhist->Scale(1/nEvents);
357 /* for(Int_t i=1;i<=multhist->GetNbinsX();i++) {
358 for(Int_t j=1;j<=multhist->GetNbinsY();j++) {
359
360 if (multhist->GetBinContent(i,j) <= 0.0001) continue;
361 //std::cout<<multhist->GetBinContent(i,j)<<std::endl;
362 fDataObject->Fill(multhist->GetXaxis()->GetBinCenter(i),
363 multhist->GetYaxis()->GetBinCenter(j),
364 v,
365 multhist->GetBinContent(i,j), multhist->GetBinError(i,j));
366 //fDataObject->SetBinContent(i,j,v,multhist->GetBinContent(i,j));
367 //fDataObject->SetBinError(i,j,v,multhist->GetBinError(i,j));
368 }
369 }
370 */
371 //if(nEvents)
372 // multhist->Scale(1/nEvents);
373 //if(ringChar == 'O')
374 // multhist->RebinY(2);
c5058a61 375 //removing edges
376
059c7c6b 377 Int_t nNonZero = 0, nNonZeroInData = 0;
378
c5058a61 379 for(Int_t i =1;i<=multhistproj->GetNbinsX();i++) {
380 if(multhistproj->GetBinContent(i) !=0)
381 nNonZero++;
382 }
7e2bf482 383 Int_t nBinsOld = fNbinsToCut;
4b8bdb60 384 //if(det == 1 && ringChar =='I') {
385 // fNbinsToCut = 0;
386 // }
df2a9c32 387 TH1F* hRingMult = (TH1F*)fMultList[what]->FindObject(Form("hRingMult_FMD%d%c_%s",det,ringChar,fAnalysisNames[what].Data()));
c5058a61 388
389 for(Int_t i=1; i<=hRingMult->GetNbinsX(); i++) {
390 if(multhistproj->GetBinContent(i)!=0) {
391 nNonZeroInData++;
392 if(nNonZeroInData <=fNbinsToCut || nNonZeroInData > (nNonZero - fNbinsToCut)) {
393 continue;
394 }
395 Float_t oldweight = 0;
396 Float_t oldwav = 0;
397 if(hRingMult->GetBinError(i)>0) {
398 oldweight = 1/TMath::Power(hRingMult->GetBinError(i),2);
399 oldwav = oldweight*hRingMult->GetBinContent(i);
7e2bf482 400 }
c5058a61 401 Float_t weight = 1/TMath::Power(multhistproj->GetBinError(i),2);
402 Float_t wav = oldwav + weight*multhistproj->GetBinContent(i);
403 Float_t sumofw = oldweight + weight;
404 if(sumofw) {
405 Float_t error = 1/TMath::Sqrt(sumofw);
406
407 hRingMult->SetBinContent(i,wav/sumofw);
408 hRingMult->SetBinError(i,error);
409
410 }
411 }
412 }
413 nNonZeroInData = 0;
414
df2a9c32 415 TH1F* sumMultHist = (TH1F*)fMultList[what]->FindObject(Form("hMult_vtxbin%d_%s",v,fAnalysisNames[what].Data()));
ae26bdd7 416
c5058a61 417 for(Int_t i =1;i<=sumMultHist->GetNbinsX();i++) {
7e2bf482 418 if(multhistproj->GetBinContent(i) != 0 ) {
419 nNonZeroInData++;
c5058a61 420
7e2bf482 421 if(nNonZeroInData <=fNbinsToCut || nNonZeroInData > (nNonZero - fNbinsToCut)) {
422 continue;
423 }
2dcdd4eb 424 if(sumMultHist->GetBinContent(i) < SMALLNUMBER && TMath::Abs(multhistproj->GetBinContent(i)) > 0){
c5058a61 425 sumMultHist->SetBinContent(i,multhistproj->GetBinContent(i));
426 sumMultHist->SetBinError(i,multhistproj->GetBinError(i));
427 }
428 else {
429
430
431 Float_t sumofw = (1/TMath::Power(multhistproj->GetBinError(i),2)) + (1/TMath::Power(sumMultHist->GetBinError(i),2));
432 Float_t wav = (multhistproj->GetBinContent(i)*(1/TMath::Power(multhistproj->GetBinError(i),2)) + sumMultHist->GetBinContent(i)*(1/TMath::Power(sumMultHist->GetBinError(i),2)))/sumofw;
433 Float_t error = 1/TMath::Sqrt(sumofw);
434 sumMultHist->SetBinContent(i, wav);
435 sumMultHist->SetBinError(i, error);
436 }
437
438
7e2bf482 439 }
c5058a61 440 }
7e2bf482 441 fNbinsToCut = nBinsOld;
c5058a61 442 }
7e2bf482 443
c5058a61 444 }
445 }
446
df2a9c32 447 TH1F* sumMult = new TH1F(Form("dNdeta_%s",fAnalysisNames[what].Data()),Form("dNdeta_%s",fAnalysisNames[what].Data()),nEtaBins,etaMin,etaMax);
c5058a61 448 sumMult->Sumw2();
70d74659 449 fMultList[what]->Add(sumMult);
df2a9c32 450 TH1F* primMult = new TH1F(Form("primary_%s",fAnalysisNames[what].Data()),Form("primary_%s",fAnalysisNames[what].Data()),nEtaBins,etaMin,etaMax);
c5058a61 451 primMult->Sumw2();
70d74659 452 fMultList[what]->Add(primMult);
c5058a61 453
454 Float_t wav = 0, sumofw=0, weight = 0, primwav = 0, primsumofw=0, primweight = 0;//, etaofbin = 0;
455 for(Int_t i =1; i<=sumMult->GetNbinsX();i++) {
456
457 for(Int_t v = 0; v<nVertexBins;v++) {
2dcdd4eb 458 if(what == kHits || what == kHitsTrVtx)
df2a9c32 459 hPrimVtx = (TH1F*)fMultList[what]->FindObject(Form("hMCHits_vtxbin%d_%s",v,fAnalysisNames[what].Data()));
c5058a61 460 else
461 hPrimVtx = (TH1F*)fList->FindObject(GetPrimName(what,0,0,v));
df2a9c32 462 TH1F* sumMultHist = (TH1F*)fMultList[what]->FindObject(Form("hMult_vtxbin%d_%s",v,fAnalysisNames[what].Data()));
c5058a61 463 //etaofbin += sumMultHist->GetBinCenter(i);
464 // if( hPrimVtx->GetBinContent(i)!=0 && hPrimVtx->GetBinError(i)>0.001 && sumMultHist->GetBinContent(i)!=0) {
2dcdd4eb 465 if( TMath::Abs(hPrimVtx->GetBinContent(i)) > 0) {
c5058a61 466 primweight = 1/TMath::Power(hPrimVtx->GetBinError(i),2);
467 primwav = primwav + primweight*hPrimVtx->GetBinContent(i);
468 primsumofw = primsumofw + primweight;
469 }
2dcdd4eb 470 if(TMath::Abs(sumMultHist->GetBinContent(i)) > 0) {
c5058a61 471
472 weight = 1/TMath::Power(sumMultHist->GetBinError(i),2);
473 wav = wav + weight*sumMultHist->GetBinContent(i);
474 sumofw = sumofw + weight;
475 }
476
477 }
478 if( primsumofw !=0 ) {// && sumofw !=0) {
479 Float_t primerror = 1/TMath::Sqrt(primsumofw);
480
481 primMult->SetBinContent(i,primwav/primsumofw);
482 primMult->SetBinError(i,primerror);
483 }
484 else {
485 primMult->SetBinContent(i,0);
486 primMult->SetBinError(i,0);
487 }
488
489 if(sumofw) {
490
491 Float_t error = 1/TMath::Sqrt(sumofw);
492
493 sumMult->SetBinContent(i,wav/sumofw);
494 sumMult->SetBinError(i,error);
495 }
496
2dcdd4eb 497
c5058a61 498 wav = 0;
499 sumofw = 0;
500 weight = 0;
501 primwav = 0;
502 primsumofw = 0;
503 primweight = 0;
504
505 }
506
507}
508//_____________________________________________________________________
541c19ed 509void AliFMDDndeta::DrawDndeta(Analysis what, Int_t rebin, Bool_t realdata, TString pythiafile) {
2dcdd4eb 510 //Draw everything
c5058a61 511 if(!fIsInit) {
512 AliWarning("Not initialised - call Init to remedy");
513 return;
514 }
515
516 if(!fIsGenerated[what]) {
517 AliWarning("Not generated - generate first then draw!");
518 return;
519 }
541c19ed 520 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
c5058a61 521 SetNames(what);
059c7c6b 522 //UA5 data NSD
523 Float_t x[19] = {0.125,0.375,0.625,0.875,1.125,1.375,1.625,1.875,2.125,2.375,
524 2.625,2.875,3.125,3.375,3.625,3.875,4.125,4.375,4.625};
525 Float_t x2[19] = {-0.125,-0.375,-0.625,-0.875,-1.125,-1.375,-1.625,-1.875,
526 -2.125,-2.375,-2.625,-2.875,-3.125,-3.375,-3.625,-3.875,
527 -4.125,-4.375,-4.625};
528
529
530 Float_t y[19] = {3.48, 3.38, 3.52, 3.68, 3.71, 3.86, 3.76, 3.66, 3.72 ,3.69,
531 3.56, 3.41, 3.15, 3.09, 2.74, 2.73, 2.32, 1.99, 1.69};
532
533 Float_t ey[19] = {0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,
534 0.07,0.07,0.08,0.08,0.09,0.09,0.1,0.13};
535
536 TGraphErrors* graph = new TGraphErrors(19,x,y,NULL,ey);
537 TGraphErrors* graph2 = new TGraphErrors(19,x2,y,NULL,ey);
538 //UA5 data INEL
539 Float_t xinel[19] = {0.125, 0.375, 0.625, 0.875, 1.125, 1.375, 1.625, 1.875, 2.125, 2.375, 2.625, 2.875, 3.125, 3.375, 3.625, 3.875, 4.125, 4.375, 4.625};
540 Float_t xinel2[19] = {-0.125, -0.375, -0.625, -0.875, -1.125, -1.375, -1.625, -1.875, -2.125, -2.375, -2.625, -2.875, -3.125, -3.375, -3.625, -3.875, -4.125, -4.375, -4.625};
541 Float_t yinel[19] = {3.14, 3.04, 3.17, 3.33, 3.33, 3.53, 3.46, 3.41, 3.45, 3.39, 3.07, 3.07, 2.93, 2.93, 2.55, 2.48, 2.18, 1.91, 1.52};
542 Float_t eyinel[19] = {0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.08, 0.08, 0.09, 0.09, 0.1, 0.13};
543 TGraphErrors* graphinel = new TGraphErrors(19,xinel,yinel,NULL,eyinel);
544 TGraphErrors* graphinel2 = new TGraphErrors(19,xinel2,yinel,NULL,eyinel);
545
546 graph->SetMarkerStyle(22);
547 graph->SetMarkerColor(kBlue);
548 graphinel->SetMarkerStyle(20);
549 graphinel->SetMarkerColor(kGreen);
550 graphinel2->SetMarkerStyle(24);
551 graphinel2->SetMarkerColor(kGreen);
552 //graph->Draw("P");
553 graph2->SetMarkerStyle(26);
554 graph2->SetMarkerColor(kBlue);
555 graphinel->GetHistogram()->SetStats(kFALSE);
556 graphinel2->GetHistogram()->SetStats(kFALSE);
557
558
559 //End UA5 data
c5058a61 560
541c19ed 561 //Published ALICE data
562 // Plot: p7742_d1x1y1
507687cd 563
564 // INEL
565
7c1a1f1d 566 TGraphAsymmErrors* p7742D1x1y1 = 0;
567 double p7742D1x1y1xval[] = { -1.3, -1.1, -0.9, -0.7, -0.5, -0.3, -0.1, 0.1, 0.3,
541c19ed 568 0.5, 0.7, 0.9, 1.1, 1.3 };
7c1a1f1d 569 double p7742D1x1y1xerrminus[] = { 0.09999999999999987, 0.09999999999999987, 0.09999999999999998, 0.10000000000000009, 0.09999999999999998, 0.10000000000000003, 0.1, 0.1, 0.09999999999999998,
541c19ed 570 0.09999999999999998, 0.09999999999999998, 0.09999999999999998, 0.10000000000000009, 0.10000000000000009 };
7c1a1f1d 571 double p7742D1x1y1xerrplus[] = { 0.10000000000000009, 0.10000000000000009, 0.09999999999999998, 0.09999999999999998, 0.09999999999999998, 0.09999999999999998, 0.1, 0.1, 0.10000000000000003,
541c19ed 572 0.09999999999999998, 0.10000000000000009, 0.09999999999999998, 0.09999999999999987, 0.09999999999999987 };
7c1a1f1d 573 double p7742D1x1y1yval[] = { 3.28, 3.28, 3.22, 3.12, 3.06, 3.02, 2.98, 3.02, 3.02,
541c19ed 574 3.05, 3.15, 3.21, 3.26, 3.33 };
7c1a1f1d 575 double p7742D1x1y1yerrminus[] = { 0.06324555320336758, 0.06324555320336758, 0.06324555320336758, 0.06324555320336758, 0.06324555320336758, 0.05385164807134505, 0.05385164807134505, 0.05385164807134505, 0.05385164807134505,
541c19ed 576 0.06324555320336758, 0.06324555320336758, 0.06324555320336758, 0.06324555320336758, 0.06324555320336758 };
7c1a1f1d 577 double p7742D1x1y1yerrplus[] = { 0.08246211251235322, 0.08246211251235322, 0.08246211251235322, 0.08246211251235322, 0.08246211251235322, 0.08246211251235322, 0.07280109889280519, 0.08246211251235322, 0.08246211251235322,
541c19ed 578 0.08246211251235322, 0.08246211251235322, 0.08246211251235322, 0.08246211251235322, 0.08246211251235322 };
7c1a1f1d 579 int p7742D1x1y1numpoints = 14;
580 p7742D1x1y1 = new TGraphAsymmErrors(p7742D1x1y1numpoints, p7742D1x1y1xval, p7742D1x1y1yval, p7742D1x1y1xerrminus, p7742D1x1y1xerrplus, p7742D1x1y1yerrminus, p7742D1x1y1yerrplus);
581 p7742D1x1y1->SetName("/HepData/7742/d1x1y1");
582 p7742D1x1y1->SetTitle("/HepData/7742/d1x1y1");
583 p7742D1x1y1->SetMarkerStyle(21);
584 p7742D1x1y1->SetMarkerColor(kRed);
585 // p7742D1x1y1->Draw("Psame");
507687cd 586
587
588 //NSD
7c1a1f1d 589 TGraphAsymmErrors* p7742D2x1y1 = 0;
590 double p7742D2x1y1xval[] = { -1.3, -1.1, -0.9, -0.7, -0.5, -0.3, -0.1, 0.1, 0.3,
507687cd 591 0.5, 0.7, 0.9, 1.1, 1.3 };
7c1a1f1d 592 double p7742D2x1y1xerrminus[] = { 0.09999999999999987, 0.09999999999999987, 0.09999999999999998, 0.10000000000000009, 0.09999999999999998, 0.10000000000000003, 0.1, 0.1, 0.09999999999999998,
507687cd 593 0.09999999999999998, 0.09999999999999998, 0.09999999999999998, 0.10000000000000009, 0.10000000000000009 };
7c1a1f1d 594 double p7742D2x1y1xerrplus[] = { 0.10000000000000009, 0.10000000000000009, 0.09999999999999998, 0.09999999999999998, 0.09999999999999998, 0.09999999999999998, 0.1, 0.1, 0.10000000000000003,
507687cd 595 0.09999999999999998, 0.10000000000000009, 0.09999999999999998, 0.09999999999999987, 0.09999999999999987 };
7c1a1f1d 596 double p7742D2x1y1yval[] = { 3.9, 3.89, 3.81, 3.7, 3.64, 3.59, 3.53, 3.58, 3.59,
507687cd 597 3.61, 3.74, 3.8, 3.87, 3.95 };
7c1a1f1d 598 double p7742D2x1y1yerrminus[] = { 0.13341664064126335, 0.13152946437965907, 0.13152946437965907, 0.1216552506059644, 0.1216552506059644, 0.1216552506059644, 0.1216552506059644, 0.1216552506059644, 0.1216552506059644,
507687cd 599 0.1216552506059644, 0.1216552506059644, 0.13152946437965907, 0.13152946437965907, 0.13341664064126335 };
7c1a1f1d 600 double p7742D2x1y1yerrplus[] = { 0.13341664064126335, 0.13152946437965907, 0.13152946437965907, 0.1216552506059644, 0.1216552506059644, 0.1216552506059644, 0.1216552506059644, 0.1216552506059644, 0.1216552506059644,
507687cd 601 0.1216552506059644, 0.1216552506059644, 0.13152946437965907, 0.13152946437965907, 0.13341664064126335 };
7c1a1f1d 602 int p7742D2x1y1numpoints = 14;
507687cd 603
7c1a1f1d 604 p7742D2x1y1 = new TGraphAsymmErrors(p7742D2x1y1numpoints, p7742D2x1y1xval, p7742D2x1y1yval, p7742D2x1y1xerrminus, p7742D2x1y1xerrplus, p7742D2x1y1yerrminus, p7742D2x1y1yerrplus);
605 p7742D2x1y1->SetName("/HepData/7742/d2x1y1");
606 p7742D2x1y1->SetTitle("/HepData/7742/d2x1y1");
607 p7742D2x1y1->SetMarkerStyle(21);
608 p7742D2x1y1->SetMarkerColor(kRed);
609 //p7742D2x1y1.Draw("AP");
507687cd 610
611
541c19ed 612
613
614
615
616 //End official data
617
507687cd 618 // CMS published NSD data
619
7c1a1f1d 620 TGraphAsymmErrors* p7743D8x1y1 = 0;
621 double p7743D8x1y1xval[] = { -2.25, -1.75, -1.25, -0.75, -0.25, 0.25, 0.75, 1.25, 1.75,
507687cd 622 2.25 };
7c1a1f1d 623 double p7743D8x1y1xerrminus[] = { 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25 };
624 double p7743D8x1y1xerrplus[] = { 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25 };
625 double p7743D8x1y1yval[] = { 3.6, 3.73, 3.62, 3.54, 3.48, 3.48, 3.54, 3.62, 3.73, 3.6 };
626 double p7743D8x1y1yerrminus[] = { 0.13, 0.14, 0.13, 0.13, 0.13, 0.13, 0.13, 0.13, 0.14,0.13 };
627 double p7743D8x1y1yerrplus[] = { 0.13, 0.14, 0.13, 0.13, 0.13, 0.13, 0.13, 0.13, 0.14, 0.13 };
628 int p7743D8x1y1numpoints = 10;
629 p7743D8x1y1 = new TGraphAsymmErrors(p7743D8x1y1numpoints, p7743D8x1y1xval, p7743D8x1y1yval, p7743D8x1y1xerrminus, p7743D8x1y1xerrplus, p7743D8x1y1yerrminus, p7743D8x1y1yerrplus);
507687cd 630
7c1a1f1d 631 p7743D8x1y1->SetMarkerStyle(20);
632 p7743D8x1y1->SetMarkerColor(kBlack);
541c19ed 633
507687cd 634 //End CMS
541c19ed 635
636
637 TH1I* hEvents = (TH1I*)fList->FindObject(fEvents.Data());
638
639 TH1I* hMCEvents = (TH1I*)fList->FindObject(fPrimEvents.Data());
640
641 //SPD part
642 TH1D* hSPDana = (TH1D*)fList->FindObject(Form("dNdeta_SPD_vtxbin%d_proj",5));
643
644 for(Int_t nn = 0; nn < pars->GetNvtxBins() ; nn++) {
645 TH1D* hSPDanalysis = (TH1D*)fList->FindObject(Form("dNdeta_SPD_vtxbin%d_proj",nn));
646 TH1D* hSPDanalysisTrVtx = (TH1D*)fList->FindObject(Form("dNdetaTrVtx_SPD_vtxbin%d_proj",nn));
647 TH1D* hSPDanalysisNSD = (TH1D*)fList->FindObject(Form("dNdetaNSD_SPD_vtxbin%d_proj",nn));
648
649 Float_t nEventsSPD = (Float_t)hEvents->GetBinContent(nn+1);
507687cd 650 if(!nEventsSPD) continue;
541c19ed 651 hSPDanalysis->Scale(1/nEventsSPD);
652 hSPDanalysisTrVtx->Scale(1/nEventsSPD);
653 hSPDanalysisNSD->Scale(1/nEventsSPD);
654 }
df2a9c32 655 Float_t lowlimits[10] = {-1,-1.25,-1.5,-1.75,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9};
656 Float_t highlimits[10] = {1.9,1.9,1.9,1.9,1.9,1.9,1.75,1.5,1.25,1};
657
541c19ed 658 TH1F* hSPDcombi = new TH1F("SPDcombi","SPDcombi",hSPDana->GetNbinsX(),hSPDana->GetXaxis()->GetXmin(),hSPDana->GetXaxis()->GetXmax());
659 TH1D* hSPDanalysis = 0;
660 for(Int_t kk = 1; kk <=hSPDana->GetNbinsX(); kk++) {
df2a9c32 661
541c19ed 662 Float_t weight = 0, wav=0,sumofw = 0;
df2a9c32 663
664 // if(hSPDana->GetBinCenter(kk) < lowlimits[nn])
665 // continue;
666 //if(hSPDana->GetBinCenter(kk) > highlimits[nn])
667 // continue;
668
669 for(Int_t nn =0; nn < pars->GetNvtxBins() ; nn++) {
670 Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
671 Float_t vtxZ1 = (delta*nn) - pars->GetVtxCutZ();
672 Float_t vtxZ2 = (delta*(nn+1)) - pars->GetVtxCutZ();
673
674 if(vtxZ1<fVtxCut1 || vtxZ2 >fVtxCut2)
675 continue;
676
677 //std::cout<<vtxZ1<<" "<<vtxZ2<<std::endl;
678
541c19ed 679 if(what == kMult)
680 hSPDanalysis = (TH1D*)fList->FindObject(Form("dNdeta_SPD_vtxbin%d_proj",nn));
681 else if(what == kMultNSD)
682 hSPDanalysis = (TH1D*)fList->FindObject(Form("dNdetaNSD_SPD_vtxbin%d_proj",nn));
683 else if(what == kMultTrVtx)
684 hSPDanalysis = (TH1D*)fList->FindObject(Form("dNdetaTrVtx_SPD_vtxbin%d_proj",nn));
685 else continue;
df2a9c32 686 if(hSPDanalysis->GetBinCenter(kk) < lowlimits[nn])
687 continue;
688 if(hSPDanalysis->GetBinCenter(kk) > highlimits[nn])
541c19ed 689 continue;
df2a9c32 690 //std::cout<<hSPDanalysis->GetBinCenter(kk)<<" "<<lowlimits[nn]<<std::endl;
691
541c19ed 692 Float_t mult = hSPDanalysis->GetBinContent(kk);
693 Float_t error = hSPDanalysis->GetBinError(kk);
df2a9c32 694 /*
541c19ed 695 if(mult > 0 && hSPDanalysis->GetBinContent(kk-1) < SMALLNUMBER)
696 continue;
697 if(mult > 0 && hSPDanalysis->GetBinContent(kk-2) < SMALLNUMBER)
698 continue;
699
700 if(mult > 0 && hSPDanalysis->GetBinContent(kk+1) < SMALLNUMBER)
701 continue;
702 if(mult > 0 && hSPDanalysis->GetBinContent(kk+2) < SMALLNUMBER)
703 continue;
df2a9c32 704 */
541c19ed 705 if(mult > 0) {
706 weight = 1/TMath::Power(error,2);
707 wav = wav + weight*mult;
708 sumofw = sumofw + weight;
709 }
710
711 }
712
713 if(sumofw && wav) {
714 Float_t errorTotal = 1/TMath::Sqrt(sumofw);
715 hSPDcombi->SetBinContent(kk,wav/sumofw);
716 hSPDcombi->SetBinError(kk,errorTotal);
717 }
718 }
719
720
721 Float_t xb1 = hSPDcombi->GetNbinsX();
722 Float_t xr1 = hSPDcombi->GetXaxis()->GetXmax() - hSPDcombi->GetXaxis()->GetXmin();
723 hSPDcombi->Scale(xb1 / xr1);
724 //hSPDcombi->Rebin(rebin);
725 //hSPDcombi->Scale(1/(Float_t)rebin);
726 //RebinHistogram(hSPDcombi,rebin);
727 hSPDcombi->SetMarkerStyle(29);
728 hSPDcombi->SetMarkerColor(kBlue);
729 //hSPDcombi->DrawCopy("same");
730
731
c5058a61 732 TCanvas* c1 = new TCanvas("cvertexbins","Overlaps",1400,800);
733 c1->Divide(5,2);
541c19ed 734 TCanvas* cCorrections = 0;
735 TCanvas* cCorrectionsPhi = 0;
059c7c6b 736 if(fDrawAll) {
737 cCorrections = new TCanvas("corrections","Correction vs Eta",1400,800);
738 cCorrections->Divide(5,2);
739
740 cCorrectionsPhi = new TCanvas("correctionsPhi","Correction vs Phi",1400,800);
741 cCorrectionsPhi->Divide(5,2);
742 }
743 //TCanvas* cphi = new TCanvas("dNdphi","dNdphi",1280,1024);
744 // cphi->Divide(3,2);
745
c5058a61 746 Int_t nVertexBins = pars->GetNvtxBins();
747
c5058a61 748
059c7c6b 749 //Int_t npadphi = 0;
c5058a61 750 TH1F* hPrimVtx = 0;
751 for(Int_t det = 1; det<=3; det++) {
752 Int_t maxRing = (det == 1 ? 0 : 1);
753 for(Int_t ring = 0;ring<=maxRing;ring++) {
059c7c6b 754 // npadphi++;
755
c5058a61 756 for(Int_t v=0; v< nVertexBins; v++) {
757 Char_t ringChar = (ring == 0 ? 'I' : 'O');
758 Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
759 Float_t vtxZ1 = (delta*v) - pars->GetVtxCutZ();
760 Float_t vtxZ2 = (delta*(v+1)) - pars->GetVtxCutZ();
761
2dcdd4eb 762 if(vtxZ1<fVtxCut1 || vtxZ2 >fVtxCut2)
c5058a61 763 continue;
764 Float_t nEvents = (Float_t)hEvents->GetBinContent(v+1);
059c7c6b 765 if(fDrawAll) {
766 TH2F* hBgCor = pars->GetBackgroundCorrection(det,ringChar,v);
767 if(hBgCor) {
768 TH1D* hBgProj = hBgCor->ProjectionX(Form("hBgProj_FMD%d%c_vtxbin%d",det,ringChar,v));
769 TH1D* hBgProjPhi = hBgCor->ProjectionY(Form("hBgProjPhi_FMD%d%c_vtxbin%d",det,ringChar,v));
770
771 hBgProjPhi->SetName(Form("FMD%d%c",det,ringChar));
772 hBgProjPhi->SetTitle("");//Form("FMD%d",det));
773 hBgProj->SetTitle("");
774 //Float_t scalefactor = (hBgProj->GetXaxis()->GetXmax() - hBgProj->GetXaxis()->GetXmin()) / (Float_t)hBgProj->GetNbinsX();
775 Float_t scalefactor = 1/(Float_t)hBgCor->GetNbinsY();
776
777 Float_t nFilledEtaBins = 0;
778
779 for(Int_t jj = 1; jj<=hBgProj->GetNbinsX(); jj++) {
780 if(hBgProj->GetBinContent(jj) > 0.01)
781 nFilledEtaBins++;
782 }
783
784 Float_t scalefactorPhi = 1/nFilledEtaBins;
785
786
787 hBgProj->Scale(scalefactor);
788 cCorrections->cd(v+1);
789
790 hBgProj->SetMarkerColor(det + 2*ring);
791 hBgProj->SetMarkerStyle(19 + det + 2*ring );
792
793 if((det + 2*ring) == 5) hBgProj->SetMarkerColor(det + 2*ring+1);
794 if((19 + det + 2*ring ) == 24) hBgProj->SetMarkerStyle(29 );
795
796 if(det == 1) {
797 hBgProj->GetYaxis()->SetRangeUser(0,4);
798 hBgProj->SetStats(kFALSE);
799 hBgProj->SetXTitle("#eta");
800 hBgProj->DrawCopy("PE");
801 TLatex* l = new TLatex(0.14,0.92,Form("Vtx range [%.1f, %.1f]",vtxZ1,vtxZ2));
802 l->SetNDC(kTRUE);
803 l->Draw();
804 }
805 else
806 hBgProj->DrawCopy("PEsame");
807
808 hBgProjPhi->Scale(scalefactorPhi);
809 cCorrectionsPhi->cd(v+1);
810 hBgProjPhi->SetMarkerColor(det + 2*ring);
811 hBgProjPhi->SetMarkerStyle(19 + det + 2*ring );
812 if((det + 2*ring) == 5) hBgProjPhi->SetMarkerColor(det + 2*ring+1);
813 if((19 + det + 2*ring ) == 24) hBgProjPhi->SetMarkerStyle(29 );
814 if(det == 1) {
815 hBgProjPhi->GetYaxis()->SetRangeUser(1,5);
816 hBgProjPhi->SetStats(kFALSE);
817
818 hBgProjPhi->SetXTitle("#Phi");
819 hBgProjPhi->DrawCopy("PE");
820
821
822 }
823 else
824 hBgProjPhi->DrawCopy("PEsame");
825
826 }
827
828 }
829
c5058a61 830 TH1F* multhistproj = (TH1F*)fList->FindObject(GetAnalysisName(what,det,ringChar,v));
831
832 c1->cd(v+1);
833
834 if(det==1) {
835 multhistproj->SetMarkerStyle(20);
836 multhistproj->SetMarkerColor(1);
837 }
838 if(det==2 && ringChar=='I') {
839 multhistproj->SetMarkerStyle(21);
840 multhistproj->SetMarkerColor(2);
841 }
842 if(det==2 && ringChar=='O') {
843 multhistproj->SetMarkerStyle(3);
844 multhistproj->SetMarkerColor(3);
845 }
846 if(det==3 && ringChar=='I') {
847 multhistproj->SetMarkerStyle(22);
848 multhistproj->SetMarkerColor(4);
849 }
850 if(det==3 && ringChar=='O') {
851 multhistproj->SetMarkerStyle(23);
852 multhistproj->SetMarkerColor(6);
853 }
854
855
856
857 if(det == 1) {
2dcdd4eb 858 if(what == kHits || what == kHitsTrVtx)
df2a9c32 859 hPrimVtx = (TH1F*)fMultList[what]->FindObject(Form("hMCHits_vtxbin%d_%s",v,fAnalysisNames[what].Data()));
c5058a61 860 else
861 hPrimVtx = (TH1F*)fList->FindObject(GetPrimName(what,det,ringChar,v));
059c7c6b 862 // std::cout<<hPrimVtx<<" "<<kHits<<" "<<kHitsTrVtx<<" "<<what<<std::endl;
df2a9c32 863 //std::cout<<Form("hMCHits_vtxbin%d_%s",v,fAnalysisNames[what].Data())<<std::endl;
c5058a61 864 hPrimVtx->SetTitle("");
865 TLatex* l = new TLatex(0.14,0.92,Form("Vtx range [%.1f, %.1f], %d events",vtxZ1,vtxZ2,(Int_t)nEvents));
866 l->SetNDC(kTRUE);
867 hPrimVtx->SetFillColor(kGray);
868 hPrimVtx->SetLabelFont(132,"xyz");
869 hPrimVtx->SetStats(0000);
870 hPrimVtx->GetXaxis()->SetTitle("#eta");
059c7c6b 871 hPrimVtx->GetYaxis()->SetRangeUser(0,6);
c5058a61 872 hPrimVtx->DrawCopy("E3");
873 l->Draw();
874 multhistproj->DrawCopy("same");
df2a9c32 875 TH1D* hSPDanavtxbin = 0;
032f0117 876 if(what == kMult || what == kMultTrVtx || what == kMultNSD) {
877
df2a9c32 878 if(what == kMult)
879 hSPDanavtxbin = (TH1D*)fList->FindObject(Form("dNdeta_SPD_vtxbin%d_proj",v));
880 if(what == kMultTrVtx)
881 hSPDanavtxbin = (TH1D*)fList->FindObject(Form("dNdetaTrVtx_SPD_vtxbin%d_proj",v));
541c19ed 882
df2a9c32 883 if(what == kMultNSD)
884 hSPDanavtxbin = (TH1D*)fList->FindObject(Form("dNdetaNSD_SPD_vtxbin%d_proj",v));
032f0117 885
541c19ed 886 hSPDanavtxbin->SetMarkerColor(kBlue);
887 hSPDanavtxbin->SetMarkerStyle(kStar);
888 hSPDanavtxbin->Scale(xb1 / xr1);
889 hSPDanavtxbin->DrawCopy("same");
890
059c7c6b 891 if(what != kMultNSD) {
892 graphinel->Draw("sameP");
893 graphinel2->Draw("sameP");
894 }
895 else
896 {
897 graph->Draw("sameP");
898 graph2->Draw("sameP");
899 }
032f0117 900 }
c5058a61 901 }
902 else
903 multhistproj->DrawCopy("same");
904
032f0117 905
c5058a61 906 }
907 }
908 }
909
059c7c6b 910 //Legends for corrections
911 if(fDrawAll) {
032f0117 912 for(Int_t v=0; v< nVertexBins; v++) {
913 TPad* pad= (TPad*)cCorrectionsPhi->cd(v+1);
914 pad->BuildLegend(0.15,0.45,0.45,0.9);
915 Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
916 Float_t vtxZ1 = (delta*v) - pars->GetVtxCutZ();
917 Float_t vtxZ2 = (delta*(v+1)) - pars->GetVtxCutZ();
918 TLatex* l = new TLatex(0.14,0.92,Form("Vtx range [%.1f, %.1f]",vtxZ1,vtxZ2));
919 l->SetNDC(kTRUE);
920 l->Draw();
921 }
059c7c6b 922 }
c5058a61 923 for(Int_t v=0; v< nVertexBins; v++) {
df2a9c32 924 TH1F* sumMultHist = (TH1F*)fMultList[what]->FindObject(Form("hMult_vtxbin%d_%s",v,fAnalysisNames[what].Data()));
c5058a61 925 c1->cd(v+1);
926 sumMultHist->SetMarkerStyle(25);
927 sumMultHist->DrawCopy("same");
928
929 }
df2a9c32 930 TH1F* primMult = (TH1F*)fMultList[what]->FindObject(Form("primary_%s",fAnalysisNames[what].Data()));
931 TH1F* sumMult = (TH1F*)fMultList[what]->FindObject(Form("dNdeta_%s",fAnalysisNames[what].Data()));
c5058a61 932 sumMult->SetMarkerStyle(23);
933 sumMult->SetMarkerColor(kRed);
934 for(Int_t det = 1; det<=3;det++) {
935 Int_t maxRing = (det == 1 ? 0 : 1);
936 for(Int_t iring = 0; iring<=maxRing; iring++) {
937 Char_t ringChar = (iring == 0 ? 'I' : 'O');
df2a9c32 938 TH1F* hRingMult= (TH1F*)fMultList[what]->FindObject(Form("hRingMult_FMD%d%c_%s",det,ringChar,fAnalysisNames[what].Data()));
c5058a61 939 hRingMult->Add(primMult,-1);
940 hRingMult->Divide(primMult);
941 }
942 }
943
032f0117 944
059c7c6b 945 // End of SPD
032f0117 946
c5058a61 947
948 //TH1F* hPrim = (TH1F*)fList->FindObject(fPrimEvents.Data());
949 // TFile testhit("/home/canute/ALICE/Simulations/TestOfAnalysis2/hitsdist.root");
950 // TH1F* hHitCor = (TH1F*)testhit.Get("selectedHits");
951 // hHitCor->Sumw2();
952 // sumMult->Divide(hHitCor);
953 TH1F* hPrim = (TH1F*)fList->FindObject(fPrimdNdeta.Data());
954 Float_t nPrimevents = hMCEvents->GetEntries();
955 //std::cout<<hMCEvents->GetEntries()<<std::endl;
956 Float_t xb = hPrim->GetNbinsX();
957 Float_t xr = hPrim->GetXaxis()->GetXmax() - hPrim->GetXaxis()->GetXmin();
958 //std::cout<<xb/xr<<std::endl;
959 hPrim->Scale(xb / xr );
960 if(nPrimevents > 0)
961 hPrim->Scale(1/nPrimevents);
962
2dcdd4eb 963
964 // primMult->Rebin(rebin);
965 // primMult->Scale(1/rebin);
966 // hPrim->Rebin(rebin);
967 // hPrim->Scale(1/rebin);
968 if(what != kHits && what != kHitsTrVtx)
c5058a61 969 primMult = hPrim;
70d74659 970
c5058a61 971 // hReldif->Add(hPrim,-1);
972 // hReldif->Divide(hPrim);
973
059c7c6b 974
975
976 /////////////*******************New thing!!!
977
978 //TH3D* hist3d = (TH3D*)fDataObject->ProjectionXYZ("hist3d");
979 // fDataObject->Sumw2();
980 //TH2F* projeta = (TH2F*)fDataObject->Project3D("yx");
981
982 // TProfile2D* projeta = fDataObject->Project3DProfile("yx");
983 // projeta->SetSumw2();
984 //TProfile* etaprofile = projeta->ProfileX("dNdeta_profile");
985 //TProfile* etaprofile = projeta->ProfileX("dNdeta_profile");
986 // TH1* etaprofile = projeta->ProjectionX("dNdeta_profile");
987 //TProfile* etaprofile = fDataObject->ProfileX();
988 /*
989 TCanvas* ctest = new TCanvas("test","test",1200,600);
990 ctest->Divide(3);
991 ctest->cd(1);
992 fDataObject->DrawCopy("box");
993 ctest->cd(2);
994 projeta->DrawCopy("colz");
995 ctest->cd(3);
996 etaprofile->DrawCopy();
997 */
70d74659 998 //sumMult->Rebin(rebin);
999 TH1F* hReldif = (TH1F*)sumMult->Clone("hReldif");
1000 if(rebin != 1) {
1001 RebinHistogram(sumMult,rebin);
ba77343d 1002 if(!realdata) {
1003 RebinHistogram(primMult,rebin);
1004 RebinHistogram(hReldif,rebin);
1005 }
70d74659 1006 }
1007 hReldif->Add(primMult,-1);
1008 hReldif->Divide(primMult);
c5058a61 1009 TCanvas* c2 = new TCanvas("dN/deta","dN/deta",640,960);
1010 c2->Divide(1,2);//,0,0);
1011 c2->cd(2);
2dcdd4eb 1012 gPad->Divide(1,2);
1013 gPad->cd(1);
c5058a61 1014 gPad->SetGridy();
1015 gPad->SetGridx();
1016 hReldif->SetTitle("");
1017 hReldif->GetYaxis()->SetTitle("Relative difference");
1018 hReldif->GetYaxis()->SetRangeUser(-0.2,0.2);
1019 hReldif->SetStats(0000);
1020 hReldif->GetXaxis()->SetTitle("#eta");
ba77343d 1021
c5058a61 1022 hReldif->DrawCopy();
059c7c6b 1023 //SPD Rel dif
1024 TH1F* hReldifSPD = (TH1F*)hSPDcombi->Clone("hReldifSPD");
1025 if(rebin != 1) {
1026 RebinHistogram(hSPDcombi,rebin);
1027 if(!realdata) {
1028 RebinHistogram(hReldifSPD,rebin);
1029 }
1030 }
1031 hReldifSPD->Add(primMult,-1);
1032 hReldifSPD->Divide(primMult);
1033 hReldifSPD->DrawCopy("same");
1034
1035 TLine* zeroLine = new TLine(-4,0,6,0);
c5058a61 1036 zeroLine->SetLineWidth(2);
1037 zeroLine->SetLineColor(kBlack);
1038 zeroLine->Draw();
1039 hPrim->SetTitle("");
1040 hPrim->SetLabelFont(132,"xyz");
1041 hPrim->SetFillColor(kGray);
1042 primMult->SetFillColor(kBlue);
70d74659 1043
1044
c5058a61 1045 c2->cd(1);
1046 gPad->SetGridy();
1047 gPad->SetGridx();
1048 hPrim->GetYaxis()->SetRangeUser(2,10);
1049 //hPrim->GetYaxis()->SetRangeUser(0,20);
1050 hPrim->SetStats(0000);
1051 hPrim->GetXaxis()->SetTitle("#eta");
ae26bdd7 1052 sumMult->SetTitle("");
1053 sumMult->SetStats(kFALSE);
059c7c6b 1054 sumMult->SetXTitle("#eta");
1055 sumMult->SetYTitle("#frac{dN}{d#eta_{ch}}");
70d74659 1056 // sumMult->GetYaxis()->SetRangeUser
2dcdd4eb 1057 //sumMult->Scale(1/(Float_t)rebin);
c5058a61 1058 sumMult->DrawCopy("PE");
059c7c6b 1059 //Syst errors
1060 TH1F* sumMultSystPos = (TH1F*)sumMult->Clone("systerrorsPos");
1061 TH1F* sumMultSystNeg = (TH1F*)sumMult->Clone("systerrorsNeg");
1062 for(Int_t jj = 1;jj <=sumMultSystPos->GetNbinsX();jj++) {
1063 if(sumMultSystPos->GetBinCenter(jj) < 0) {
1064 sumMultSystPos->SetBinError(jj,0);
1065 sumMultSystPos->SetBinContent(jj,0);
1066 continue;
1067 }
1068
1069 if(sumMultSystPos->GetBinContent(jj) > 0)
1070 sumMultSystPos->SetBinError(jj,TMath::Sqrt(TMath::Power(sumMultSystPos->GetBinError(jj),2)+TMath::Power(0.1*sumMultSystPos->GetBinContent(jj),2)));
1071 }
1072 for(Int_t jj = 1;jj <=sumMultSystNeg->GetNbinsX();jj++) {
1073 if(sumMultSystNeg->GetBinCenter(jj) > 0) {
1074 sumMultSystNeg->SetBinError(jj,0);
1075 sumMultSystNeg->SetBinContent(jj,0);
1076 continue;
1077 }
1078
1079 if(sumMultSystNeg->GetBinContent(jj) > 0)
1080 sumMultSystNeg->SetBinError(jj,TMath::Sqrt(TMath::Power(sumMultSystNeg->GetBinError(jj),2)+TMath::Power(0.1*sumMultSystNeg->GetBinContent(jj),2)));
1081 }
1082
1083
1084
1085 sumMultSystPos->SetFillColor(18);
1086 sumMultSystNeg->SetFillColor(18);
541c19ed 1087 //sumMultSystPos->DrawCopy("sameE5");
1088 //sumMultSystNeg->DrawCopy("sameE5");
059c7c6b 1089 //End syst errors
1090
1091
1092
2dcdd4eb 1093 if(what != kHits && what != kHitsTrVtx && hPrim->GetEntries())
ae26bdd7 1094 hPrim->DrawCopy("E3same");
2dcdd4eb 1095 if(what == kHits || what == kHitsTrVtx)
ae26bdd7 1096 primMult->DrawCopy("E3same");
c5058a61 1097 hPrim->GetXaxis()->SetTitle("#eta");
1098
2dcdd4eb 1099 TH1F* hMirror = new TH1F("hMirror","mirrored",sumMult->GetNbinsX(),sumMult->GetXaxis()->GetXmin(),sumMult->GetXaxis()->GetXmax());
1100
65f0c9c8 1101 for(Int_t i= (Int_t)0.5*sumMult->GetNbinsX(); i<=sumMult->GetNbinsX(); i++) {
2dcdd4eb 1102 Float_t eta = sumMult->GetBinCenter(i);
1103 Float_t value = sumMult->GetBinContent(i);
1104 Float_t error = sumMult->GetBinError(i);
1105 if(value > 0) {
1106 hMirror->SetBinContent(hMirror->FindBin(-1*eta),value);
1107 hMirror->SetBinError(hMirror->FindBin(-1*eta),error);
1108 }
1109 }
1110 TH1F* hReldifLR = (TH1F*)sumMult->Clone("hReldifLeftRight");
1111 hReldifLR->Add(hMirror,-1);
1112 hReldifLR->Divide(hMirror);
1113 c2->cd(2);
1114 gPad->cd(1);
059c7c6b 1115 hReldifLR->GetYaxis()->SetRangeUser(-0.2,0.2);
1116
1117 hReldifLR->SetYTitle("Relative difference left-right");
1118 hReldifLR->SetLabelSize(2*hReldifLR->GetLabelSize("X"),"X");
1119 hReldifLR->SetLabelSize(2*hReldifLR->GetLabelSize("Y"),"Y");
1120 hReldifLR->SetTitleSize(2*hReldifLR->GetTitleSize("X"),"X");
1121 hReldifLR->SetTitleSize(2*hReldifLR->GetTitleSize("Y"),"Y");
1122 hReldifLR->SetTitleOffset(0.7*hReldifLR->GetTitleOffset("X"),"X");
1123 hReldifLR->SetTitleOffset(0.5*hReldifLR->GetTitleOffset("Y"),"Y");
1124
1125
1126 if(realdata)
1127 hReldifLR->DrawCopy();
1128 zeroLine->Draw();
2dcdd4eb 1129 c2->cd(1);
1130 hMirror->SetMarkerStyle(25);
1131
c5058a61 1132 hPrim->SetMarkerStyle(8);
2dcdd4eb 1133
2dcdd4eb 1134 //graph2->Draw("P");
1135 TH1F* hPythiaMB = 0;
059c7c6b 1136
541c19ed 1137 TFile* fpyt = 0;
4b8bdb60 1138
541c19ed 1139 if(!pythiafile.Contains("none"))
1140 fpyt = TFile::Open(pythiafile);
4b8bdb60 1141
059c7c6b 1142 if(realdata ) {
541c19ed 1143 if(fpyt && !pythiafile.Contains("none")) {
059c7c6b 1144 hPythiaMB = (TH1F*)fpyt->Get("hPythiaMB");
1145 }
1146 else
1147 std::cout<<"no pythia for this energy"<<std::endl;
2dcdd4eb 1148 }
1149 if(what != kHits && what != kHitsTrVtx) {
1150 if(hPrim->GetEntries() != 0 && !realdata)
1151 hPrim->DrawCopy("PE3same");
1152 if(realdata) {
1153 //graph->Draw("PEsame");
1154 //graph2->Draw("PEsame");
059c7c6b 1155 if(what != kMultNSD) {
1156 graphinel->Draw("PEsame");
1157 graphinel2->Draw("PEsame");
7c1a1f1d 1158 p7742D1x1y1->Draw("PEsame");
059c7c6b 1159 }
1160 else{
1161 graph->Draw("PEsame");
1162 graph2->Draw("PEsame");
7c1a1f1d 1163 p7742D2x1y1->Draw("PEsame");
1164 p7743D8x1y1->Draw("PEsame");
059c7c6b 1165
1166 }
1167
2dcdd4eb 1168 }
1169
1170 }
2dcdd4eb 1171
c5058a61 1172
059c7c6b 1173 sumMult->DrawCopy("PEsame");
c5058a61 1174
1175
059c7c6b 1176
1177 hSPDcombi->DrawCopy("same");
1178 if(realdata)
1179 hMirror->DrawCopy("PEsame");
1180
1181 TLegend* leg = new TLegend(0.3,0.2,0.7,0.45,"");
2dcdd4eb 1182 if(!realdata) {
1183 if(what != kHits && what != kHitsTrVtx)
1184 leg->AddEntry(hPrim,"Primary","pf");
1185 else
1186 leg->AddEntry(primMult,"Hits","pf");
1187 }
059c7c6b 1188
1189 if(what == kMult)
1190 leg->AddEntry(sumMult,"FMD INEL","p");
1191 else if(what == kMultTrVtx)
1192 leg->AddEntry(sumMult,"FMD TrVtx","p");
1193 else if(what == kMultNSD)
1194 leg->AddEntry(sumMult,"FMD NSD","p");
2dcdd4eb 1195
1196
2dcdd4eb 1197 if(realdata) {
059c7c6b 1198 if(what == kMult)
541c19ed 1199 leg->AddEntry(hMirror,"Mirror FMD INEL","p");
059c7c6b 1200 else if(what == kMultTrVtx)
541c19ed 1201 leg->AddEntry(hMirror,"Mirror FMD TrVtx","p");
059c7c6b 1202 else if(what == kMultNSD)
541c19ed 1203 leg->AddEntry(hMirror,"Mirror FMD NSD","p");
059c7c6b 1204
1205 if(what == kMultNSD) {
1206 leg->AddEntry(graph,"UA5 NSD","p");
1207 leg->AddEntry(graph2,"Mirror UA5 NSD","p"); }
1208 else if(what == kMult) {
2dcdd4eb 1209 leg->AddEntry(graphinel,"UA5 INEL","p");
1210 leg->AddEntry(graphinel2,"Mirror UA5 INEL","p");
059c7c6b 1211 }
1212 //leg->AddEntry(hPythiaMB,"Pythia MB","l");
541c19ed 1213 leg->AddEntry(hSPDcombi,"HHD SPD clusters","p");
507687cd 1214 if(what == kMult)
7c1a1f1d 1215 leg->AddEntry(p7742D1x1y1,"ALICE INEL published","p");
507687cd 1216 if(what == kMultNSD) {
7c1a1f1d 1217 leg->AddEntry(p7742D2x1y1,"ALICE NSD published","p");
1218 leg->AddEntry(p7743D8x1y1, "CMS NSD published","p");
507687cd 1219 }
2dcdd4eb 1220 }
c5058a61 1221 leg->Draw();
2dcdd4eb 1222
1223 c2->cd(2);
1224 gPad->cd(2);
1225 TH1F* hRatioMultPythia = 0;
1226 TH1F* hRatioMultUA5 = 0;
541c19ed 1227 TH1F* hRatioMultUA5_rebin5 = new TH1F("hRatioMultUA5_rebin5","hRatioMultUA5_rebin5",sumMult->GetNbinsX(),sumMult->GetXaxis()->GetXmin(),sumMult->GetXaxis()->GetXmax());
1228
1229 Double_t xval=0, yval=0;
1230 Int_t ipos=0;
1231 for(Int_t bb =hRatioMultUA5_rebin5->FindBin(0); bb<=hRatioMultUA5_rebin5->GetNbinsX();bb++) {
1232
1233
1234
1235 xval=yval=0;
1236
1237 if(hRatioMultUA5_rebin5->GetBinCenter(bb) >= 0) {
507687cd 1238 if(what == kMultNSD)
1239 graph->GetPoint(ipos,xval,yval);
1240 else
1241 graphinel->GetPoint(ipos,xval,yval);
541c19ed 1242
1243 if(yval>0) {
1244 hRatioMultUA5_rebin5->SetBinContent(bb,yval);
1245 hRatioMultUA5_rebin5->SetBinError(bb,graphinel->GetErrorY(ipos));
1246 if(hRatioMultUA5_rebin5->GetBinCenter(bb) < 4) {
1247 hRatioMultUA5_rebin5->SetBinContent(hRatioMultUA5_rebin5->FindBin(-1*hRatioMultUA5_rebin5->GetBinCenter(bb)),yval);
507687cd 1248 hRatioMultUA5_rebin5->SetBinError(hRatioMultUA5_rebin5->FindBin(-1*hRatioMultUA5_rebin5->GetBinCenter(bb)),(what == kMultNSD ? graph->GetErrorY(ipos) : graph->GetErrorY(ipos)));
541c19ed 1249
1250 }
1251 ipos++;
1252 }
1253 }
1254 }
1255
1256
1257
2dcdd4eb 1258 if(realdata) {
1259
1260 Float_t ratio = 1;
541c19ed 1261 if(hPythiaMB) {
1262 if(hPythiaMB->GetNbinsX() != sumMult->GetNbinsX())
1263 ratio = (Float_t)sumMult->GetNbinsX() / (Float_t)hPythiaMB->GetNbinsX();
1264 }
1265
2dcdd4eb 1266 hRatioMultPythia = (TH1F*)sumMult->Clone("hRatioMultPythia");
1267 hRatioMultUA5 = (TH1F*)sumMult->Clone("hRatioMultUA5");
1268 if(ratio > 1) {
65f0c9c8 1269 hRatioMultPythia->Rebin((Int_t)ratio);
2dcdd4eb 1270 hRatioMultPythia->Scale(1/ratio);
1271 }
541c19ed 1272 if(ratio < 1 && hPythiaMB) {
65f0c9c8 1273 hPythiaMB->Rebin((Int_t)(1/ratio));
2dcdd4eb 1274 hPythiaMB->Scale(ratio);
1275 }
1276
541c19ed 1277 if(rebin !=5) {
507687cd 1278 TGraphErrors* tmp1;
1279 if(what == kMultNSD)
059c7c6b 1280 tmp1 = (TGraphErrors*)graph->Clone("UA5tmp");
1281 else
1282 tmp1 = (TGraphErrors*)graphinel->Clone("UA5tmp");
1283
1284 tmp1->Fit("pol8","Q0","Q0",1.5,6);
1285 TF1* hFit = tmp1->GetFunction("pol8");
1286 for(Int_t ii = hRatioMultUA5->GetNbinsX() / 2; ii<hRatioMultUA5->GetNbinsX();ii++) {
1287 if(hRatioMultUA5->GetBinContent(ii) > 0) {
1288 Float_t errorratio = hRatioMultUA5->GetBinError(ii) / hRatioMultUA5->GetBinContent(ii);
1289 hRatioMultUA5->SetBinContent(ii,
1290 hRatioMultUA5->GetBinContent(ii) /
1291 ((hRatioMultUA5->GetBinCenter(ii) > 4.8) ? hFit->Eval(hRatioMultUA5->GetBinCenter(ii-1)) : hFit->Eval(hRatioMultUA5->GetBinCenter(ii))));
1292 hRatioMultUA5->SetBinError(ii,hRatioMultUA5->GetBinContent(ii) * TMath::Sqrt(0.02*0.02 + TMath::Power(errorratio,2)));
1293
1294
1295 }
1296
1297 }
1298
1299 TGraphErrors* tmp2;
1300 if(what == kMultNSD)
1301 tmp2 = (TGraphErrors*)graph2->Clone("UA5tmp2");
1302 else
1303 tmp2 = (TGraphErrors*)graphinel2->Clone("UA5tmp2");
2dcdd4eb 1304
059c7c6b 1305 //tmp2->Fit("pol8","Q0","Q0",-3.7,-1.5);
1306 tmp2->Fit("pol7","Q0","Q0",-3.7,0);
1307 hFit = tmp2->GetFunction("pol7");
1308 for(Int_t ii = 0; ii<hRatioMultUA5->GetNbinsX()/2;ii++) {
1309 if(hRatioMultUA5->GetBinContent(ii) > 0) {
1310 Float_t errorratio = hRatioMultUA5->GetBinError(ii) / hRatioMultUA5->GetBinContent(ii);
1311 hRatioMultUA5->SetBinContent(ii,hRatioMultUA5->GetBinContent(ii) / hFit->Eval(hRatioMultUA5->GetBinCenter(ii)) );
1312 hRatioMultUA5->SetBinError(ii,hRatioMultUA5->GetBinContent(ii) * TMath::Sqrt(0.02*0.02 + TMath::Power(errorratio,2)));
1313
1314 }
1315
1316
1317 graphinel->GetHistogram()->SetStats(kFALSE);
1318 graphinel2->GetHistogram()->SetStats(kFALSE);
1319
1320 }
541c19ed 1321 }
1322 else hRatioMultUA5->Divide(hRatioMultUA5_rebin5);
059c7c6b 1323 //}
1324
1325
1326
1327 gPad->SetGridx();
1328 gPad->SetGridy();
1329 TLine* oneLine = new TLine(-4,1,6,1);
1330 oneLine->SetLineWidth(2);
2dcdd4eb 1331
2dcdd4eb 1332 hRatioMultUA5->SetMarkerColor(kGreen);
1333 hRatioMultUA5->SetMarkerStyle(20);
059c7c6b 1334 hRatioMultUA5->GetYaxis()->SetRangeUser(0.75,1.25);
1335
1336 hRatioMultUA5->SetYTitle("Ratio FMD/UA5");
1337 hRatioMultUA5->SetLabelSize(2*hRatioMultUA5->GetLabelSize("X"),"X");
1338 hRatioMultUA5->SetLabelSize(2*hRatioMultUA5->GetLabelSize("Y"),"Y");
1339 hRatioMultUA5->SetTitleSize(2*hRatioMultUA5->GetTitleSize("X"),"X");
1340 hRatioMultUA5->SetTitleSize(2*hRatioMultUA5->GetTitleSize("Y"),"Y");
1341 hRatioMultUA5->SetTitleOffset(0.7*hRatioMultUA5->GetTitleOffset("X"),"X");
1342 hRatioMultUA5->SetTitleOffset(0.5*hRatioMultUA5->GetTitleOffset("Y"),"Y");
1343
1344
1345
1346 hRatioMultUA5->DrawCopy();
1347 //hRatioMultPythia->DrawCopy("same");
2dcdd4eb 1348 oneLine->Draw("same");
1349 }
059c7c6b 1350
1351 TCanvas* cPaper = new TCanvas("FMD_SPD_UA5","FMD_SPD_UA5",800,600);
1352 cPaper->cd();
1353
1354 sumMult->DrawCopy();
541c19ed 1355 //sumMultSystPos->DrawCopy("sameE5");
1356 // sumMultSystNeg->DrawCopy("sameE5");
059c7c6b 1357 //hTestdNdeta->DrawCopy("same");
541c19ed 1358 //hSPDcombi->DrawCopy("same");
507687cd 1359 if(what == kMult)
7c1a1f1d 1360 p7742D1x1y1->Draw("PEsame");
507687cd 1361 if(what == kMultNSD) {
7c1a1f1d 1362 p7742D2x1y1->Draw("PEsame");
1363 p7743D8x1y1->Draw("PEsame");
507687cd 1364 }
1365
541c19ed 1366 TLegend* leg2 = new TLegend(0.3,0.2,0.7,0.45,"");
1367 if(what == kMult)
1368 leg2->AddEntry(sumMult,"FMD INEL","p");
1369 else if(what == kMultTrVtx)
1370 leg2->AddEntry(sumMult,"FMD TrVtx","p");
507687cd 1371 else if(what == kMultNSD)
541c19ed 1372 leg2->AddEntry(sumMult,"FMD NSD","p");
1373
1374 if(realdata) {
1375
1376 /* if(what == kMult)
1377 leg2->AddEntry(hMirror,"Mirror FMD INEL","p");
1378 else if(what == kMultTrVtx)
1379 leg2->AddEntry(hMirror,"FMD TrVtx","p");
1380 else if(what == kMultNSD)
1381 leg2->AddEntry(hMirror,"FMD NSD","p");*/
1382
1383 if(what == kMultNSD) {
1384 leg2->AddEntry(graph,"UA5 NSD","p");
1385 leg2->AddEntry(graph2,"Mirror UA5 NSD","p"); }
1386 else if(what == kMult) {
1387 leg2->AddEntry(graphinel,"UA5 INEL","p");
1388 leg2->AddEntry(graphinel2,"Mirror UA5 INEL","p");
1389 }
1390 //leg2->AddEntry(hPythiaMB,"Pythia MB","l");
1391 //leg2->AddEntry(hSPDcombi,"HHD SPD clusters","p");
507687cd 1392 if(what == kMult)
7c1a1f1d 1393 leg2->AddEntry(p7742D1x1y1,"ALICE INEL published","p");
507687cd 1394 if(what == kMultNSD) {
7c1a1f1d 1395 leg2->AddEntry(p7742D2x1y1,"ALICE NSD published","p");
1396 leg2->AddEntry(p7743D8x1y1,"CMS NSD published","p");
507687cd 1397 }
541c19ed 1398 }
1399 leg2->Draw();
1400
1401
1402
059c7c6b 1403 if(what != kMultNSD) {
1404 graphinel->Draw("PEsame");
1405 graphinel2->Draw("PEsame");
1406 }
1407 else {
1408 graph->Draw("PEsame");
1409 graph2->Draw("PEsame");
1410 }
1411
1412 sumMult->DrawCopy("PEsame");
2dcdd4eb 1413 c2->cd(1);
059c7c6b 1414
1415 c2->Print("fmdana.png");
c5058a61 1416 TString species;
1417
1418 switch(what) {
1419 case kMult:
541c19ed 1420 species = "mult";
c5058a61 1421 break;
1422 case kMultTrVtx:
541c19ed 1423 species = "mult_TrVtx";
c5058a61 1424 break;
1425 case kHits:
541c19ed 1426 species = "hits";
c5058a61 1427 break;
2dcdd4eb 1428 case kHitsTrVtx:
541c19ed 1429 species = "hits_TrVtx";
2dcdd4eb 1430 break;
059c7c6b 1431 case kMultNSD:
541c19ed 1432 species = "mult_NSD";
059c7c6b 1433 break;
1434
c5058a61 1435 default:
1436 AliWarning("no valid Analysis entry!");
1437 break;
1438 }
1439 TFile fout(Form("fmd_dNdeta_%s.root",species.Data()),"recreate");
2dcdd4eb 1440 if(hRatioMultPythia)
1441 hRatioMultPythia->Write();
1442 hPrim->Write();
059c7c6b 1443 graph->Write("UA5nsd");
1444 graph2->Write("UA5mirrornsd");
1445 graphinel->Write("UA5inel");
1446 graphinel2->Write("UA5mirrorinel");
c5058a61 1447 sumMult->Write();
059c7c6b 1448 hSPDcombi->Write();
c5058a61 1449 hReldif->Write();
70d74659 1450 fMultList[what]->Write();
2dcdd4eb 1451 c2->Write();
059c7c6b 1452 //fDataObject->Write();
c5058a61 1453 fout.Close();
1454
7e2bf482 1455}
1456//_____________________________________________________________________
1457void AliFMDDndeta::CreateSharingEfficiency(const Char_t* filename, Bool_t store) {
2dcdd4eb 1458 //Generate sharing efficiency
7e2bf482 1459 Init(filename);
1460 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
1461 //pars->Init(kTRUE,AliFMDAnaParameters::kBackgroundCorrection);
1462 Int_t nVertexBins = pars->GetNvtxBins();
1463
1464 SetNames(kHits);
3d4a1473 1465 // "nEvents";
7e2bf482 1466 TH1I* hEvents = (TH1I*)fList->FindObject(fEvents.Data());
3d4a1473 1467 // "nMCEventsNoCuts"
7e2bf482 1468 TH1I* hPrimEvents = (TH1I*)fList->FindObject(fPrimEvents.Data());
1469
1470 SetNames(kHitsTrVtx);
3d4a1473 1471 // "nEvents";
7e2bf482 1472 TH1I* hEventsTrVtx = (TH1I*)fList->FindObject(fEvents.Data());
3d4a1473 1473 // "nMCEvents"
7e2bf482 1474 TH1I* hPrimEventsTrVtx = (TH1I*)fList->FindObject(fPrimEvents.Data());
1475
1476 AliFMDAnaCalibSharingEfficiency* sharEff = new AliFMDAnaCalibSharingEfficiency();
1477
1478 for(Int_t det = 1; det<=3; det++) {
1479 Int_t maxRing = (det == 1 ? 0 : 1);
1480 for(Int_t ring = 0;ring<=maxRing;ring++) {
1481 Char_t ringChar = (ring == 0 ? 'I' : 'O');
1482 for(Int_t v=0; v< nVertexBins; v++) {
3d4a1473 1483 // Get histograms like hits_NoCuts_FMD%d%c_vtxbin%d_proj
1484 // - from AliFMDAnalysisTaskBackgroundCorrection.cxx
7e2bf482 1485 TH1F* hHits = (TH1F*)fList->FindObject(GetAnalysisName(kHits,det, ringChar, v));
3d4a1473 1486 // Get histograms like hits_FMD%d%c_vtxbin%d_proj
1487 // - from AliFMDAnalysisTaskBackgroundCorrection.cxx
7e2bf482 1488 TH1F* hHitsTrVtx = (TH1F*)fList->FindObject(GetAnalysisName(kHitsTrVtx,det, ringChar, v));
3d4a1473 1489 // Get histograms like "hMCHits_nocuts_FMD%d%c_vtxbin%d"
1490 // - from AliFMDAnalysisTaskSharing.cxx
7e2bf482 1491 TH1F* hMCHits = (TH1F*)fList->FindObject(GetPrimName(kHits,det, ringChar, v));
3d4a1473 1492 // Get histograms like "hMCHits_FMD%d%c_vtxbin%d"
1493 // - from AliFMDAnalysisTaskSharing.cxx
7e2bf482 1494 TH1F* hMCHitsTrVtx = (TH1F*)fList->FindObject(GetPrimName(kHitsTrVtx,det, ringChar, v));
1495
541c19ed 1496 TH1F* hCorrection = (TH1F*)hHits->Clone(Form("hCorrection_FMD%d%c_vtx%d",det, ringChar, v));
1497 TH1F* hCorrectionTrVtx = (TH1F*)hHitsTrVtx->Clone(Form("hCorrection_FMD%d%c_vtx%d",det, ringChar, v));
059c7c6b 1498 if(hEvents->GetBinContent(v+1))
1499 hCorrection->Scale(1./(Float_t)hEvents->GetBinContent(v+1));
1500 if(hPrimEvents->GetBinContent(v+1))
1501 hMCHits->Scale(1./(Float_t)hPrimEvents->GetBinContent(v+1));
1502 if(hEventsTrVtx->GetBinContent(v+1))
1503 hCorrectionTrVtx->Scale(1./(Float_t)hEventsTrVtx->GetBinContent(v+1));
1504 if(hPrimEventsTrVtx->GetBinContent(v+1))
1505 hMCHitsTrVtx->Scale(1./(Float_t)hPrimEventsTrVtx->GetBinContent(v+1));
7e2bf482 1506
7e2bf482 1507 hCorrection->Divide(hMCHits);
1508 hCorrectionTrVtx->Divide(hMCHitsTrVtx);
1509
059c7c6b 1510 //sharEff->SetSharingEff(det,ringChar,v,hCorrection);
7e2bf482 1511 sharEff->SetSharingEff(det,ringChar,v,hCorrection);
1512 sharEff->SetSharingEffTrVtx(det,ringChar,v,hCorrectionTrVtx);
059c7c6b 1513
7e2bf482 1514 // std::cout<<hHits<<" "<<hHitsTrVtx<<" "<<hPrim<<" "<<hPrimTrVtx<<std::endl;
1515
1516 }
1517 }
1518 }
1519
1520 if(store) {
1521 TFile fsharing(pars->GetPath(pars->GetSharingEffID()),"RECREATE");
1522 sharEff->Write(AliFMDAnaParameters::GetSharingEffID());
1523 fsharing.Close();
1524 }
1525
c5058a61 1526}
2dcdd4eb 1527//_____________________________________________________________________
1528void AliFMDDndeta::RebinHistogram(TH1F* hist, Int_t rebin) {
1529 //Clever rebinner
1530
1531 if(hist->GetNbinsX()%rebin) {
1532 std::cout<<"cannot rebin "<<hist->GetNbinsX()%rebin<< "is not zero"<<std::endl;
1533 return;
c5058a61 1534
2dcdd4eb 1535 }
1536 TH1F* cloneHist = (TH1F*)hist->Clone();
1537 cloneHist->Rebin(rebin);
1538
1539 Int_t nBinsNew = hist->GetNbinsX() / rebin;
1540
1541 for(Int_t i =1;i<=nBinsNew; i++) {
1542 Float_t bincontent = 0;
1543 Float_t sumofw = 0;
1544 Float_t sumformean = 0;
1545 Int_t nBins = 0;
1546 for(Int_t j=1; j<=rebin;j++) {
1547 if(hist->GetBinContent((i-1)*rebin + j) > 0) {
1548 bincontent = bincontent + hist->GetBinContent((i-1)*rebin + j);
1549 nBins++;
1550 Float_t weight = 1 / TMath::Power(hist->GetBinError((i-1)*rebin + j),2);
1551 sumofw = sumofw + weight;
1552 sumformean = sumformean + weight*hist->GetBinContent((i-1)*rebin + j);
1553
1554 }
1555 }
1556
1557 if(bincontent > 0 ) {
1558 cloneHist->SetBinContent(i,sumformean / sumofw);
1559 cloneHist->SetBinError(i,TMath::Sqrt(sumformean));//TMath::Sqrt(1/sumofw));
1560 }
1561 }
1562 hist->Rebin(rebin);
1563 for(Int_t i =1;i<=nBinsNew; i++) {
1564 hist->SetBinContent(i,cloneHist->GetBinContent(i));
1565 }
1566
1567
1568}
c5058a61 1569//_____________________________________________________________________
1570//
1571// EOF
1572//
5be76c19 1573// EOF