Coding violation corrected (G.Galazka)
[u/mrichter/AliRoot.git] / HBTAN / AliHBTWeightasCorrFctn.cxx
1 #include "AliHBTWeightasCorrFctn.h"
2 #include <TH1.h>
3 #include <Riostream.h>
4
5 ///////////////////////////////////////////////////////
6 //                                                   //
7 // AliHBTWeightasCorrFctn.h                             //
8 //                                                   //
9 // Class for calculating 3D Weightas correlation        //
10 // functions                                         //
11 //                                                   //
12 ///////////////////////////////////////////////////////
13
14 ClassImp(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 /******************************************************************/
28 AliHBTWeightasCorrFctn::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 /******************************************************************/
40 AliHBTWeightasCorrFctn::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
54 AliHBTWeightasCorrFctn::~AliHBTWeightasCorrFctn()
55 {
56  //dtor
57
58      delete fNum;
59      delete fDen;
60      delete fRat;
61      
62 }
63
64 /******************************************************************/
65 void 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 //-------------------------------------
89 void 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 /****************************************************************/
123 void AliHBTWeightasCorrFctn::Init()
124 {
125      BuildHistos();
126 }
127 /****************************************************************/
128
129
130 void AliHBTWeightasCorrFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
131 {
132      //process particles from different events 
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());
136      trackpair = CheckPair(trackpair);
137      
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      {
147      Double_t qout=GetValue(trackpair, partpair);       
148              den.Fill(qout);    
149      }
150 }
151
152
153 /******************************************************************/
154
155
156 void AliHBTWeightasCorrFctn::SetParams(Int_t nbins, Float_t maxXval, Float_t minXval)
157 {
158 //set histogram parameters
159      fnbins=nbins;
160      fmaxXval= maxXval;
161      fminXval=minXval;
162 }
163 TH1* AliHBTWeightasCorrFctn::GetResult()
164 {
165 //just for compatibility       
166      TH1D *den = ((TH1D*)fDen->UncheckedAt(1));
167      return den;
168  }
169
170
171
172
173 ClassImp(AliHBTQOutWeightasCorrFctn)
174      
175      AliHBTQOutWeightasCorrFctn::AliHBTQOutWeightasCorrFctn(const char* name, const char* title, Int_t nbins, Float_t maxXval, Float_t minXval):
176 AliHBTWeightasCorrFctn(name,title,nbins,maxXval,minXval)
177
178 {
179 //ct0r
180 }
181
182 void AliHBTQOutWeightasCorrFctn::BuildHistos()
183 {
184 //buils histograms (allocates memory etc)
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      
196      TH1D *num;
197      TH1D *den;
198      TH1D *rat;
199      
200      TString nameNum = "WeightNumOut";
201      TString nameDen = "WeightDenOut";
202      TString nameRat = "WeightRatOut";
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           
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);
217           
218           num->Sumw2();
219           den->Sumw2();
220           rat->Sumw2();
221           
222           num->Reset();
223           den->Reset();
224           rat->Reset();
225           
226           fNum->Add(num);
227           fDen->Add(den);
228           fRat->Add(rat);
229           
230           nameNum = TString("WeightNumOut");
231           nameDen = TString("WeightDenOut");
232           nameRat = TString("WeightRatOut");
233           
234      }
235
236      
237  }
238
239 ClassImp(AliHBTQSideWeightasCorrFctn)
240      
241      AliHBTQSideWeightasCorrFctn::AliHBTQSideWeightasCorrFctn(const char* name, const char* title, Int_t nbins, Float_t maxXval, Float_t minXval):
242 AliHBTWeightasCorrFctn(name,title,nbins,maxXval,minXval)
243
244 {
245 //ct0r
246 }
247
248 void AliHBTQSideWeightasCorrFctn::BuildHistos()
249 {
250 //builds histograms
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      
261      TH1D *num;
262      TH1D *den;
263      TH1D *rat;
264      
265      
266      TString nameNum = "WeightNumSide";
267      TString nameDen = "WeightDenSide";
268      TString nameRat = "WeightRatSide";
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           
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);
283           
284           num->Sumw2();
285           den->Sumw2();
286           rat->Sumw2();
287           
288           num->Reset();
289           den->Reset();
290           rat->Reset();
291           
292           fNum->Add(num);
293           fDen->Add(den);
294           fRat->Add(rat);
295           
296           nameNum = TString("WeightNumSide");
297           nameDen = TString("WeightDenSide");
298           nameRat = TString("WeightRatSide");
299           
300      }
301
302      
303  }
304
305 ClassImp(AliHBTQLongWeightasCorrFctn)
306      
307      AliHBTQLongWeightasCorrFctn::AliHBTQLongWeightasCorrFctn(const char* name, const char* title, Int_t nbins, Float_t maxXval, Float_t minXval):
308 AliHBTWeightasCorrFctn(name,title,nbins,maxXval,minXval)
309
310 {
311 //ct0r
312 }
313
314 void AliHBTQLongWeightasCorrFctn::BuildHistos()
315 {
316 //builds histograms
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      
327      TH1D *num;
328      TH1D *den;
329      TH1D *rat;
330      
331      TString nameNum = "WeightNumLong";
332      TString nameDen = "WeightDenLong";
333      TString nameRat = "WeightRatLong";
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           
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);
348           
349           num->Sumw2();
350           den->Sumw2();
351           rat->Sumw2();
352           
353           num->Reset();
354           den->Reset();
355           rat->Reset();
356           
357           fNum->Add(num);
358           fDen->Add(den);
359           fRat->Add(rat);
360           
361           nameNum = TString("WeightNumLong");
362           nameDen = TString("WeightDenLong");
363           nameRat = TString("WeightRatLong");
364           
365      }
366
367      
368  }
369
370 /********************************************************************/
371 ClassImp(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 /******************************************************************/
386 AliHBTasWeightQOSLCorrFctn::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 /******************************************************************/
400 AliHBTasWeightQOSLCorrFctn::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
414 AliHBTasWeightQOSLCorrFctn::~AliHBTasWeightQOSLCorrFctn()
415 {
416  //dtor
417
418      delete fNum;
419      delete fDen;
420      delete fRat;
421      
422 }
423
424 /******************************************************************/
425 void 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++){
433 Info("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 //-------------------------------------
450 void 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 /****************************************************************/
488 void AliHBTasWeightQOSLCorrFctn::Init()
489 {
490      BuildHistos();
491 }
492 /****************************************************************/
493
494
495 void 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
522 void 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 }
538 TH1* AliHBTasWeightQOSLCorrFctn::GetResult()
539 {
540  //just for compatibility      
541      TH1D *den = ((TH1D*)fDen->UncheckedAt(1));
542      return den;
543  }
544
545 void 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 }
602 Double_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 :)
605 if (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