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