Coding violation corrected (G.Galazka)
[u/mrichter/AliRoot.git] / HBTAN / AliHBTWeightasCorrFctn.cxx
CommitLineData
94a709e1 1#include "AliHBTWeightasCorrFctn.h"
2#include <TH1.h>
265f04f7 3#include <Riostream.h>
94a709e1 4
5///////////////////////////////////////////////////////
6// //
7// AliHBTWeightasCorrFctn.h //
8// //
9// Class for calculating 3D Weightas correlation //
10// functions //
11// //
12///////////////////////////////////////////////////////
13
14ClassImp(AliHBTWeightasCorrFctn)
15
16 AliHBTWeightasCorrFctn::AliHBTWeightasCorrFctn(const char* name, const char* title):
17 AliHBTTwoPairFctn1D(name,title),
18
19 fNum(0x0),
20 fDen(0x0),
21 fRat(0x0)
22
23{
24//ctor
25}
26
27/******************************************************************/
28AliHBTWeightasCorrFctn::AliHBTWeightasCorrFctn(const char* name, const char* title, Int_t nbins, Float_t maxXval, Float_t minXval):
29 AliHBTTwoPairFctn1D(name,title,nbins,maxXval,minXval),
30
31
32 fNum(new TObjArray()),
33 fDen(new TObjArray()),
34 fRat(new TObjArray())
35{
36 SetParams(nbins,maxXval, minXval);
37}
38
39/******************************************************************/
40AliHBTWeightasCorrFctn::AliHBTWeightasCorrFctn(const AliHBTWeightasCorrFctn& in):
41 AliHBTTwoPairFctn1D(in),
42
43
44
45 fNum((in.fNum)?(TObjArray*)in.fNum->Clone():0x0),
46 fDen((in.fDen)?(TObjArray*)in.fDen->Clone():0x0),
47 fRat((in.fRat)?(TObjArray*)in.fRat->Clone():0x0)
48 {
49//ctor
50}
51
52/******************************************************************/
53
54AliHBTWeightasCorrFctn::~AliHBTWeightasCorrFctn()
55{
56 //dtor
57
58 delete fNum;
59 delete fDen;
60 delete fRat;
61
62}
63
64/******************************************************************/
65void AliHBTWeightasCorrFctn::Write()
66{
67//out
68 Int_t i;
69// Int_t n=GetNumberOfIntervals();
70 Double_t scale;
71
72 for(i=0;i<fNumberOfIntervals;i++){
73 TH1D *num = ((TH1D*)fNum->At(i));
74 TH1D *den = ((TH1D*)fDen->At(i));
75 TH1D &rat = *((TH1D*)fRat->At(i));
76 scale = Scale(num,den);
77 Info("Write():","Scale in interval %d = %lf",i,scale);
78 rat.Divide(num,den,scale);
79
80
81 num->Write();
82 den->Write();
83 rat.Write();
84 }
85
86}
87
88//-------------------------------------
89void AliHBTWeightasCorrFctn::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
90{
91 //Fills the numerator using pair from the same event
92 trackpair = CheckPair(trackpair);
93 if(partpair == 0x0) return;
94 if(trackpair == 0x0) return;
95
96 Double_t weight = 1.0;
97
98 int n = fNumberOfIntervals;
99 Double_t rplane=0.; //reaction plane angle - 2 B determined
100 Double_t phi=(trackpair->Particle1()->Phi()+trackpair->Particle2()->Phi())/2.-rplane; //deltaphi bo nie mam nic innego pod reka
101
102 phi=phi*360/(2*TMath::Pi());
103 Double_t q=GetValue(trackpair, partpair);
104 Int_t ntv;
105 ntv = (int)(phi*n/(360.));
106
107 TH1D *num = ((TH1D*)fNum->At(ntv));
108
109
110 if ( trackpair && partpair)
111 {
112 if ( ( trackpair->Particle1()->GetPdgCode() == partpair->Particle1()->GetPdgCode()) &&
113 ( trackpair->Particle2()->GetPdgCode() == partpair->Particle2()->GetPdgCode()) )
114 {
115 weight=partpair->GetWeight();
116// Info("ProcessSameEvent","weight=%lf",weight);
117 }
118 num->Fill(q,weight);
119 }
120}
121
122/****************************************************************/
123void AliHBTWeightasCorrFctn::Init()
124{
125 BuildHistos();
126}
127/****************************************************************/
128
129
130void AliHBTWeightasCorrFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
131{
265f04f7 132 //process particles from different events
94a709e1 133 Double_t rplane=0.; //reaction plane angle - 2 B determined
134 Double_t phi=(trackpair->Particle1()->Phi()+trackpair->Particle2()->Phi())/2.-rplane; //deltaphi bo nie mam nic innego pod reka
135 phi=phi*360/(2*TMath::Pi());
265f04f7 136 trackpair = CheckPair(trackpair);
137
94a709e1 138
139 int n = fNumberOfIntervals;
140
141 Int_t ntv;
142 ntv = (int)(phi*n/(360.));
143
144 TH1D &den = *((TH1D*)fDen->At(ntv));
145 if ( trackpair && partpair)
146 {
265f04f7 147 Double_t qout=GetValue(trackpair, partpair);
148 den.Fill(qout);
94a709e1 149 }
150}
151
152
153/******************************************************************/
154
155
265f04f7 156void AliHBTWeightasCorrFctn::SetParams(Int_t nbins, Float_t maxXval, Float_t minXval)
157{
158//set histogram parameters
94a709e1 159 fnbins=nbins;
160 fmaxXval= maxXval;
161 fminXval=minXval;
162}
163TH1* AliHBTWeightasCorrFctn::GetResult()
164{
265f04f7 165//just for compatibility
94a709e1 166 TH1D *den = ((TH1D*)fDen->UncheckedAt(1));
167 return den;
168 }
169
170
171
172
173ClassImp(AliHBTQOutWeightasCorrFctn)
174
175 AliHBTQOutWeightasCorrFctn::AliHBTQOutWeightasCorrFctn(const char* name, const char* title, Int_t nbins, Float_t maxXval, Float_t minXval):
176AliHBTWeightasCorrFctn(name,title,nbins,maxXval,minXval)
177
178{
179//ct0r
180}
181
182void AliHBTQOutWeightasCorrFctn::BuildHistos()
183{
265f04f7 184//buils histograms (allocates memory etc)
94a709e1 185 Int_t i;
186 int n=GetNumberOfIntervals();
187
188 int nbins=Getnbins();
189
190 double max = GetmaxXval();
191 double min = GetminXval();
192 char buff[10];
193
194
195
265f04f7 196 TH1D *num;
197 TH1D *den;
198 TH1D *rat;
94a709e1 199
265f04f7 200 TString nameNum = "WeightNumOut";
201 TString nameDen = "WeightDenOut";
202 TString nameRat = "WeightRatOut";
94a709e1 203
204 for(i=0;i<n;i++){
205
206 sprintf(buff,"%d",i);
207
208 nameNum +=TString(buff);
209
210 nameDen +=TString(buff);
211 nameRat +=TString(buff);
212
213
265f04f7 214 num = new TH1D(nameNum.Data(),nameNum.Data(),nbins,min,max);
215 den = new TH1D(nameDen.Data(),nameDen.Data(),nbins,min,max);
216 rat = new TH1D(nameRat.Data(),nameRat.Data(),nbins,min,max);
94a709e1 217
265f04f7 218 num->Sumw2();
219 den->Sumw2();
220 rat->Sumw2();
94a709e1 221
265f04f7 222 num->Reset();
223 den->Reset();
224 rat->Reset();
94a709e1 225
265f04f7 226 fNum->Add(num);
227 fDen->Add(den);
228 fRat->Add(rat);
94a709e1 229
265f04f7 230 nameNum = TString("WeightNumOut");
231 nameDen = TString("WeightDenOut");
232 nameRat = TString("WeightRatOut");
94a709e1 233
234 }
235
236
237 }
238
239ClassImp(AliHBTQSideWeightasCorrFctn)
240
241 AliHBTQSideWeightasCorrFctn::AliHBTQSideWeightasCorrFctn(const char* name, const char* title, Int_t nbins, Float_t maxXval, Float_t minXval):
242AliHBTWeightasCorrFctn(name,title,nbins,maxXval,minXval)
243
244{
245//ct0r
246}
247
248void AliHBTQSideWeightasCorrFctn::BuildHistos()
249{
265f04f7 250//builds histograms
94a709e1 251 Int_t i;
252 int n=GetNumberOfIntervals();
253 int nbins=Getnbins();
254
255 double max = GetmaxXval();
256 double min = GetminXval();
257 char buff[10];
258
259
260
265f04f7 261 TH1D *num;
262 TH1D *den;
263 TH1D *rat;
94a709e1 264
265f04f7 265
266 TString nameNum = "WeightNumSide";
267 TString nameDen = "WeightDenSide";
268 TString nameRat = "WeightRatSide";
94a709e1 269
270 for(i=0;i<n;i++){
271
272 sprintf(buff,"%d",i);
273
274 nameNum +=TString(buff);
275
276 nameDen +=TString(buff);
277 nameRat +=TString(buff);
278
279
265f04f7 280 num = new TH1D(nameNum.Data(),nameNum.Data(),nbins,min,max);
281 den = new TH1D(nameDen.Data(),nameDen.Data(),nbins,min,max);
282 rat = new TH1D(nameRat.Data(),nameRat.Data(),nbins,min,max);
94a709e1 283
265f04f7 284 num->Sumw2();
285 den->Sumw2();
286 rat->Sumw2();
94a709e1 287
265f04f7 288 num->Reset();
289 den->Reset();
290 rat->Reset();
94a709e1 291
265f04f7 292 fNum->Add(num);
293 fDen->Add(den);
294 fRat->Add(rat);
94a709e1 295
265f04f7 296 nameNum = TString("WeightNumSide");
297 nameDen = TString("WeightDenSide");
298 nameRat = TString("WeightRatSide");
94a709e1 299
300 }
301
302
303 }
304
305ClassImp(AliHBTQLongWeightasCorrFctn)
306
307 AliHBTQLongWeightasCorrFctn::AliHBTQLongWeightasCorrFctn(const char* name, const char* title, Int_t nbins, Float_t maxXval, Float_t minXval):
308AliHBTWeightasCorrFctn(name,title,nbins,maxXval,minXval)
309
310{
311//ct0r
312}
313
314void AliHBTQLongWeightasCorrFctn::BuildHistos()
315{
265f04f7 316//builds histograms
94a709e1 317 Int_t i;
318 int n=GetNumberOfIntervals();
319 int nbins=Getnbins();
320
321 double max = GetmaxXval();
322 double min = GetminXval();
323 char buff[10];
324
325
326
265f04f7 327 TH1D *num;
328 TH1D *den;
329 TH1D *rat;
94a709e1 330
265f04f7 331 TString nameNum = "WeightNumLong";
332 TString nameDen = "WeightDenLong";
333 TString nameRat = "WeightRatLong";
94a709e1 334
335 for(i=0;i<n;i++){
336
337 sprintf(buff,"%d",i);
338
339 nameNum +=TString(buff);
340
341 nameDen +=TString(buff);
342 nameRat +=TString(buff);
343
344
265f04f7 345 num = new TH1D(nameNum.Data(),nameNum.Data(),nbins,min,max);
346 den = new TH1D(nameDen.Data(),nameDen.Data(),nbins,min,max);
347 rat = new TH1D(nameRat.Data(),nameRat.Data(),nbins,min,max);
94a709e1 348
265f04f7 349 num->Sumw2();
350 den->Sumw2();
351 rat->Sumw2();
94a709e1 352
265f04f7 353 num->Reset();
354 den->Reset();
355 rat->Reset();
94a709e1 356
265f04f7 357 fNum->Add(num);
358 fDen->Add(den);
359 fRat->Add(rat);
94a709e1 360
265f04f7 361 nameNum = TString("WeightNumLong");
362 nameDen = TString("WeightDenLong");
363 nameRat = TString("WeightRatLong");
94a709e1 364
365 }
366
367
368 }
265f04f7 369
370/********************************************************************/
371ClassImp(AliHBTasWeightQOSLCorrFctn)
372 AliHBTasWeightQOSLCorrFctn::AliHBTasWeightQOSLCorrFctn(const char* name, const char* title):
373 AliHBTTwoPairFctn3D(name,title),
374
375 fNum(0x0),
376 fDen(0x0),
377 fRat(0x0)
378// fNumberOfIntervals(1)
379
380
381{
382//ctor
383}
384
385/******************************************************************/
386AliHBTasWeightQOSLCorrFctn::AliHBTasWeightQOSLCorrFctn(const char* name, const char* title, Int_t nXbins, Double_t maxXval, Double_t minXval, Int_t nYbins, Double_t maxYval, Double_t minYval,Int_t nZbins, Double_t maxZval, Double_t minZval):
387 AliHBTTwoPairFctn3D(name,title,nXbins,maxXval,minXval,nYbins,maxYval,minYval,nZbins,maxZval,minZval),
388
389
390 fNum(new TObjArray()),
391 fDen(new TObjArray()),
392 fRat(new TObjArray())
393// fNumberOfIntervals(1)
394
395{
396 SetParams(nXbins,maxXval,minXval,nYbins,maxYval,minYval,nZbins,maxZval,minZval);
397}
398
399/******************************************************************/
400AliHBTasWeightQOSLCorrFctn::AliHBTasWeightQOSLCorrFctn(const AliHBTasWeightQOSLCorrFctn& in):
401 AliHBTTwoPairFctn3D(in),
402
403
404
405 fNum((in.fNum)?(TObjArray*)in.fNum->Clone():0x0),
406 fDen((in.fDen)?(TObjArray*)in.fDen->Clone():0x0),
407 fRat((in.fRat)?(TObjArray*)in.fRat->Clone():0x0)
408 {
409//ctor
410}
411
412/******************************************************************/
413
414AliHBTasWeightQOSLCorrFctn::~AliHBTasWeightQOSLCorrFctn()
415{
416 //dtor
417
418 delete fNum;
419 delete fDen;
420 delete fRat;
421
422}
423
424/******************************************************************/
425void AliHBTasWeightQOSLCorrFctn::Write()
426{
427//out
428 Int_t i;
429// Int_t n=GetNumberOfIntervals();
430 // Double_t scale;
431
432 for(i=0;i<fNumberOfIntervals;i++){
433Info("Write()","%d",i);
434 TH3D *num = ((TH3D*)fNum->UncheckedAt(i));
435 TH3D *den = ((TH3D*)fDen->UncheckedAt(i));
436 TH3D &rat = *((TH3D*)fRat->At(i));
437 // scale = Scale(num,den);
438// Info("Write():","Scale in interval %d = %lf",i,scale);
439 rat.Divide(num,den);//,scale);
440
441
442 num->Write();
443 den->Write();
444 rat.Write();
445 }
446
447}
448
449//-------------------------------------
450void AliHBTasWeightQOSLCorrFctn::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
451{
452 //Fills the numerator using pair from the same event
453 trackpair = CheckPair(trackpair);
454 if(partpair == 0x0) return;
455 if(trackpair == 0x0) return;
456
457 Double_t weight = 1.0;
458
459 int n = fNumberOfIntervals;
460 Double_t rplane=0.; //reaction plane angle - 2 B determined
461 Double_t phi=(trackpair->Particle1()->Phi()+trackpair->Particle2()->Phi())/2.-rplane; //deltaphi bo nie mam nic innego pod reka
462
463 phi=phi*360/(2*TMath::Pi());
464 Double_t out = TMath::Abs(trackpair->GetQOutLCMS());
465 Double_t side = TMath::Abs(trackpair->GetQSideLCMS());
466 Double_t lon = TMath::Abs(trackpair->GetQLongLCMS());
467
468
469 Int_t ntv;
470 ntv = (int)(phi*n/(360.));
471
472 TH3D *num = ((TH3D*)fNum->At(ntv));
473
474
475 if ( trackpair && partpair)
476 {
477 if ( ( trackpair->Particle1()->GetPdgCode() == partpair->Particle1()->GetPdgCode()) &&
478 ( trackpair->Particle2()->GetPdgCode() == partpair->Particle2()->GetPdgCode()) )
479 {
480 weight=partpair->GetWeight();
481// Info("ProcessSameEvent","weight=%lf",weight);
482 }
483 num->Fill(out,side,lon,weight);
484 }
485}
486
487/****************************************************************/
488void AliHBTasWeightQOSLCorrFctn::Init()
489{
490 BuildHistos();
491}
492/****************************************************************/
493
494
495void AliHBTasWeightQOSLCorrFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
496{
497//fill denominator
498 Double_t rplane=0.; //reaction plane angle - 2 B determined
499 Double_t phi=(trackpair->Particle1()->Phi()+trackpair->Particle2()->Phi())/2.-rplane; //deltaphi bo nie mam nic innego pod reka
500 phi=phi*360/(2*TMath::Pi());
501
502 Double_t out = TMath::Abs(trackpair->GetQOutLCMS());
503 Double_t side = TMath::Abs(trackpair->GetQSideLCMS());
504 Double_t lon = TMath::Abs(trackpair->GetQLongLCMS());
505
506 int n = fNumberOfIntervals;
507
508 Int_t ntv;
509 ntv = (int)(phi*n/(360.));
510
511 TH3D &den = *((TH3D*)fDen->At(ntv));
512 if ( trackpair && partpair)
513 {
514 den.Fill(out,side,lon);
515 }
516}
517
518
519/******************************************************************/
520
521
522void AliHBTasWeightQOSLCorrFctn::SetParams(Int_t nXbins, Float_t maxXval, Float_t minXval,Int_t nYbins, Float_t maxYval, Float_t minYval,Int_t nZbins, Float_t maxZval, Float_t minZval)
523{
524//set histogram parametera
525 fnXbins=nXbins;
526 fmaxXval= maxXval;
527 fminXval=minXval;
528 fnYbins=nYbins;
529 fmaxYval= maxYval;
530 fminYval=minYval;
531
532 fnZbins=nZbins;
533 fmaxZval= maxZval;
534 fminZval=minZval;
535
536
537}
538TH1* AliHBTasWeightQOSLCorrFctn::GetResult()
539{
540 //just for compatibility
541 TH1D *den = ((TH1D*)fDen->UncheckedAt(1));
542 return den;
543 }
544
545void AliHBTasWeightQOSLCorrFctn::BuildHistos()
546{
547 //builds histograms
548 Int_t i;
549 int n=GetNumberOfIntervals();
550 int nXbins=GetnXbins();
551 double maxX = GetmaxXval();
552 double minX = GetminXval();
553
554 int nYbins=GetnYbins();
555 double maxY = GetmaxXval();
556 double minY = GetminXval();
557 int nZbins=GetnZbins();
558 double maxZ = GetmaxZval();
559 double minZ = GetminZval();
560 char buff[10];
561
562 TH3D *num;
563 TH3D *den;
564 TH3D *rat;
565
566 TString nameNum = "OSLWeightNum";
567 TString nameDen = "OSLWeightDen";
568 TString nameRat = "OSLWeightRat";
569
570 for(i=0;i<n;i++){
571
572 sprintf(buff,"%d",i);
573
574 nameNum +=TString(buff);
575
576 nameDen +=TString(buff);
577 nameRat +=TString(buff);
578
579
580 num = new TH3D(nameNum.Data(),nameNum.Data(),nXbins,minX,maxX,nYbins,minY,maxY,nZbins,minZ,maxZ);
581 den = new TH3D(nameDen.Data(),nameDen.Data(),nXbins,minX,maxX,nYbins,minY,maxY,nZbins,minZ,maxZ);
582 rat = new TH3D(nameRat.Data(),nameRat.Data(),nXbins,minX,maxX,nYbins,minY,maxY,nZbins,minZ,maxZ);
583
584 num->Sumw2();
585 den->Sumw2();
586 rat->Sumw2();
587
588 num->Reset();
589 den->Reset();
590 rat->Reset();
591
592 fNum->Add(num);
593 fDen->Add(den);
594 fRat->Add(rat);
595
596 nameNum = TString("OSLWeightNumLong");
597 nameDen = TString("OSLWeightDenLong");
598 nameRat = TString("OSLWeightRatLong");
599
600 }
601}
602Double_t AliHBTasWeightQOSLCorrFctn::Scale(TH3D *num, TH3D *den)
603{
604//scale histograms method, the default one didn't like me so I've done this one :)
605if (AliVAODParticle::GetDebug()>0) Info("Scale","Enetered Scale()");
606 if(!num)
607 {
608 Error("Scale","No numerator");
609 return 0.0;
610 }
611 if(!den)
612 {
613 Error("Scale","No denominator");
614 return 0.0;
615 }
616
617 if( (fNBinsToScaleX < 1) || (fNBinsToScaleY < 1) || (fNBinsToScaleZ < 1))
618 {
619 Error("Scale","Number of bins for scaling is smaller thnan 1");
620 return 0.0;
621 }
622 UInt_t nbinsX = num->GetNbinsX();
623 if (fNBinsToScaleX > nbinsX)
624 {
625 Error("Scale","Number of X bins for scaling is bigger thnan number of bins in histograms");
626 return 0.0;
627 }
628
629 UInt_t nbinsY = num->GetNbinsX();
630 if (fNBinsToScaleY > nbinsY)
631 {
632 Error("Scale","Number of Y bins for scaling is bigger thnan number of bins in histograms");
633 return 0.0;
634 }
635
636 UInt_t nbinsZ = num->GetNbinsZ();
637 if (fNBinsToScaleZ > nbinsZ)
638 {
639 Error("Scale","Number of Z bins for scaling is bigger thnan number of bins in histograms");
640 return 0.0;
641 }
642
643 if (AliVAODParticle::GetDebug()>0) Info("Scale","No errors detected");
644 Int_t offsetX = nbinsX - fNBinsToScaleX - 1; //bin that we start loop over bins in axis X
645 Int_t offsetY = nbinsY - fNBinsToScaleY - 1; //bin that we start loop over bins in axis Y
646 Int_t offsetZ = nbinsZ - fNBinsToScaleZ - 1; //bin that we start loop over bins in axis Z
647 Double_t densum = 0.0;
648 Double_t numsum = 0.0;
649
650 for (UInt_t k = offsetZ; k<nbinsZ; k++)
651 for (UInt_t j = offsetY; j<nbinsY; j++)
652 for (UInt_t i = offsetX; i<nbinsX; i++)
653 {
654 if ( num->GetBinContent(i,j,k) > 0.0 )
655 {
656 densum += den->GetBinContent(i,j,k);
657 numsum += num->GetBinContent(i,j,k);
658 }
659 }
660 if(AliVAODParticle::GetDebug() > 0)
661 Info("Scale","numsum=%f densum=%f fNBinsToScaleX=%d fNBinsToScaleY=%d fNBinsToScaleZ=%d",
662 numsum,densum,fNBinsToScaleX,fNBinsToScaleY,fNBinsToScaleZ);
663 if (numsum == 0) return 0.0;
664 Double_t ret = densum/numsum;
665 if(AliVAODParticle::GetDebug() > 0) Info("Scale","returning %f",ret);
666 return ret;
667}
668