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