]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EVGEN/AliGenHBTosl.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update
[u/mrichter/AliRoot.git] / EVGEN / AliGenHBTosl.cxx
CommitLineData
e087ab36 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
695dc56f 18#include "AliGenHBTosl.h"
e087ab36 19#include "AliLog.h"
20
695dc56f 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
695dc56f 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
5fdc214b 88
695dc56f 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),
5fdc214b 100 fQRange(0.06),
695dc56f 101 fQNBins(60),
102 fGenerator(0x0),
1c56e311 103 fStackBuffer(0),
695dc56f 104 fBufferSize(5),
5fdc214b 105 fNBinsToScale(Int_t(fQNBins*0.1)),
f5b253b3 106 fDebug(0),
1c56e311 107 fSignalShapeCreated(0),
5fdc214b 108 fMaxIterations(1000),
695dc56f 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),
5fdc214b 116 fMinFill(1000),
1c56e311 117 fSwapped(0),
5fdc214b 118 fLogFile(0x0)
695dc56f 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),
5fdc214b 132 fQRange(0.06),
695dc56f 133 fQNBins(60),
134 fGenerator(0x0),
1c56e311 135 fStackBuffer(0),
695dc56f 136 fBufferSize(5),
5fdc214b 137 fNBinsToScale(Int_t(fQNBins*0.1)),
f5b253b3 138 fDebug(0),
695dc56f 139 fSignalShapeCreated(kFALSE),
5fdc214b 140 fMaxIterations(1000),
695dc56f 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),
5fdc214b 148 fMinFill(1000),
1c56e311 149 fSwapped(0),
5fdc214b 150 fLogFile(0x0)
695dc56f 151{
152//default constructor
153}
1c56e311 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}
695dc56f 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;
5fdc214b 198 delete fLogFile;
695dc56f 199}
5fdc214b 200/***********************************************************/
695dc56f 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]);
5fdc214b 237 static const Double_t kDegToRadCF = 180./TMath::Pi();
695dc56f 238 fGenerator->SetMomentumRange(fPtMin,fPtMax);
5fdc214b 239 fGenerator->SetPhiRange(kDegToRadCF*fPhiMin,kDegToRadCF*fPhiMax);
695dc56f 240 fGenerator->SetYRange(fYMin,fYMax);
5fdc214b 241 fGenerator->SetThetaRange(kDegToRadCF*fThetaMin,kDegToRadCF*fThetaMax);
695dc56f 242 fGenerator->Init();
243
244 }
245
5fdc214b 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
695dc56f 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
695dc56f 257 fQSignal->Sumw2();
258 fQBackground->Sumw2();
259 fQSecondSignal->Sumw2();
260 fQSecondBackground->Sumw2();
5fdc214b 261
262 fLogFile = new ofstream("BadEvent",ios::out);
263
695dc56f 264}
265/***********************************************************/
266
267void AliGenHBTosl::Generate()
268{
269 //the main method
5fdc214b 270
271 ofstream& logfile = *fLogFile;
272 logfile<<"Generate"<<"Attempts to generate "<<fNpart<<" particles.";
273
274
695dc56f 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
5fdc214b 345 {
346 fQCoarseBackground = new TH3D("fQCoarseBackground","",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
347 fQCoarseBackground->Sumw2();
348
695dc56f 349 FillCoarse(); //create coarse background - just to know the spectrum
350 }
5fdc214b 351
352 fQCoarseSignal = new TH3D("fQCoarseSignal","fQCoarseSignal",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
353 fQCoarseSignal->Sumw2();
695dc56f 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();
ea79897e 376 Int_t ntr = 0;
695dc56f 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);
5fdc214b 386 work->SetDirectory(0x0);
387 work->Sumw2();
695dc56f 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 }
5fdc214b 426
427 char msg[1000];
428 logfile<<endl;
ddf405a9 429 snprintf(msg,1000, "\n");
695dc56f 430 Int_t middlebin = fQNBins/2;
5fdc214b 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);
ddf405a9 437 snprintf(msg,1000, "% 6.5f ",GetQOutQSideQLongCorrTheorValue(tx,ty,tz));
5fdc214b 438 logfile<<msg;
439 }
440 logfile<<endl;
441
695dc56f 442 for (Int_t k = middlebin-5; k < middlebin+5; k++)
443 {
ddf405a9 444 snprintf(msg,1000, "% 6.5f ",work->GetBinContent(30,30,k));
5fdc214b 445 logfile<<msg;
695dc56f 446 }
5fdc214b 447 logfile<<endl;
695dc56f 448
449 for (Int_t k = middlebin-5; k < middlebin+5; k++)
450 {
ddf405a9 451 snprintf(msg, 1000, "% 6.5f ",chiarray[30][30][k]);
5fdc214b 452 logfile<<msg;
695dc56f 453 }
5fdc214b 454 logfile<<endl;
695dc56f 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 }
5fdc214b 518 Double_t qlongc = work->GetZaxis()->GetBinCenter(zmax);
519 Double_t qsidec = work->GetYaxis()->GetBinCenter(ymax);
520 Double_t qoutc = work->GetXaxis()->GetBinCenter(xmax);
695dc56f 521
5fdc214b 522
ddf405a9 523 snprintf(msg,1000, "Generate Fill bin chi2(%d,%d,%d)=%f",xmax,ymax,zmax,chiarray[xmax][ymax][zmax]);
5fdc214b 524 logfile<<msg;
695dc56f 525
5fdc214b 526 qout = gRandom->Uniform(qoutc -binwdh, qoutc +binwdh);
527 qside = gRandom->Uniform(qsidec-binwdh, qsidec+binwdh);
528 qlong = gRandom->Uniform(qlongc-binwdh, qlongc+binwdh);
695dc56f 529
530 TParticle* first = 0;
ea79897e 531 Int_t jj = 0;
532
533 while (jj < genstack->GetNtrack())
695dc56f 534 {
ea79897e 535 TParticle* tmpp = genstack->Particle(jj++);
695dc56f 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
5fdc214b 580 Double_t y = GetQOutQSideQLongCorrTheorValue(qoutc,qsidec,qlongc);
695dc56f 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
5fdc214b 590
591 delete work;
592
695dc56f 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
e087ab36 651 AliDebug(1,Form("Randomized qinv = %f, obtained = %f",qinv,GetQInv(first,second)));
695dc56f 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////////////
e087ab36 775 if ( AliDebugLevel() > 3 )
695dc56f 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
5fdc214b 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
695dc56f 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
5fdc214b 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;
695dc56f 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
695dc56f 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
5fdc214b 888 logfile<<"StartSignal\n";
889 logfile<<" Row 1 Theory, 2 current value, 3 Chi2 \n";
695dc56f 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];
ddf405a9 901 snprintf(buff,50, "QTWorkPass2.%3d.root",counter);
695dc56f 902 TFile* file = TFile::Open(buff,"update");
903 work->Write();
904 work->SetDirectory(0x0);
905 file->Close();
906 delete file;
907
ddf405a9 908 snprintf(buff,50, "QTBackgroundPass2.%3d.root",counter);
695dc56f 909 file = TFile::Open(buff,"update");
910 fQBackground->Write();
911 fQBackground->SetDirectory(0x0);
912 file->Close();
913 delete file;
914
ddf405a9 915 snprintf(buff,50, "QTSignalPass2.%3d.root",counter);
695dc56f 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
5fdc214b 954 char msg[1000];
695dc56f 955
5fdc214b 956 printf("\n");
957
958 for (Int_t k = 25; k < 36; k++)
695dc56f 959 {
5fdc214b 960 Double_t tx = work->GetXaxis()->GetBinCenter(30);
961 Double_t ty = work->GetYaxis()->GetBinCenter(30);
695dc56f 962 Double_t tz = work->GetZaxis()->GetBinCenter(k);
ddf405a9 963 snprintf(msg,1000, "% 6.5f ",GetQOutQSideQLongCorrTheorValue(tx,ty,tz));
5fdc214b 964 logfile<<msg;
695dc56f 965 }
5fdc214b 966 logfile<<endl;
695dc56f 967
5fdc214b 968 for (Int_t k = 25; k < 36; k++)
695dc56f 969 {
ddf405a9 970 snprintf(msg, 1000, "%6.5f ",work->GetBinContent(30,30,k));
5fdc214b 971 logfile<<msg;
695dc56f 972 }
5fdc214b 973 logfile<<endl;
695dc56f 974
5fdc214b 975 for (Int_t k = 25; k < 36; k++)
695dc56f 976 {
ddf405a9 977 snprintf(msg,1000, "% 6.5f ",chiarray[30][30][k]);
5fdc214b 978 logfile<<msg;
695dc56f 979 }
5fdc214b 980 logfile<<endl;
695dc56f 981
982 chisqrchange = TMath::Abs(chisqrnew - chisqrold)/chisqrnew;
983 chisqrold = chisqrnew;
984
5fdc214b 985 chisqrPerDF = chisqrnew/kNDF;
695dc56f 986
987 Info("StartSignal","Iteration %d Chi-sq change %f%%",niterations,chisqrchange*100.0);
5fdc214b 988 Info("StartSignal","ChiSq = %f, NDF = %f, ChiSq/NDF = %f",chisqrnew, kNDF, chisqrPerDF );
695dc56f 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
5fdc214b 1158
1159
1160// if ( (chisqrPerDF < 2.0) && (fSwapped == kFALSE) )
1161// {
1162// SwapGeneratingHistograms();
1163// }
695dc56f 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 {
ea79897e 1295 for (j = 1; j<=fQNBins; j++)
695dc56f 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
e087ab36 1339 //if (AliDebugLevel())
695dc56f 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
5fdc214b 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)
695dc56f 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;
60e55aee 1568 Double_t gamma = 1.0/TMath::Sqrt((1.-beta)*(1.+beta));
695dc56f 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
e087ab36 1584 AliDebug(1,"Entered");
695dc56f 1585 if(!num)
1586 {
e087ab36 1587 AliError("No numerator");
695dc56f 1588 return 0.0;
1589 }
1590 if(!den)
1591 {
e087ab36 1592 AliError("No denominator");
695dc56f 1593 return 0.0;
1594 }
1595
1596 if(fNBinsToScale < 1)
1597 {
02ea70d8 1598 AliError("Number of bins for scaling is smaller than 1");
695dc56f 1599 return 0.0;
695dc56f 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 {
e087ab36 1608 AliError("Number of X bins for scaling is bigger thnan number of bins in histograms");
695dc56f 1609 return 0.0;
1610 }
1611
1612 Int_t nbinsY = num->GetNbinsX();
1613 if (fNBinsToScaleY > nbinsY)
1614 {
e087ab36 1615 AliError("Number of Y bins for scaling is bigger thnan number of bins in histograms");
695dc56f 1616 return 0.0;
1617 }
1618
1619 Int_t nbinsZ = num->GetNbinsZ();
1620 if (fNBinsToScaleZ > nbinsZ)
1621 {
e087ab36 1622 AliError("Number of Z bins for scaling is bigger thnan number of bins in histograms");
695dc56f 1623 return 0.0;
1624 }
1625
e087ab36 1626 AliDebug(1,"No errors detected");
695dc56f 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
e087ab36 1647 AliDebug(1,Form("numsum=%f densum=%f fNBinsToScaleX=%d fNBinsToScaleY=%d fNBinsToScaleZ=%d",
1648 numsum,densum,fNBinsToScaleX,fNBinsToScaleY,fNBinsToScaleZ));
695dc56f 1649
1650 if (numsum == 0) return 0.0;
1651 Double_t ret = densum/numsum;
1652
e087ab36 1653 AliDebug(1,Form("returning %f",ret));
695dc56f 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
5fdc214b 1713void AliGenHBTosl::SetTrack(TParticle* p, Int_t& ntr)
695dc56f 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
5fdc214b 1746void AliGenHBTosl::SetTrack(TParticle* p, Int_t& ntr, AliStack* stack) const
695dc56f 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;
e087ab36 1772 if (AliDebugLevel()>=1)
695dc56f 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
e087ab36 1789 if (AliDebugLevel()>5)
695dc56f 1790 {
1791 TVector3 test(firstPx,0.0,0.0);
1792 test.RotateY(rotAngleY);
1793 test.RotateZ(rotAngleZ);
e087ab36 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()));
695dc56f 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
e087ab36 1815 if (AliDebugLevel()>5)
5fdc214b 1816 {
e087ab36 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)));
5fdc214b 1820 }
1821
695dc56f 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
e087ab36 1830 AliInfo(Form("Calculated rot %f, modulus %f",xbis,TMath::Sqrt(x*x+y*y+z*z)));
695dc56f 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)
5fdc214b 1858 static const Double_t kFactorsqrd = 0.197*0.197;//squared conversion factor SI<->eV
695dc56f 1859
5fdc214b 1860 return 1.0 + 0.5*TMath::Exp(-qinv*qinv*fQRadius*fQRadius/kFactorsqrd);
695dc56f 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
5fdc214b 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);
695dc56f 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
02ea70d8 1876 //PH return kFALSE;
695dc56f 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}
1c56e311 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