]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - EVGEN/AliGenHBTosl.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
fDebug is Int_t instead of Bool_t
[u/mrichter/AliRoot.git] / EVGEN / AliGenHBTosl.cxx
... / ...
CommitLineData
1#include "AliGenHBTosl.h"
2//__________________________________________________________
3/////////////////////////////////////////////////////////////
4// //
5// class AliGenHBTosl //
6// //
7// Genarator simulating particle correlations //
8// //
9// The main idea of the generator is to produce particles //
10// according to some distribution of two particle //
11// property. In HBT they are qout,qsie and qlong. //
12// In order to be able to generate signal that produces //
13// given two particle correlation background must be //
14// known before in order to produce the shape of signal //
15// to randomize given distribution from. //
16// //
17// The generator works as follows: //
18// 1. Coarse Background (fQCoarseBackground) is generated //
19// ade from the particles //
20// given by the external generator (variable //
21// fGenerator) by the mixing technique. //
22// 2. Coarse signal is prduced by multiplying Coarse //
23// background by a required function //
24// See method FillCoarseSignal //
25// 3. Signal is randomized out of the coarse signal //
26// histogram (two particle property). First particle //
27// is taken from the external generator, and the //
28// second one is CALCULATED on the basis of the first //
29// one and the two particle property (qout,qside,qlong)//
30// Background is made by the mixing out of the //
31// genereted signal events. //
32// This step is cotinued up to the moment signal //
33// histogram has enough statistics (data member //
34// fMinFill) //
35// See method StartSignalPass1() //
36// 4. chi is calculated for each bin (chiarray variqable) //
37// (not the chi2 because sign is important) //
38// Two particle prioperty //
39// (qout,qside,qlong) is chosen at the points that //
40// chi is the smallest. First particle is taken from //
41// the the external generator (fGenerator) and second's /
42// momenta are caclulated out of their momenta and //
43// (qout,qside,qlong). Background is updated //
44// continuesely for all the events. This step is //
45// continued until stability conditions are fullfiled //
46// or maximum number of iteration is reached. //
47// 5. The same as step 4 but events are stored. //
48// //
49////////////////////////////////////////////////////////////
50
51#include <Riostream.h>
52#include <TCanvas.h>
53
54
55#include <TH3D.h>
56#include <TList.h>
57#include <TPDGCode.h>
58#include <TParticle.h>
59#include <AliStack.h>
60#include <TMath.h>
61#include <TVector3.h>
62#include <TStopwatch.h>
63#include <TFile.h>
64
65#include "AliGenCocktailAfterBurner.h"
66#include "AliGeVSimParticle.h"
67#include "AliGenGeVSim.h"
68#include "AliGenHIJINGpara.h"
69
70/***********************************************************/
71ClassImp(AliGenHBTosl)
72
73AliGenHBTosl::AliGenHBTosl():
74 AliGenerator(),
75 fQCoarseBackground(0x0),
76 fQCoarseSignal(0x0),
77 fQSignal(0x0),
78 fQBackground(0x0),
79 fQSecondSignal(0x0),
80 fQSecondBackground(0x0),
81 fQRange(0.15),
82 fQNBins(60),
83 fGenerator(0x0),
84 fBufferSize(5),
85 fNBinsToScale(Int_t(fQNBins*0.2)),
86 fDebug(0),
87 fMaxIterations(300),
88 fMaxChiSquereChange(0.01),
89 fMaxChiSquerePerNDF(1.5),
90 fQRadius(8.0),
91 fPID(kPiPlus),
92 fSamplePhiMin(-0.01),
93 fSamplePhiMax(TMath::TwoPi()+0.01),
94 fSignalRegion(0.0),
95 fMinFill(300)
96{
97//default constructor
98}
99/***********************************************************/
100
101AliGenHBTosl::AliGenHBTosl(Int_t n, Int_t pid):
102 AliGenerator(n),
103 fQCoarseBackground(0x0),
104 fQCoarseSignal(0x0),
105 fQSignal(0x0),
106 fQBackground(0x0),
107 fQSecondSignal(0x0),
108 fQSecondBackground(0x0),
109 fQRange(0.15),
110 fQNBins(60),
111 fGenerator(0x0),
112 fBufferSize(5),
113 fNBinsToScale(Int_t(fQNBins*0.2)),
114 fDebug(0),
115 fSignalShapeCreated(kFALSE),
116 fMaxIterations(300),
117 fMaxChiSquereChange(0.01),
118 fMaxChiSquerePerNDF(1.5),
119 fQRadius(8.0),
120 fPID(pid),
121 fSamplePhiMin(-0.01),
122 fSamplePhiMax(TMath::TwoPi()+0.01),
123 fSignalRegion(0.0),
124 fMinFill(300)
125{
126//default constructor
127}
128/***********************************************************/
129
130AliGenHBTosl::~AliGenHBTosl()
131{
132//destructor
133 delete fQCoarseSignal;
134 delete fQCoarseBackground;
135 delete fQSignal;
136 delete fQBackground;
137 delete fGenerator;
138 delete fQSecondSignal;
139 delete fQSecondBackground;
140
141}
142
143void AliGenHBTosl::Init()
144{
145 //Initializes generator
146 if (fGenerator == 0x0)
147 {
148
149 AliGenHIJINGpara* bkggen = new AliGenHIJINGpara(fNpart*4);
150 fGenerator = bkggen;
151
152/*
153 AliGenGeVSim * gevsim = new AliGenGeVSim(0.0);
154 AliGeVSimParticle* kplus = new AliGeVSimParticle(fPID,1,fNpart, 0.17, 0.9);
155 gevsim->AddParticleType(kplus);
156
157 fGenerator = gevsim;
158*/
159
160/*
161
162 AliMevSimConfig *c = new AliMevSimConfig(1);
163 c->SetRectPlane(1); // reaction plane control, model 4
164 c->SetGrid(80,80);
165
166 AliGenMevSim *mevsim = new AliGenMevSim(c);
167 mevsim->SetPtRange(0.001, 3);
168 mevsim->SetMomentumRange(0.1, 3);
169 mevsim->SetTrackingFlag(0);
170 mevsim->SetOrigin(0.0, 0.0, 0.0);
171 mevsim->SetSigma(0.0, 0.0, 0.0);
172 AliMevSimParticle *kplus = new AliMevSimParticle(kKPlus, fNpart, 0, 0.25, 0.0, 2, 0.15, 0.0, 0.0 );
173 mevsim->AddParticleType(kplus);
174 fGenerator = mevsim;
175*/
176
177 fGenerator->SetOrigin(fOrigin[0],fOrigin[1],fOrigin[2]);
178 static const Double_t degtoradcf = 180./TMath::Pi();
179 fGenerator->SetMomentumRange(fPtMin,fPtMax);
180 fGenerator->SetPhiRange(degtoradcf*fPhiMin,degtoradcf*fPhiMax);
181 fGenerator->SetYRange(fYMin,fYMax);
182 fGenerator->SetThetaRange(degtoradcf*fThetaMin,degtoradcf*fThetaMax);
183 fGenerator->Init();
184
185 }
186
187 fQCoarseBackground = new TH3D("fQCoarseBackground","",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
188 fQCoarseSignal = new TH3D("fQCoarseSignal","fQCoarseSignal",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
189 fQSignal = new TH3D("fQSignal1","fQSignal",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
190 fQBackground = new TH3D("fQBackground1","fQBackground",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
191
192 fQSecondSignal = new TH3D("fQSignal2","fQSignal",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
193 fQSecondBackground = new TH3D("fQBackground2","fQBackground",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
194
195 fQCoarseBackground->Sumw2();
196 fQCoarseSignal->Sumw2();
197 fQSignal->Sumw2();
198 fQBackground->Sumw2();
199 fQSecondSignal->Sumw2();
200 fQSecondBackground->Sumw2();
201}
202/***********************************************************/
203
204void AliGenHBTosl::Generate()
205{
206 //the main method
207 Info("Generate","Attempts to generate %d particles.",fNpart);
208 if (fStackBuffer == 0x0) fStackBuffer = new TList();
209 //Here is initialization level
210 if (fSignalShapeCreated == kFALSE)
211 {
212 TH3D *hs = 0x0, *hb = 0x0;
213 TFile* file;
214
215 file = TFile::Open("QTSignal.root");
216 if (file)
217 {
218 hs = (TH3D*)file->Get("fQSignal1");
219 if (hs) hs->SetDirectory(0x0);
220 }
221 delete file;
222
223 file = TFile::Open("QTBackground.root");
224 if (file)
225 {
226 hb = (TH3D*)file->Get("fQBackground1");
227 if (hb) hb->SetDirectory(0x0);
228 }
229 delete file;
230
231 if (hs && hb)
232 {
233 Info("Generate","**********************************");
234 Info("Generate","Found starting histograms in files");
235 Info("Generate","**********************************");
236 delete fQSignal;
237 delete fQBackground;
238 fQSignal = hs;
239 fQBackground = hb;
240 }
241 else
242 {
243 TH3D *cs = 0x0, *cb = 0x0;
244 file = TFile::Open("QTCoarseBackground.root");
245 if (file)
246 {
247 cb = (TH3D*)file->Get("fQCoarseBackground");
248 if (cb) cb->SetDirectory(0x0);
249 }
250 delete file;
251
252 file = TFile::Open("QTCoarseSignal.root");
253 if (file)
254 {
255 cs = (TH3D*)file->Get("fQCoarseSignal");
256 if (cs) cs->SetDirectory(0x0);
257 }
258 delete file;
259
260 if (cs && cb)
261 {
262
263 Info("Generate","Got Coarse signal and bkg from files");
264 delete fQCoarseBackground;
265 delete fQCoarseSignal;
266 fQCoarseSignal = cs;
267 fQCoarseBackground = cb;
268 }
269 else
270 {
271 if (cb)
272 {
273 Info("Generate","Got Coarse bkg from file");
274 delete fQCoarseBackground;
275 fQCoarseBackground = cb;
276 }
277 else
278 {
279 FillCoarse(); //create coarse background - just to know the spectrum
280 }
281 FillCoarseSignal();//create first coarse signal by brutal multplication coarse background and required function shape
282 }
283
284 StartSignal(); //Initilizes the stack that is used for generation
285 }
286 fSignalShapeCreated = kTRUE;
287 }
288
289 AliStack* stack = RotateStack();
290
291 AliStack* genstack = fGenerator->GetStack();
292 if (genstack == 0x0)
293 {
294 genstack = new AliStack(fNpart);
295 fGenerator->SetStack(genstack);
296 }
297 else
298 {
299 genstack->Reset();
300 }
301
302 fGenerator->Generate();
303 Int_t j = 0, ntr = 0;
304 if ( genstack->GetNtrack() < fNpart/2)
305 {
306 Warning("Generate","************************************************************");
307 Warning("Generate","Generator generated (%d) less particles then expected (%d).",
308 stack->GetNtrack(),fNpart/2);
309 Warning("Generate","************************************************************");
310 }
311
312 TH3D* work = new TH3D("work","work",fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange);
313
314
315 Double_t*** chiarray = new Double_t** [fQNBins+1];
316 Double_t*** sigarray = new Double_t** [fQNBins+1];
317
318 for (Int_t i = 1; i<=fQNBins; i++)
319 {
320 chiarray[i] = new Double_t* [fQNBins+1];
321 sigarray[i] = new Double_t* [fQNBins+1];
322
323 for (Int_t k = 1; k<=fQNBins; k++)
324 {
325 chiarray[i][k] = new Double_t [fQNBins+1];
326 sigarray[i][k] = new Double_t [fQNBins+1];
327 }
328 }
329
330
331 Double_t scale = Scale(fQSignal,fQBackground);
332 work->Divide(fQSignal,fQBackground,scale);
333
334 Double_t binwdh = work->GetBinWidth(1)/2.;
335
336 for (Int_t k = 1; k<=fQNBins; k++)
337 {
338 Double_t z = work->GetZaxis()->GetBinCenter(k);
339 for (Int_t j = 1; j<=fQNBins; j++)
340 {
341 Double_t y = work->GetYaxis()->GetBinCenter(j);
342 for (Int_t i = 1; i<=fQNBins; i++)
343 {
344 sigarray[i][j][k] = fQSignal->GetBinContent(i,j,k);//store current value of signal histogram
345 Double_t x = work->GetXaxis()->GetBinCenter(i);//get center value of a bin (qinv)
346 Double_t v = GetQOutQSideQLongCorrTheorValue(x,y,z);//get expected value of CF in that qinv
347 Double_t diff = v - work->GetBinContent(i,j,k);//store difference betweeon current value, and desired value
348 chiarray[i][j][k] = diff; // no-x x is a weight to get good distribution
349 }
350 }
351 }
352
353 printf("\n");
354 Int_t middlebin = fQNBins/2;
355 for (Int_t k = middlebin-5; k < middlebin+5; k++)
356 {
357 printf("% 6.5f ",work->GetBinContent(1,1,k));
358 }
359 printf("\n");
360
361 for (Int_t k = middlebin-5; k < middlebin+5; k++)
362 {
363 printf("% 6.5f ",chiarray[1][1][k]);
364 }
365 printf("\n");
366
367 TParticle particle(fPID,0,-1,-1,-1,-1,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);
368 TParticle* second = &particle;
369
370 Bool_t shortloop = kTRUE;
371 Int_t sc = 0;//security check against infinite loop
372
373 while ( (ntr+1) < fNpart)
374 {
375 Int_t xmax = 1;
376 Int_t ymax = 1;
377 Int_t zmax = 1;
378 Double_t qout;
379 Double_t qside;
380 Double_t qlong;
381
382
383 Int_t loopmin;
384 Int_t loopmax;
385
386 if (shortloop)
387 {
388 shortloop = kFALSE;
389 loopmax = fQNBins;
390 loopmin = 1;
391 }
392 else
393 {
394 shortloop = kTRUE;
395 loopmax = fQNBins/2+fQNBins/4;
396 loopmin = fQNBins/2-fQNBins/4;
397 }
398
399
400 for (Int_t k = loopmin; k <=loopmax; k++ )
401 {
402 qlong = work->GetZaxis()->GetBinCenter(k);
403 for (Int_t j = loopmin; j<=loopmax; j++)
404 {
405 qside = work->GetYaxis()->GetBinCenter(j);
406 for (Int_t i = loopmin; i<=loopmax; i++)
407 {
408 qout = work->GetXaxis()->GetBinCenter(i);
409 if (chiarray[xmax][ymax][zmax] < chiarray[i][j][k])
410 {
411 xmax = i;
412 ymax = j;
413 zmax = k;
414 }
415
416// Double_t qdist = TMath::Sqrt(qout*qout + qside*qside + qlong*qlong);
417
418// Double_t fact = chiarray[i][j][k];//chiarray is chi2
419// if (fact > work->GetBinError(i,j,k))//if differece between what we want and
420// { //what we have is bigger than stat. error
421// xmax = i; //we force to fill that bin
422// ymax = j;
423// zmax = k;
424// break;
425// }
426 }
427 }
428 }
429 qlong = work->GetZaxis()->GetBinCenter(zmax);
430 qside = work->GetYaxis()->GetBinCenter(ymax);
431 qout = work->GetXaxis()->GetBinCenter(xmax);
432
433 Info("Generate","Fill bin (%d,%d,%d)",xmax,ymax,zmax);
434
435 qout = gRandom->Uniform(qout-binwdh,qout+binwdh);
436 qside = gRandom->Uniform(qside-binwdh,qside+binwdh);
437 qlong = gRandom->Uniform(qlong-binwdh,qlong+binwdh);
438
439 TParticle* first = 0;
440 while (j < genstack->GetNtrack())
441 {
442 TParticle* tmpp = genstack->Particle(j++);
443 if (tmpp->GetPdgCode() == fPID)
444 {
445 if (CheckParticle(tmpp,0x0,stack) == kFALSE)
446 {
447 first = tmpp;
448 break;
449 }
450 }
451 }
452
453 if (first == 0x0)
454 {
455 if ( fDebug > 2 ) Info("StartSignal","No more particles of that type");
456 break;
457 }
458
459 //Here put the check if this particle do not fall into signal region with other partticle
460
461 Int_t retval = GetThreeD(first,second,qout,qside,qlong);
462 if (retval)
463 {
464 //Info("StartSignal","Can not find momenta for this OSL and particle");
465 continue;
466 }
467 //in case this particle is falling into signal area with another
468 //particle we take a another pair
469 //it can intruduce artificial correlations
470 Bool_t checkresult = CheckParticle(second,first,stack);
471 if ( checkresult && (sc < 10) )
472 {
473 sc++;
474 continue;
475 }
476 sc = 0;
477
478 //Put on output stack
479 SetTrack(first,ntr);
480 SetTrack(second,ntr);
481
482 //Put on internal stack
483 Int_t etmp;
484 SetTrack(first,etmp,stack);
485 SetTrack(second,etmp,stack);
486
487 Double_t y = GetQOutQSideQLongCorrTheorValue(qout,qside,qlong);
488
489 sigarray[xmax][ymax][zmax] ++;
490 chiarray[xmax][ymax][zmax] = scale*sigarray[xmax][ymax][zmax]/fQBackground->GetBinContent(xmax,ymax,zmax);
491 chiarray[xmax][ymax][zmax] = (y - chiarray[xmax][ymax][zmax]);
492
493 }
494
495 Mix(fStackBuffer,fQBackground,fQSecondSignal); //upgrate background
496 Mix(stack,fQSignal,fQSecondBackground); //upgrate signal
497
498 for (Int_t i = 1; i<=fQNBins; i++)
499 {
500 for (Int_t k = 1; k<=fQNBins; k++)
501 {
502 delete [] chiarray[i][k];
503 delete [] sigarray[i][k];
504 }
505 delete [] chiarray[i];
506 delete [] sigarray[i];
507 }
508 delete [] chiarray;
509 delete [] sigarray;
510}
511/***********************************************************/
512
513void AliGenHBTosl::GetOneD(TParticle* first, TParticle* second,Double_t qinv)
514{
515//deprecated method that caclulates momenta of the second particle
516// out of qinv and the first particle
517//first particle is rotated that only X is non-zero
518
519
520 Double_t m = first->GetMass();
521 Double_t msqrd = m*m;
522 Double_t fourmassSquered = 4.*msqrd;
523
524 //Condition that R must fullfill to be possible to have qinv less smaller then randomized
525// Double_t rRange = qinv*TMath::Sqrt(qinv*qinv + fourmassSquered)/fourmassSquered;
526// Double_t r = gRandom->Uniform(rRange);
527
528 Double_t r = gRandom->Uniform(qinv);
529 Double_t phi = gRandom->Uniform(TMath::TwoPi());
530
531 Double_t firstPx = first->P();//first particle is rotated that only X is non-zero thus P==Px
532 Double_t px = 2.*msqrd*firstPx + firstPx*qinv*qinv;
533 Double_t temp = qinv*qinv*qinv*qinv + fourmassSquered * (qinv*qinv - r*r );
534 if (temp < 0.0)
535 {
536 Error("GetOneD","temp is less then 0: %f",temp);
537 return;
538 }
539 temp = temp*(msqrd+firstPx*firstPx);
540
541 px = (px - TMath::Sqrt(temp))/(2.*msqrd);
542
543 Double_t py = r*TMath::Sin(phi);
544 Double_t pz = r*TMath::Cos(phi);
545
546 TVector3 firstpvector(first->Px(),first->Py(),first->Pz());
547 TVector3 vector(px,py,pz);
548 Rotate(firstpvector,vector);
549
550 Double_t e = TMath::Sqrt(msqrd + vector.X()*vector.X() + vector.Y()*vector.Y() + vector.Z()*vector.Z());
551 second->SetMomentum(vector.X(),vector.Y(),vector.Z(),e);
552// TParticle* f = new TParticle(first->GetPdgCode(),0,-1,-1,-1,-1, firstPx,0,0,e=TMath::Sqrt(msqrd+firstPx*firstPx),0.0,0.0,0.0,0.0);
553// TParticle(pdg, is, parent, -1, kFirstDaughter, kLastDaughter,
554// px, py, pz, e, vx, vy, vz, tof);
555
556 if (GetDebug()) Info("GetOneD","Randomized qinv = %f, obtained = %f",qinv,GetQInv(first,second));
557
558}
559/***********************************************************/
560
561Int_t AliGenHBTosl::GetThreeD(TParticle* first,TParticle* second, Double_t qout, Double_t qside, Double_t qlong)
562{
563//deprecated method that caclulates momenta of the second particle
564//out of qout qside and qlong and the first particle
565 Double_t m = first->GetMass();
566 Double_t m2 = m*m;
567
568 Double_t px = first->P();//first particle is rotated that only X is non-zero thus P==Px
569 Double_t px2 = px*px;
570
571
572 Double_t qout2 = qout*qout;
573 Double_t qside2 = qside*qside;
574 Double_t qlong2 = qlong*qlong;
575
576
577 Double_t util1 = 4.*px2 - qside2;//4*P1x^2 - Y^2
578 if (util1 < 0)
579 {
580 Info("GetThreeD","4.*px2* - qside2 is negative px: %f, qside: %f",px,qside);
581 return 1;
582 }
583 Double_t util2 = TMath::Sqrt(px2*qout2*util1);
584
585
586 Double_t p2x,p2y,p2z;
587
588// if ( (qside >= 0) && (qout >= 0) && (qlong >= 0))
589 if (qout >= 0)
590 {
591 //p2x
592 Double_t tmp = px*(2.*px2 - qside2);
593 tmp -= util2;
594 p2x = tmp/(2.*px2);
595
596 //p2y
597 tmp = qout - TMath::Sqrt(util1);
598 p2y = - (tmp*qside)/(2.*px);
599
600 //p2z
601 tmp = 4.*m2 + 2.*qout2+qlong2;
602 tmp *= px;
603 tmp -= 2.*util2;//!!!
604 tmp += 4.*px*px2;
605 tmp *= qlong2;
606
607 Double_t m2px2 = m2+px2;
608 Double_t tmp2 = m2px2*tmp;
609
610 tmp = 4.*(m2px2+qout2) + qlong2;
611 tmp *= px;
612 tmp -= 4.*util2;
613 tmp *= 4.*(m2px2) + qlong2;
614 tmp *= qlong2*qlong2;
615 tmp *= m2px2*m2px2;
616 tmp *= px;
617 if (tmp < 0)
618 {
619 Error("","Argument of sqrt is negative");
620 return 1;
621 }
622
623 tmp2 += TMath::Sqrt(tmp);
624
625 tmp = 8.0*px*m2px2*m2px2;
626 p2z = -TMath::Sqrt(tmp2/tmp);
627 if (qlong < 0) p2z = -p2z;
628 }
629 else
630 {
631 //p2x
632 Double_t tmp = px*(2.*px2 - qside2);
633 tmp += util2;
634 p2x = tmp/(2.*px2);
635
636 //p2y
637 tmp = qout - TMath::Sqrt(util1);
638 p2y = - (tmp*qside)/(2.*px);
639
640 //p2z
641 tmp = 4.*m2 + 2.*qout2+qlong2;
642 tmp *= px;
643 tmp += 2.*util2;//!!!
644 tmp += 4.*px*px2;
645 tmp *= qlong2;
646
647 Double_t m2px2 = m2+px2;
648 Double_t tmp2 = m2px2*tmp;
649
650 tmp = 4.*(m2px2+qout2) + qlong2;
651 tmp *= px;
652 tmp += 4.*util2;
653 tmp *= 4.*(m2px2) + qlong2;
654 tmp *= qlong2*qlong2;
655 tmp *= m2px2*m2px2;
656 tmp *= px;
657 if (tmp < 0)
658 {
659 Error("","Argument of sqrt is negative");
660 return 1;
661 }
662
663 tmp2 += TMath::Sqrt(tmp);
664
665 tmp = 8.0*px*m2px2*m2px2;
666 p2z = -TMath::Sqrt(tmp2/tmp);
667 if (qlong < 0) p2z = -p2z;
668 }
669
670// if ( (qside >= 0) && (qout >= 0) && (qlong >= 0)) p2z = -p2z;
671
672 TVector3 firstpvector(first->Px(),first->Py(),first->Pz());
673 TVector3 vector(p2x,p2y,p2z);
674 Rotate(firstpvector,vector);
675
676 Double_t e = TMath::Sqrt(m2 + vector.X()*vector.X() + vector.Y()*vector.Y() + vector.Z()*vector.Z());
677 second->SetMomentum(vector.X(),vector.Y(),vector.Z(),e);
678
679////////////
680 if ( GetDebug() > 3 )
681 {
682 e=TMath::Sqrt(m2+px*px);
683 TParticle* f = new TParticle(first->GetPdgCode(),0,-1,-1,-1,-1, px , 0.0, 0.0, e,0.0,0.0,0.0,0.0);
684
685 e = TMath::Sqrt(m2 + p2x*p2x + p2y*p2y + p2z*p2z);
686 TParticle* s = new TParticle(first->GetPdgCode(),0,-1,-1,-1,-1, p2x, p2y, p2z, e, 0.0, 0.0, 0.0, 0.0);
687
688 Double_t qo, qs, ql;
689 GetQOutQSideQLong(f,s,qo, qs, ql);
690
691 Info("GetThreeD","TEST");
692 f->Print();
693 s->Print();
694 Info("GetThreeD","Required %f %f %f",qout,qside,qlong);
695 Info("GetThreeD","Got %f %f %f",qo,qs,ql);
696 if ( qout == qo)
697 if (qside == qs)
698 if (qlong == ql)
699 Info("GetThreeD","!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
700 }
701/////////////
702 return 0;
703}
704/***********************************************************/
705
706void AliGenHBTosl::StartSignal()
707{
708//Starts the signal histograms
709 AliStack* stack;
710
711 fSwapped = kFALSE;
712
713 TParticle particle(fPID,0,-1,-1,-1,-1,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);
714 TParticle* second = &particle;
715
716 TIter next(fStackBuffer);
717 while(( stack=(AliStack*)next() ))
718 {
719 stack->Reset();
720 }
721
722 AliStack* genstack = fGenerator->GetStack();
723 if (genstack == 0x0)
724 {
725 genstack = new AliStack(fNpart);
726 fGenerator->SetStack(genstack);
727 }
728 else
729 {
730 genstack->Reset();
731 }
732
733 const Double_t ndf = fQNBins*fQNBins*fQNBins;
734
735 TH3D* work = new TH3D("work","work",fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange);
736 work->Sumw2();
737 work->SetDirectory(0x0);
738
739
740 Double_t binwdh = work->GetBinWidth(1)/2.;
741
742 Double_t*** chiarray = new Double_t** [fQNBins+1];
743 Double_t*** sigarray = new Double_t** [fQNBins+1];
744
745 for (Int_t i = 1; i<=fQNBins; i++)
746 {
747 chiarray[i] = new Double_t* [fQNBins+1];
748 sigarray[i] = new Double_t* [fQNBins+1];
749
750 for (Int_t k = 1; k<=fQNBins; k++)
751 {
752 chiarray[i][k] = new Double_t [fQNBins+1];
753 sigarray[i][k] = new Double_t [fQNBins+1];
754 }
755 }
756
757
758 Float_t chisqrchange = fMaxChiSquereChange + 1.;
759 Float_t chisqrPerDF = fMaxChiSquerePerNDF;
760 Float_t chisqrold = 0.0;
761
762 Int_t counter = 1;
763 Int_t niterations = 1;
764 Int_t rotaxisorder = 1;//defines order of looping over 3D histo (X,Y,Z or Y,Z,X or Z,X,Y)
765
766 Bool_t flag = kTRUE;
767 Bool_t shortloop = kTRUE;
768 TCanvas* c1 = new TCanvas();
769
770 StartSignalPass1();
771
772 printf("\n");
773 Info("StartSignal","\n\n\n\nSecond Pass\n\n\n\n");
774
775 while ( ( (chisqrPerDF > fMaxChiSquereChange) || flag) && (niterations++ < fMaxIterations) )
776 {
777
778 printf("StartSignal\n");
779 printf(" Row 1 Theory, 2 current value, 3 Chi2 \n");
780
781 Double_t chisqrnew = 0.0;
782 flag = kFALSE;
783
784 Double_t scale = Scale(fQSignal,fQBackground);
785 work->Divide(fQSignal,fQBackground,scale);
786
787 if ( (counter%100) == 0)
788 {
789 c1->cd();
790 char buff[50];
791 sprintf(buff,"QTWorkPass2.%3d.root",counter);
792 TFile* file = TFile::Open(buff,"update");
793 work->Write();
794 work->SetDirectory(0x0);
795 file->Close();
796 delete file;
797
798 sprintf(buff,"QTBackgroundPass2.%3d.root",counter);
799 file = TFile::Open(buff,"update");
800 fQBackground->Write();
801 fQBackground->SetDirectory(0x0);
802 file->Close();
803 delete file;
804
805 sprintf(buff,"QTSignalPass2.%3d.root",counter);
806 file = TFile::Open(buff,"update");
807 fQSignal->Write();
808 fQSignal->SetDirectory(0x0);
809 file->Close();
810 delete file;
811 }
812
813 counter++;
814 Int_t novertresh = 0;
815 for (Int_t k = 1; k<=fQNBins; k++)
816 {
817 Double_t z = work->GetZaxis()->GetBinCenter(k);
818 for (Int_t j = 1; j<=fQNBins; j++)
819 {
820 Double_t y = work->GetYaxis()->GetBinCenter(j);
821 for (Int_t i = 1; i<=fQNBins; i++)
822 {
823 Double_t x = work->GetXaxis()->GetBinCenter(i);//get center value of a bin (qout)
824 sigarray[i][j][k] = fQSignal->GetBinContent(i,j,k);//store current value of signal histogram
825 Double_t v = GetQOutQSideQLongCorrTheorValue(x,y,z);//get expected value of CF in that qinv
826 Double_t diff = v - work->GetBinContent(i,j,k);//store difference betweeon current value, and desired value
827 chiarray[i][j][k] = diff; // no-x x is a weight to get good distribution
828 Double_t be = work->GetBinError(i,j,k);
829 chisqrnew += diff*diff/(be*be);//add up chisq
830
831 //even if algorithm is stable (chi sqr change less then threshold)
832 //and any bin differs more then 5% from expected value we continue
833 Double_t fact = diff;
834 if (TMath::Abs(fact) > 0.1)
835 {
836 flag = kTRUE;
837 novertresh++;
838 }
839 }
840 }
841 }
842
843
844 printf("\n");
845
846 for (Int_t k = 1; k < 11; k++)
847 {
848 Double_t tx = work->GetXaxis()->GetBinCenter(1);
849 Double_t ty = work->GetYaxis()->GetBinCenter(1);
850 Double_t tz = work->GetZaxis()->GetBinCenter(k);
851 printf("% 6.5f ",GetQOutQSideQLongCorrTheorValue(tx,ty,tz));
852 }
853 printf("\n");
854
855 for (Int_t k = 1; k < 11; k++)
856 {
857 printf("% 6.5f ",work->GetBinContent(1,1,k));
858 }
859 printf("\n");
860
861 for (Int_t k = 1; k < 11; k++)
862 {
863 printf("% 6.5f ",chiarray[1][1][k]);
864 }
865 printf("\n");
866
867 chisqrchange = TMath::Abs(chisqrnew - chisqrold)/chisqrnew;
868 chisqrold = chisqrnew;
869
870 chisqrPerDF = chisqrnew/ndf;
871
872 Info("StartSignal","Iteration %d Chi-sq change %f%%",niterations,chisqrchange*100.0);
873 Info("StartSignal","ChiSq = %f, NDF = %f, ChiSq/NDF = %f",chisqrnew, ndf, chisqrPerDF );
874 Info("StartSignal","novertresh = %d",novertresh);
875
876
877 stack = RotateStack();
878 genstack->Reset();
879 fGenerator->Generate();
880 Int_t ninputparticle = 0, ntr = 0;
881 if ( genstack->GetNtrack() < fNpart/2)
882 {
883 Warning("StartSignal","**********************************");
884 Warning("StartSignal","Generator generated (%d) less particles then expected (%d).",
885 genstack->GetNtrack(),fNpart/2);
886 Warning("StartSignal","**********************************");
887 }
888
889 Int_t sc = 0; //security check against infinite loop
890 while ( (ntr+1) < fNpart)//ntr is number of track generated up to now
891 {
892 Int_t xmax = 1;
893 Int_t ymax = 1;
894 Int_t zmax = 1;
895 Double_t qout;
896 Double_t qside;
897 Double_t qlong;
898
899 Int_t i,j,k;
900
901 Int_t* cx = 0x0;
902 Int_t* cy = 0x0;
903 Int_t* cz = 0x0;
904
905 switch (rotaxisorder)
906 {
907 case 1:
908 cx = &i;
909 cy = &j;
910 cz = &k;
911 break;
912 case 2:
913 cx = &j;
914 cy = &k;
915 cz = &i;
916 break;
917 case 3:
918 cx = &k;
919 cy = &i;
920 cz = &j;
921 break;
922 }
923 rotaxisorder++;
924 if (rotaxisorder > 3) rotaxisorder = 1;
925 Int_t nrange;
926
927 if (shortloop)
928 {
929 shortloop = kFALSE;
930 nrange = fQNBins;
931 }
932 else
933 {
934 shortloop = kTRUE;
935 nrange = fQNBins/4;
936 }
937
938// Bool_t force = kFALSE;
939 for ( k = 1; k <=nrange;k++ )
940 {
941 for ( j = 1; j<=nrange; j++)
942 {
943 for ( i = 1; i<=nrange; i++)
944 {
945 if ( (chiarray[*cx][*cy][*cz]) > (chiarray[xmax][ymax][zmax]) )
946 {
947 xmax = *cx;
948 ymax = *cy;
949 zmax = *cz;
950 }
951
952// Double_t fact = chiarray[*cx][*cy][*cz];//chiarray is chi2*qinv
953// if (fact > work->GetBinError(*cx,*cy,*cz))//if differece between what we want and
954// { //what we have is bigger than stat. error
955// //we force to fill that bin
956// qout = work->GetXaxis()->GetBinCenter(*cx);
957// qside = work->GetYaxis()->GetBinCenter(*cy);
958// qlong = work->GetZaxis()->GetBinCenter(*cz);
959
960// Info("StartSignal"," bin: (%d,%d,%d) loop status (%d,%d,%d) \nUsing Force: chiarray: %f \nq(o,s,l): (%f,%f,%f) signal: %d background: %d binerror: %f",
961// *cx,*cy,*cz,i,j,k,fact,qout,qside,qlong,
962// (Int_t)sigarray[*cx][*cy][*cz],(Int_t)fQBackground->GetBinContent(*cx,*cy,*cz),work->GetBinError(*cx,*cy,*cz));
963// force = kTRUE;
964// break;
965// }
966
967 }
968// if (force) break;
969 }
970// if (force) break;
971 }
972
973
974 qout = work->GetXaxis()->GetBinCenter(xmax);
975 qside = work->GetYaxis()->GetBinCenter(ymax);
976 qlong = work->GetZaxis()->GetBinCenter(zmax);
977
978// Info("StartSignal"," bin: (%d,%d,%d) chiarray: %f \nq(o,s,l): (%f,%f,%f) signal: %d background: %d binerror: %f",
979// xmax,ymax,zmax,chiarray[xmax][ymax][zmax],qout,qside,qlong,
980// (Int_t)sigarray[xmax][ymax][zmax],
981// (Int_t)fQBackground->GetBinContent(xmax,ymax,zmax),
982// work->GetBinError(xmax,ymax,zmax));
983
984 qout = gRandom->Uniform(qout-binwdh,qout+binwdh);
985 qside = gRandom->Uniform(qside-binwdh,qside+binwdh);
986 qlong = gRandom->Uniform(qlong-binwdh,qlong+binwdh);
987
988 TParticle* first = 0;
989 while (ninputparticle < genstack->GetNtrack())
990 {
991 TParticle* tmpp = genstack->Particle(ninputparticle++);
992 if (tmpp->GetPdgCode() == fPID)
993 {
994 if (CheckParticle(tmpp,0x0,stack) == kFALSE)
995 {
996 first = tmpp;
997 break;
998 }
999 }
1000 }
1001
1002 if (first == 0x0)
1003 {
1004 if ( fDebug > 2 ) Info("StartSignal","No more particles of that type");
1005 break;
1006 }
1007
1008 Int_t retval = GetThreeD(first,second,qout,qside,qlong);
1009 if (retval)
1010 {
1011 Info("StartSignal","Can not find momenta for this OSL and particle OSL = %f %f %f",qout,qside,qlong);
1012 first->Print();
1013 second->Print();
1014
1015 continue;
1016 }
1017 //in case this particle is falling into signal area with another
1018 //particle we take a another pair
1019 //it can intruduce artificial correlations
1020 Bool_t checkresult = CheckParticle(second,first,stack);
1021 if ( checkresult && (sc < 10) )
1022 {
1023 sc++;
1024 continue;
1025 }
1026 sc = 0;
1027
1028 //Put on output stack
1029 SetTrack(first,ntr,stack);
1030 SetTrack(second,ntr,stack);
1031
1032 Double_t y = GetQOutQSideQLongCorrTheorValue(qout,qside,qlong);
1033
1034 sigarray[xmax][ymax][zmax] ++;
1035 chiarray[xmax][ymax][zmax] = scale*sigarray[xmax][ymax][zmax]/fQBackground->GetBinContent(xmax,ymax,zmax);
1036 chiarray[xmax][ymax][zmax] = (y - chiarray[xmax][ymax][zmax]);
1037
1038 }
1039 Info("StartSignal","Mixing background...");
1040 Mix(fStackBuffer,fQBackground,fQSecondBackground); //upgrate background
1041 Info("StartSignal","Mixing signal...");
1042 Mix(stack,fQSignal,fQSecondSignal); //upgrate background
1043
1044 if ( (chisqrPerDF < 2.0) && (fSwapped == kFALSE) )
1045 {
1046 SwapGeneratingHistograms();
1047 }
1048
1049 }
1050 TFile* filef = TFile::Open("QTBackground.root","recreate");
1051 fQBackground->Write();
1052 fQBackground->SetDirectory(0x0);
1053 filef->Close();
1054 delete filef;
1055
1056 filef = TFile::Open("QTSignal.root","recreate");
1057 fQSignal->Write();
1058 fQSignal->SetDirectory(0x0);
1059 filef->Close();
1060 delete filef;
1061
1062
1063 delete c1;
1064 delete work;
1065
1066 for (Int_t i = 1; i<=fQNBins; i++)
1067 {
1068 for (Int_t k = 1; k<=fQNBins; k++)
1069 {
1070 delete [] chiarray[i][k];
1071 delete [] sigarray[i][k];
1072 }
1073 delete [] chiarray[i];
1074 delete [] sigarray[i];
1075 }
1076 delete [] chiarray;
1077 delete [] sigarray;
1078
1079}
1080/***********************************************************/
1081
1082void AliGenHBTosl::StartSignalPass1()
1083{
1084 //This method makes first part of the initialization of working histograms
1085 //It randomizes qout, qside and qlong from the coarse signal histogram
1086
1087 Bool_t flag = kTRUE;
1088 TParticle particle(fPID,0,-1,-1,-1,-1,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);
1089 TParticle* second = &particle;
1090 Double_t qout;
1091 Double_t qside;
1092 Double_t qlong;
1093
1094 Info("StartSignalPass1","\n\nFirst Pass\n\n");
1095
1096 while (flag)
1097 {
1098 Info("StartSignalPass1","NextEvent");
1099 AliStack* stack = RotateStack();
1100 AliStack* genstack = fGenerator->GetStack();
1101 genstack->Reset();
1102 fGenerator->Generate();
1103 Int_t j = 0, ntr = 0;
1104 if ( genstack->GetNtrack() < fNpart/2)
1105 {
1106 Warning("StartSignalPass1","**********************************");
1107 Warning("StartSignalPass1","Generator generated (%d) less particles then expected (%d).",
1108 genstack->GetNtrack(),fNpart/2);
1109 Warning("StartSignalPass1","**********************************");
1110 }
1111
1112 Int_t sc = 0;//security check against infinite loop
1113 while ((ntr+1)<fNpart)
1114 {
1115
1116// Info("StartSignal","Number of track on output stack: = %d", ntr);
1117// Info("StartSignal","Number of track on input stack: = %d\n", j);
1118
1119 TParticle* first = 0;
1120 while (j < genstack->GetNtrack())
1121 {
1122 TParticle* tmpp = genstack->Particle(j++);
1123 if (tmpp->GetPdgCode() == fPID)
1124 {
1125 if (CheckParticle(tmpp,0x0,stack) == kFALSE)
1126 {
1127 first = tmpp;
1128 break;
1129 }
1130 else
1131 {
1132 Info("StartSignalPass1","Particle did not pass the safety check 1");
1133 tmpp->Print();
1134 }
1135 }
1136 }
1137
1138 if (first == 0x0)
1139 {
1140 if ( fDebug > 2 ) Info("StartSignalPass1","No more particles of that type");
1141
1142 break;
1143 }
1144
1145 SetTrack(first,ntr,stack);
1146
1147 fQCoarseSignal->GetRandom3(qout,qside,qlong);
1148
1149 Int_t retval = GetThreeD(first,second,qout,qside,qlong);
1150 if (retval)
1151 {
1152 //Info("StartSignal","Can not find momenta for this OSL and particle");
1153 continue;
1154 }
1155 //in case this particle is falling into signal area with another
1156 //particle we take a another pair
1157 //it can intruduce artificial correlations
1158 Bool_t checkresult = CheckParticle(second,first,stack);
1159 if ( checkresult && (sc < 10) )
1160 {
1161 sc++;
1162 Info("StartSignalPass1","Particle did not pass the safety check 2");
1163 second->Print();
1164 continue;
1165 }
1166
1167 sc = 0;
1168
1169 SetTrack(second,ntr,stack);
1170 }
1171
1172 Mix(stack,fQSignal,fQSecondSignal);
1173 Mix(fStackBuffer,fQBackground,fQSecondBackground);
1174
1175 flag = kFALSE;
1176
1177 for (Int_t k = 1; k<=fQNBins; k++)
1178 {
1179 for (Int_t j = 1; j<=fQNBins; j++)
1180 {
1181 for (Int_t i = 1; i<=fQNBins; i++)
1182 {
1183 if ( (fQBackground->GetBinContent(i,j,k) < fMinFill) )
1184 {
1185 //(fQSignal->GetBinContent(i,j,k) < fMinFill) ||
1186 Info("StartSignalPass1","bin (%d,%d,%d): signal=%f background=%f",i,j,k,
1187 fQSignal->GetBinContent(i,j,k),fQBackground->GetBinContent(i,j,k));
1188 flag = kTRUE;//continue while
1189 break;//breakes for not while
1190 }
1191 }
1192 if (flag == kTRUE) break;
1193 }
1194 if (flag == kTRUE) break;
1195 }
1196
1197 }//while (flag)
1198
1199
1200}
1201/***********************************************************/
1202
1203void AliGenHBTosl::FillCoarseSignal()
1204{
1205 //Makes coarse signal by multiplying the coarse background and required function
1206 Info("FillCoarseSignal","START");
1207 for (Int_t k = 1; k <=fQNBins ;k++ )
1208 {
1209 Double_t z = fQCoarseBackground->GetZaxis()->GetBinCenter(k);
1210 for (Int_t j = 1; j <=fQNBins; j++)
1211 {
1212 Double_t y = fQCoarseBackground->GetYaxis()->GetBinCenter(j);
1213 for (Int_t i = 1; i <=fQNBins; i++)
1214 {
1215 Double_t x = fQCoarseBackground->GetXaxis()->GetBinCenter(i);
1216 Double_t v = GetQOutQSideQLongCorrTheorValue(x,y,z);
1217 Info("FillCoarseSignal","Bin (%d,%d,%d): osl(%f,%f,%f)=%f",i,j,k,x,y,z,v);
1218 fQCoarseSignal->SetBinContent(i,j,k,v*fQCoarseBackground->GetBinContent(i,j,k));
1219 }
1220 }
1221 }
1222
1223 //if (GetDebug())
1224 TestCoarseSignal();
1225
1226 Info("FillCoarseSignal","DONE");
1227}
1228/***********************************************************/
1229
1230void AliGenHBTosl::FillCoarse()
1231{
1232 //creates the statistical background histogram on the base of input from
1233 //fGenerator
1234 Info("FillCoarse","START");
1235
1236 AliStack* stack;
1237 Int_t niter = 0;
1238
1239 Bool_t cont;
1240 TH3D tmph("tmph","tmph",2,0,1,2,0,1,2,0,1);
1241 printf("\n");
1242
1243 do
1244 {
1245// if (niter > 20) break;
1246
1247 cout<<niter++<<" bincont "<<fQCoarseBackground->GetBinContent(1,1,1)
1248 <<" "<<fQCoarseBackground->GetBinContent(1,1,2)
1249 <<" "<<fQCoarseBackground->GetBinContent(1,1,3)
1250 <<" "<<fQCoarseBackground->GetBinContent(1,1,4)
1251 <<" "<<fQCoarseBackground->GetBinContent(1,1,5)
1252 <<"\n";
1253 fflush(0);
1254
1255 stack = RotateStack();
1256 fGenerator->SetStack(stack);
1257 fGenerator->Init();
1258 fGenerator->Generate();
1259
1260 Mix(fStackBuffer,fQCoarseBackground,&tmph);
1261
1262 cont = kFALSE;
1263
1264 Info("FillCoarse","fMinFill = %d",fMinFill);
1265 for (Int_t k = 1; k<=fQNBins; k++)
1266 {
1267 for (Int_t j = 1; j<=fQNBins; j++)
1268 {
1269 for (Int_t i = 1; i<=fQNBins; i++)
1270 {
1271 if ( fQCoarseBackground->GetBinContent(i,j,k) < fMinFill )
1272 {
1273 cont = kTRUE;
1274 Info("FillCoarse","bin (%d,%d,%d)=%f",i,j,k,fQCoarseBackground->GetBinContent(i,j,k));
1275 break;
1276 }
1277
1278 }
1279 if (cont) break;
1280 }
1281 if (cont) break;
1282 }
1283 }while(cont);
1284
1285 printf("\n");
1286 fGenerator->SetStack(0x0);
1287 Info("FillCoarse","DONE");
1288
1289}
1290/***********************************************************/
1291
1292void AliGenHBTosl::Mix(TList* eventbuffer,TH3D* denominator,TH3D* denominator2)
1293{
1294 //Fills denominators
1295 //Mixes events stored in the eventbuffer and fills the background histograms
1296 static TStopwatch stoper;
1297
1298 if (eventbuffer == 0x0)
1299 {
1300 Error("Mix","Buffer List is null.");
1301 return;
1302 }
1303
1304 if (denominator == 0x0)
1305 {
1306 Error("Mix","Denominator histogram is null.");
1307 return;
1308 }
1309
1310 if (denominator2 == 0x0)
1311 {
1312 Error("Mix","Denominator2 histogram is null.");
1313 return;
1314 }
1315
1316 Info("Mix","%s",denominator->GetName());
1317 stoper.Start();
1318
1319 TIter next(eventbuffer);
1320 AliStack* firstevent;
1321 AliStack* secondevent = 0x0;
1322
1323 while(( firstevent=(AliStack*)next() ))
1324 {
1325 if (secondevent == 0x0)
1326 {
1327 secondevent = firstevent;
1328 continue;
1329 }
1330// Info("Mix","Mixing %#x with %#x",firstevent,secondevent);
1331 for(Int_t j = 0; j < firstevent->GetNtrack(); j++ )
1332 {
1333 TParticle* firstpart = firstevent->Particle(j);
1334
1335 Float_t phi = firstpart->Phi();
1336 if ( (phi < fSamplePhiMin) || ( phi > fSamplePhiMax)) continue;
1337
1338// Info("Mix","Mixing %d phi %f min %f max %f",j,phi,fSamplePhiMin,fSamplePhiMax);
1339
1340 for(Int_t i = 0; i < secondevent->GetNtrack(); i++ )
1341 {
1342 TParticle* secondpart = secondevent->Particle(i);
1343 phi = secondpart->Phi();
1344 if ( (phi < fSamplePhiMin) || ( phi > fSamplePhiMax)) continue;
1345
1346 Double_t qout;
1347 Double_t qside;
1348 Double_t qlong;
1349 GetQOutQSideQLong(firstpart,secondpart,qout,qside,qlong);
1350 denominator->Fill(qout,qside,qlong);
1351 denominator2->Fill(qout,qside,qlong);
1352 }
1353 }
1354
1355 secondevent = firstevent;
1356 }
1357 stoper.Stop();
1358 stoper.Print();
1359
1360}
1361/***********************************************************/
1362
1363void AliGenHBTosl::Mix(AliStack* stack, TH3D* numerator, TH3D* numerator2)
1364{
1365//fils numerator with particles from stack
1366 static TStopwatch stoper;
1367 if (stack == 0x0)
1368 {
1369 Error("Mix","Stack is null.");
1370 return;
1371 }
1372
1373 if ( (numerator == 0x0) || (numerator2 == 0x0) )
1374 {
1375 Error("Mix","Numerator histogram is null.");
1376 return;
1377 }
1378
1379 Info("Mix","%s",numerator->GetName());
1380 stoper.Start();
1381
1382 for(Int_t j = 0; j < stack->GetNtrack(); j++ )
1383 {
1384 TParticle* firstpart = stack->Particle(j);
1385 Float_t phi = firstpart->Phi();
1386 if ( (phi < fSamplePhiMin) || ( phi > fSamplePhiMax)) continue;
1387
1388 for(Int_t i = j+1; i < stack->GetNtrack(); i++ )
1389 {
1390 TParticle* secondpart = stack->Particle(i);
1391 phi = secondpart->Phi();
1392 if ( (phi < fSamplePhiMin) || ( phi > fSamplePhiMax)) continue;
1393 Double_t qout;
1394 Double_t qside;
1395 Double_t qlong;
1396 GetQOutQSideQLong(firstpart,secondpart,qout,qside,qlong);
1397 numerator->Fill(qout,qside,qlong);
1398 numerator2->Fill(qout,qside,qlong);
1399 }
1400 }
1401 stoper.Stop();
1402 stoper.Print();
1403
1404}
1405/***********************************************************/
1406
1407Double_t AliGenHBTosl::GetQInv(TParticle* f, TParticle* s)
1408{
1409//calculates qinv
1410// cout<<f->Px()<<" "<<s->Px()<<endl;
1411 Double_t pxdiff = f->Px() - s->Px();
1412 Double_t pydiff = f->Py() - s->Py();
1413 Double_t pzdiff = f->Pz() - s->Pz();
1414 Double_t ediff = f->Energy() - s->Energy();
1415
1416 Double_t qinvl = ediff*ediff - ( pxdiff*pxdiff + pydiff*pydiff + pzdiff*pzdiff );
1417 Double_t qinv = TMath::Sqrt(TMath::Abs(qinvl));
1418 return qinv;
1419}
1420/***********************************************************/
1421
1422void AliGenHBTosl::GetQOutQSideQLong(TParticle* f, TParticle* s,Double_t& out, Double_t& side, Double_t& lon)
1423{
1424 //returns qout,qside and qlong of the pair of particles
1425 out = side = lon = 10e5;
1426
1427 Double_t pxsum = f->Px() + s->Px();
1428 Double_t pysum = f->Py() + s->Py();
1429 Double_t pzsum = f->Pz() + s->Pz();
1430 Double_t esum = f->Energy() + s->Energy();
1431 Double_t pxdiff = f->Px() - s->Px();
1432 Double_t pydiff = f->Py() - s->Py();
1433 Double_t pzdiff = f->Pz() - s->Pz();
1434 Double_t ediff = f->Energy() - s->Energy();
1435 Double_t kt = 0.5*TMath::Hypot(pxsum,pysum);
1436
1437 Double_t k2 = pxsum*pxdiff+pysum*pydiff;
1438
1439 if (kt == 0.0)
1440 {
1441 f->Print();
1442 s->Print();
1443 kt = 10e5;
1444 }
1445 else
1446 {
1447 out = 0.5*k2/kt;
1448 side = (f->Px()*s->Py()-s->Px()*f->Py())/kt;
1449 }
1450
1451 Double_t beta = pzsum/esum;
1452 Double_t gamma = 1.0/TMath::Sqrt(1.0 - beta*beta);
1453
1454 lon = gamma * ( pzdiff - beta*ediff );
1455
1456// out = TMath::Abs(out);
1457// side = TMath::Abs(side);
1458// lon = TMath::Abs(lon);
1459}
1460
1461/***********************************************************/
1462
1463Double_t AliGenHBTosl::Scale(TH3D* num, TH3D* den)
1464{
1465 //Calculates the factor that should be used to scale
1466 //quatience of num and den to 1 at tail
1467
1468 if (GetDebug()) Info("Scale","Entered Scale()");
1469 if(!num)
1470 {
1471 Error("Scale","No numerator");
1472 return 0.0;
1473 }
1474 if(!den)
1475 {
1476 Error("Scale","No denominator");
1477 return 0.0;
1478 }
1479
1480 if(fNBinsToScale < 1)
1481 {
1482 return 0.0;
1483 Error("Scale","Number of bins for scaling is smaller thnan 1");
1484 }
1485 Int_t fNBinsToScaleX = fNBinsToScale;
1486 Int_t fNBinsToScaleY = fNBinsToScale;
1487 Int_t fNBinsToScaleZ = fNBinsToScale;
1488
1489 Int_t nbinsX = num->GetNbinsX();
1490 if (fNBinsToScaleX > nbinsX)
1491 {
1492 Error("Scale","Number of X bins for scaling is bigger thnan number of bins in histograms");
1493 return 0.0;
1494 }
1495
1496 Int_t nbinsY = num->GetNbinsX();
1497 if (fNBinsToScaleY > nbinsY)
1498 {
1499 Error("Scale","Number of Y bins for scaling is bigger thnan number of bins in histograms");
1500 return 0.0;
1501 }
1502
1503 Int_t nbinsZ = num->GetNbinsZ();
1504 if (fNBinsToScaleZ > nbinsZ)
1505 {
1506 Error("Scale","Number of Z bins for scaling is bigger thnan number of bins in histograms");
1507 return 0.0;
1508 }
1509
1510 if (GetDebug()>0) Info("Scale","No errors detected");
1511
1512 Int_t offsetX = nbinsX - fNBinsToScaleX - 1; //bin that we start loop over bins in axis X
1513 Int_t offsetY = nbinsY - fNBinsToScaleY - 1; //bin that we start loop over bins in axis Y
1514 Int_t offsetZ = nbinsZ - fNBinsToScaleZ - 1; //bin that we start loop over bins in axis Z
1515
1516 Double_t densum = 0.0;
1517 Double_t numsum = 0.0;
1518
1519 for (Int_t k = offsetZ; k<nbinsZ; k++)
1520 for (Int_t j = offsetY; j<nbinsY; j++)
1521 for (Int_t i = offsetX; i<nbinsX; i++)
1522 {
1523 if ( num->GetBinContent(i,j,k) > 0.0 )
1524 {
1525
1526 densum += den->GetBinContent(i,j,k);
1527 numsum += num->GetBinContent(i,j,k);
1528 }
1529 }
1530
1531 if(GetDebug() > 0)
1532 Info("Scale","numsum=%f densum=%f fNBinsToScaleX=%d fNBinsToScaleY=%d fNBinsToScaleZ=%d",
1533 numsum,densum,fNBinsToScaleX,fNBinsToScaleY,fNBinsToScaleZ);
1534
1535 if (numsum == 0) return 0.0;
1536 Double_t ret = densum/numsum;
1537
1538 if(GetDebug() > 0) Info("Scale","returning %f",ret);
1539 return ret;
1540
1541}
1542/***********************************************************/
1543
1544void AliGenHBTosl::TestCoarseSignal()
1545{
1546//Tests how works filling from generated histogram shape
1547 TH3D* work = new TH3D("work","work",fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange);
1548
1549// for (Int_t i = 0; i < fQCoarseBackground->GetEntries() ;i++)
1550// {
1551// Double_t x,y,z;
1552// fQCoarseSignal->GetRandom3(x,y,z);
1553// work->Fill(x,y,z);
1554// }
1555
1556 TCanvas* c1 = new TCanvas();
1557 c1->cd();
1558 work->Draw();
1559 c1->SaveAs("QTwork.root");
1560 TFile* file = TFile::Open("QTwork.root","update");
1561// work->Write();
1562 work->SetDirectory(0x0);
1563 file->Close();
1564
1565 fQCoarseSignal->Draw();
1566 c1->SaveAs("QTCoarseSignal.root");
1567 file = TFile::Open("QTCoarseSignal.root","update");
1568 fQCoarseSignal->Write();
1569 fQCoarseSignal->SetDirectory(0x0);
1570 file->Close();
1571
1572 fQCoarseBackground->Draw();
1573 c1->SaveAs("QTCoarseBackground.root");
1574 file = TFile::Open("QTCoarseBackground.root","update");
1575 fQCoarseBackground->Write();
1576 fQCoarseBackground->SetDirectory(0x0);
1577 file->Close();
1578
1579 TH1 *result = (TH1*)fQCoarseBackground->Clone("ratio");
1580 result->SetTitle("ratio");
1581 Float_t normfactor = Scale(work,fQCoarseBackground);
1582 result->Divide(work,fQCoarseBackground,normfactor);//work
1583
1584
1585 c1->cd();
1586 result->Draw();
1587 c1->SaveAs("QTresult.root");
1588 file = TFile::Open("QTresult.root","update");
1589 result->Write();
1590 result->SetDirectory(0x0);
1591 file->Close();
1592
1593 delete work;
1594 delete c1;
1595}
1596/***********************************************************/
1597
1598void AliGenHBTosl::SetTrack(TParticle* p, Int_t& ntr)
1599{
1600//Shortcut to PushTrack(bla,bla,bla,bla.............)
1601 if (p->P() == 0.0)
1602 {
1603 Error("SetTrack(TParticle*,Int_t&)","Particle has zero momentum");
1604 return;
1605 }
1606
1607
1608 Int_t pdg = p->GetPdgCode();
1609 Double_t px = p->Px();
1610 Double_t py = p->Py();
1611 Double_t pz = p->Pz();
1612 Double_t e = p->Energy();
1613 Double_t vx = p->Vx();
1614 Double_t vy = p->Vy();
1615 Double_t vz = p->Vz();
1616 Double_t tof = p->T();
1617
1618 TVector3 pol;
1619 p->GetPolarisation(pol);
1620
1621 Double_t polx = pol.X();
1622 Double_t poly = pol.Y();
1623 Double_t polz = pol.Z();
1624 TMCProcess mech = AliGenCocktailAfterBurner::IntToMCProcess(p->GetUniqueID());
1625 Float_t weight = p->GetWeight();
1626
1627 AliGenerator::PushTrack(fTrackIt, -1, pdg, px, py, pz, e, vx, vy, vz, tof,polx, poly, polz, mech, ntr, weight);
1628}
1629/***********************************************************/
1630
1631void AliGenHBTosl::SetTrack(TParticle* p, Int_t& ntr, AliStack* stack)
1632{
1633//Shortcut to SetTrack(bla,bla,bla,bla.............)
1634 if (p->P() == 0.0)
1635 {
1636 Error("SetTrack(TParticle*,Int_t&,AliStack*)","Particle has zero momentum");
1637 return;
1638 }
1639
1640 Int_t pdg = p->GetPdgCode();
1641 Double_t px = p->Px();
1642 Double_t py = p->Py();
1643 Double_t pz = p->Pz();
1644 Double_t e = p->Energy();
1645
1646 stack->PushTrack(fTrackIt, -1, pdg, px, py, pz, e, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, kPPrimary, ntr,1,0);
1647}
1648/***********************************************************/
1649
1650void AliGenHBTosl::Rotate(TVector3& relvector, TVector3& vector)
1651{
1652//This method rotates vector about the angeles that are needed to rotate
1653//relvector from postion (firstPx,0,0) to its actual positon
1654//In other words: To make equations easier
1655
1656 static TVector3 first;
1657 if (GetDebug())
1658 {
1659 first.SetXYZ(relvector.x(),relvector.y(),relvector.z());
1660 }
1661
1662 Double_t firstPx = TMath::Sqrt( relvector.x()*relvector.x() +
1663 relvector.y()*relvector.y() +
1664 relvector.z()*relvector.z() );
1665
1666 Double_t rotAngleZ = -TMath::ATan2(relvector.y(),relvector.x());//calculating rot angles
1667 relvector.RotateZ(rotAngleZ);
1668 rotAngleZ = -rotAngleZ;
1669 Double_t rotAngleY = -TMath::ATan2(relvector.z(),relvector.x());
1670
1671 vector.RotateY(rotAngleY);
1672 vector.RotateZ(rotAngleZ);
1673
1674 if (GetDebug()>5)
1675 {
1676 TVector3 test(firstPx,0.0,0.0);
1677 test.RotateY(rotAngleY);
1678 test.RotateZ(rotAngleZ);
1679 ::Info("Rotate","Rotation test: px %f %f",first.x(),test.x());
1680 ::Info("Rotate","Rotation test: py %f %f",first.y(),test.y());
1681 ::Info("Rotate","Rotation test: pz %f %f",first.z(),test.z());
1682 }
1683}
1684/***********************************************************/
1685
1686Double_t AliGenHBTosl::Rotate(Double_t x,Double_t y,Double_t z)
1687{
1688//Rotates vector to base where only x - coordinate is no-zero, and returns that
1689
1690 Double_t xylength = TMath::Hypot(x,y);
1691 Double_t sinphi = -y/xylength;
1692 Double_t cosphi = x/xylength;
1693
1694 Double_t xprime = cosphi*x - sinphi*y;
1695 Double_t yprime = sinphi*x + cosphi*y;
1696
1697 TVector3 v(x,y,z);
1698 Double_t a1 = -TMath::ATan2(v.Y(),v.X());
1699
1700 ::Info("Rotate","Xpr = %f Ypr = %f",xprime,yprime);
1701 ::Info("Rotate","Calc sin = %f, and %f",sinphi,TMath::Sin(a1));
1702 ::Info("Rotate","Calc cos = %f, and %f",cosphi,TMath::Cos(a1));
1703
1704
1705 Double_t xprimezlength = TMath::Hypot(xprime,z);
1706
1707 Double_t sintheta = z/xprimezlength;
1708 Double_t costheta = xprime/xprimezlength;
1709
1710
1711 Double_t xbis = sintheta*z + costheta*(cosphi*x - sinphi*y);
1712
1713 ::Info("Rotate","Calculated rot %f, modulus %f",xbis,TMath::Sqrt(x*x+y*y+z*z));
1714 return xbis;
1715}
1716/***********************************************************/
1717
1718AliStack* AliGenHBTosl::RotateStack()
1719{
1720//swaps to next stack last goes to first and is reseted
1721
1722 AliStack* stack;
1723 if ( fStackBuffer->GetSize() >= fBufferSize )
1724 {
1725 stack = (AliStack*)fStackBuffer->Remove(fStackBuffer->Last());
1726 }
1727 else
1728 {
1729 stack = new AliStack(fNpart);
1730 }
1731
1732 fStackBuffer->AddFirst(stack);
1733 stack->Reset();
1734 return stack;
1735}
1736/***********************************************************/
1737
1738Double_t AliGenHBTosl::GetQInvCorrTheorValue(Double_t qinv) const
1739{
1740//Function (deprecated)
1741 static const Double_t factorsqrd = 0.197*0.197;//squared conversion factor SI<->eV
1742
1743 return 1.0 + 0.5*TMath::Exp(-qinv*qinv*fQRadius*fQRadius/factorsqrd);
1744}
1745/***********************************************************/
1746
1747Double_t AliGenHBTosl::GetQOutQSideQLongCorrTheorValue(Double_t& out, Double_t& side, Double_t& lon) const
1748{
1749 //Theoretical function. Wa want to get correlation of the shape of this function
1750 static const Double_t factorsqrd = 0.197*0.197;//squared conversion factor SI<->eV
1751 return 1.0 + 0.7*TMath::Exp(-fQRadius*fQRadius*(out*out+side*side+lon*lon)/factorsqrd);
1752}
1753/***********************************************************/
1754
1755Bool_t AliGenHBTosl::CheckParticle(TParticle* p, TParticle* aupair ,AliStack* stack)
1756{
1757 //Checks if a given particle is falling into signal region with any other particle
1758 //already existing on stack
1759 return kFALSE;
1760
1761 if (fSignalRegion <=0) return kFALSE;
1762
1763 for (Int_t i = 0; i < stack->GetNtrack(); i++)
1764 {
1765 TParticle* part = stack->Particle(i);
1766 if (part == aupair) continue;
1767 Double_t qout = 10e5;
1768 Double_t qside= 10e5;
1769 Double_t qlong= 10e5;
1770 GetQOutQSideQLong(p,part,qout,qside,qlong);
1771
1772 if (TMath::Abs(qout) < fSignalRegion)
1773 if (TMath::Abs(qside) < fSignalRegion)
1774 if (TMath::Abs(qlong) < fSignalRegion)
1775 return kTRUE;
1776 }
1777 return kFALSE;
1778}
1779/***********************************************************/
1780
1781void AliGenHBTosl::SwapGeneratingHistograms()
1782{
1783 //Checks if it is time to swap signal and background histograms
1784 //if yes it swaps them
1785 Int_t threshold = fMinFill;
1786 for (Int_t k = 1; k<=fQNBins; k++)
1787 {
1788 for (Int_t j = 1; j<=fQNBins; j++)
1789 {
1790 for (Int_t i = 1; i<=fQNBins; i++)
1791 {
1792 if ( fQSecondBackground->GetBinContent(i,j,k) < threshold) return;
1793 }
1794 }
1795
1796 }
1797
1798
1799 Info("SwapGeneratingHistograms","*******************************************");
1800 Info("SwapGeneratingHistograms","*******************************************");
1801 Info("SwapGeneratingHistograms","*******************************************");
1802 Info("SwapGeneratingHistograms","**** SWAPPING HISTOGRAMS ****");
1803 Info("SwapGeneratingHistograms","*******************************************");
1804 Info("SwapGeneratingHistograms","*******************************************");
1805 Info("SwapGeneratingHistograms","*******************************************");
1806
1807
1808 TH3D* h = fQSignal;
1809 fQSignal = fQSecondSignal;
1810 fQSecondSignal = h;
1811 fQSecondSignal->Reset();
1812 fQSecondSignal->SetDirectory(0x0);
1813
1814 h = fQBackground;
1815 fQBackground = fQSecondBackground;
1816 fQSecondBackground = h;
1817 fQSecondBackground->Reset();
1818 fQSecondBackground->SetDirectory(0x0);
1819
1820 fSwapped = kTRUE;
1821
1822}