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