]> 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
An additional label
[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/***********************************************************/
83cb58d0 90using std::cout;
91using std::endl;
92using std::ofstream;
93using std::ios;
695dc56f 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),
5fdc214b 104 fQRange(0.06),
695dc56f 105 fQNBins(60),
106 fGenerator(0x0),
1c56e311 107 fStackBuffer(0),
695dc56f 108 fBufferSize(5),
5fdc214b 109 fNBinsToScale(Int_t(fQNBins*0.1)),
f5b253b3 110 fDebug(0),
1c56e311 111 fSignalShapeCreated(0),
5fdc214b 112 fMaxIterations(1000),
695dc56f 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),
5fdc214b 120 fMinFill(1000),
1c56e311 121 fSwapped(0),
5fdc214b 122 fLogFile(0x0)
695dc56f 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),
5fdc214b 136 fQRange(0.06),
695dc56f 137 fQNBins(60),
138 fGenerator(0x0),
1c56e311 139 fStackBuffer(0),
695dc56f 140 fBufferSize(5),
5fdc214b 141 fNBinsToScale(Int_t(fQNBins*0.1)),
f5b253b3 142 fDebug(0),
695dc56f 143 fSignalShapeCreated(kFALSE),
5fdc214b 144 fMaxIterations(1000),
695dc56f 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),
5fdc214b 152 fMinFill(1000),
1c56e311 153 fSwapped(0),
5fdc214b 154 fLogFile(0x0)
695dc56f 155{
156//default constructor
157}
1c56e311 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}
695dc56f 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;
5fdc214b 202 delete fLogFile;
695dc56f 203}
5fdc214b 204/***********************************************************/
695dc56f 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]);
5fdc214b 241 static const Double_t kDegToRadCF = 180./TMath::Pi();
695dc56f 242 fGenerator->SetMomentumRange(fPtMin,fPtMax);
5fdc214b 243 fGenerator->SetPhiRange(kDegToRadCF*fPhiMin,kDegToRadCF*fPhiMax);
695dc56f 244 fGenerator->SetYRange(fYMin,fYMax);
5fdc214b 245 fGenerator->SetThetaRange(kDegToRadCF*fThetaMin,kDegToRadCF*fThetaMax);
695dc56f 246 fGenerator->Init();
247
248 }
249
5fdc214b 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
695dc56f 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
695dc56f 261 fQSignal->Sumw2();
262 fQBackground->Sumw2();
263 fQSecondSignal->Sumw2();
264 fQSecondBackground->Sumw2();
5fdc214b 265
266 fLogFile = new ofstream("BadEvent",ios::out);
267
695dc56f 268}
269/***********************************************************/
270
271void AliGenHBTosl::Generate()
272{
273 //the main method
5fdc214b 274
275 ofstream& logfile = *fLogFile;
276 logfile<<"Generate"<<"Attempts to generate "<<fNpart<<" particles.";
277
278
695dc56f 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
5fdc214b 349 {
350 fQCoarseBackground = new TH3D("fQCoarseBackground","",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
351 fQCoarseBackground->Sumw2();
352
695dc56f 353 FillCoarse(); //create coarse background - just to know the spectrum
354 }
5fdc214b 355
356 fQCoarseSignal = new TH3D("fQCoarseSignal","fQCoarseSignal",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
357 fQCoarseSignal->Sumw2();
695dc56f 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();
ea79897e 380 Int_t ntr = 0;
695dc56f 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);
5fdc214b 390 work->SetDirectory(0x0);
391 work->Sumw2();
695dc56f 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 }
5fdc214b 430
431 char msg[1000];
432 logfile<<endl;
ddf405a9 433 snprintf(msg,1000, "\n");
695dc56f 434 Int_t middlebin = fQNBins/2;
5fdc214b 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);
ddf405a9 441 snprintf(msg,1000, "% 6.5f ",GetQOutQSideQLongCorrTheorValue(tx,ty,tz));
5fdc214b 442 logfile<<msg;
443 }
444 logfile<<endl;
445
695dc56f 446 for (Int_t k = middlebin-5; k < middlebin+5; k++)
447 {
ddf405a9 448 snprintf(msg,1000, "% 6.5f ",work->GetBinContent(30,30,k));
5fdc214b 449 logfile<<msg;
695dc56f 450 }
5fdc214b 451 logfile<<endl;
695dc56f 452
453 for (Int_t k = middlebin-5; k < middlebin+5; k++)
454 {
ddf405a9 455 snprintf(msg, 1000, "% 6.5f ",chiarray[30][30][k]);
5fdc214b 456 logfile<<msg;
695dc56f 457 }
5fdc214b 458 logfile<<endl;
695dc56f 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 }
5fdc214b 522 Double_t qlongc = work->GetZaxis()->GetBinCenter(zmax);
523 Double_t qsidec = work->GetYaxis()->GetBinCenter(ymax);
524 Double_t qoutc = work->GetXaxis()->GetBinCenter(xmax);
695dc56f 525
5fdc214b 526
ddf405a9 527 snprintf(msg,1000, "Generate Fill bin chi2(%d,%d,%d)=%f",xmax,ymax,zmax,chiarray[xmax][ymax][zmax]);
5fdc214b 528 logfile<<msg;
695dc56f 529
5fdc214b 530 qout = gRandom->Uniform(qoutc -binwdh, qoutc +binwdh);
531 qside = gRandom->Uniform(qsidec-binwdh, qsidec+binwdh);
532 qlong = gRandom->Uniform(qlongc-binwdh, qlongc+binwdh);
695dc56f 533
534 TParticle* first = 0;
ea79897e 535 Int_t jj = 0;
536
537 while (jj < genstack->GetNtrack())
695dc56f 538 {
ea79897e 539 TParticle* tmpp = genstack->Particle(jj++);
695dc56f 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
5fdc214b 584 Double_t y = GetQOutQSideQLongCorrTheorValue(qoutc,qsidec,qlongc);
695dc56f 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
5fdc214b 594
595 delete work;
596
695dc56f 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
e087ab36 655 AliDebug(1,Form("Randomized qinv = %f, obtained = %f",qinv,GetQInv(first,second)));
695dc56f 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////////////
e087ab36 779 if ( AliDebugLevel() > 3 )
695dc56f 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
5fdc214b 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
695dc56f 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
5fdc214b 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;
695dc56f 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
695dc56f 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
5fdc214b 892 logfile<<"StartSignal\n";
893 logfile<<" Row 1 Theory, 2 current value, 3 Chi2 \n";
695dc56f 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];
ddf405a9 905 snprintf(buff,50, "QTWorkPass2.%3d.root",counter);
695dc56f 906 TFile* file = TFile::Open(buff,"update");
907 work->Write();
908 work->SetDirectory(0x0);
909 file->Close();
910 delete file;
911
ddf405a9 912 snprintf(buff,50, "QTBackgroundPass2.%3d.root",counter);
695dc56f 913 file = TFile::Open(buff,"update");
914 fQBackground->Write();
915 fQBackground->SetDirectory(0x0);
916 file->Close();
917 delete file;
918
ddf405a9 919 snprintf(buff,50, "QTSignalPass2.%3d.root",counter);
695dc56f 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
5fdc214b 958 char msg[1000];
695dc56f 959
5fdc214b 960 printf("\n");
961
962 for (Int_t k = 25; k < 36; k++)
695dc56f 963 {
5fdc214b 964 Double_t tx = work->GetXaxis()->GetBinCenter(30);
965 Double_t ty = work->GetYaxis()->GetBinCenter(30);
695dc56f 966 Double_t tz = work->GetZaxis()->GetBinCenter(k);
ddf405a9 967 snprintf(msg,1000, "% 6.5f ",GetQOutQSideQLongCorrTheorValue(tx,ty,tz));
5fdc214b 968 logfile<<msg;
695dc56f 969 }
5fdc214b 970 logfile<<endl;
695dc56f 971
5fdc214b 972 for (Int_t k = 25; k < 36; k++)
695dc56f 973 {
ddf405a9 974 snprintf(msg, 1000, "%6.5f ",work->GetBinContent(30,30,k));
5fdc214b 975 logfile<<msg;
695dc56f 976 }
5fdc214b 977 logfile<<endl;
695dc56f 978
5fdc214b 979 for (Int_t k = 25; k < 36; k++)
695dc56f 980 {
ddf405a9 981 snprintf(msg,1000, "% 6.5f ",chiarray[30][30][k]);
5fdc214b 982 logfile<<msg;
695dc56f 983 }
5fdc214b 984 logfile<<endl;
695dc56f 985
986 chisqrchange = TMath::Abs(chisqrnew - chisqrold)/chisqrnew;
987 chisqrold = chisqrnew;
988
5fdc214b 989 chisqrPerDF = chisqrnew/kNDF;
695dc56f 990
991 Info("StartSignal","Iteration %d Chi-sq change %f%%",niterations,chisqrchange*100.0);
5fdc214b 992 Info("StartSignal","ChiSq = %f, NDF = %f, ChiSq/NDF = %f",chisqrnew, kNDF, chisqrPerDF );
695dc56f 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
5fdc214b 1162
1163
1164// if ( (chisqrPerDF < 2.0) && (fSwapped == kFALSE) )
1165// {
1166// SwapGeneratingHistograms();
1167// }
695dc56f 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 {
ea79897e 1299 for (j = 1; j<=fQNBins; j++)
695dc56f 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
e087ab36 1343 //if (AliDebugLevel())
695dc56f 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
5fdc214b 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)
695dc56f 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;
60e55aee 1572 Double_t gamma = 1.0/TMath::Sqrt((1.-beta)*(1.+beta));
695dc56f 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
e087ab36 1588 AliDebug(1,"Entered");
695dc56f 1589 if(!num)
1590 {
e087ab36 1591 AliError("No numerator");
695dc56f 1592 return 0.0;
1593 }
1594 if(!den)
1595 {
e087ab36 1596 AliError("No denominator");
695dc56f 1597 return 0.0;
1598 }
1599
1600 if(fNBinsToScale < 1)
1601 {
02ea70d8 1602 AliError("Number of bins for scaling is smaller than 1");
695dc56f 1603 return 0.0;
695dc56f 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 {
e087ab36 1612 AliError("Number of X bins for scaling is bigger thnan number of bins in histograms");
695dc56f 1613 return 0.0;
1614 }
1615
1616 Int_t nbinsY = num->GetNbinsX();
1617 if (fNBinsToScaleY > nbinsY)
1618 {
e087ab36 1619 AliError("Number of Y bins for scaling is bigger thnan number of bins in histograms");
695dc56f 1620 return 0.0;
1621 }
1622
1623 Int_t nbinsZ = num->GetNbinsZ();
1624 if (fNBinsToScaleZ > nbinsZ)
1625 {
e087ab36 1626 AliError("Number of Z bins for scaling is bigger thnan number of bins in histograms");
695dc56f 1627 return 0.0;
1628 }
1629
e087ab36 1630 AliDebug(1,"No errors detected");
695dc56f 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
e087ab36 1651 AliDebug(1,Form("numsum=%f densum=%f fNBinsToScaleX=%d fNBinsToScaleY=%d fNBinsToScaleZ=%d",
1652 numsum,densum,fNBinsToScaleX,fNBinsToScaleY,fNBinsToScaleZ));
695dc56f 1653
1654 if (numsum == 0) return 0.0;
1655 Double_t ret = densum/numsum;
1656
e087ab36 1657 AliDebug(1,Form("returning %f",ret));
695dc56f 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
5fdc214b 1717void AliGenHBTosl::SetTrack(TParticle* p, Int_t& ntr)
695dc56f 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
5fdc214b 1750void AliGenHBTosl::SetTrack(TParticle* p, Int_t& ntr, AliStack* stack) const
695dc56f 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;
e087ab36 1776 if (AliDebugLevel()>=1)
695dc56f 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
e087ab36 1793 if (AliDebugLevel()>5)
695dc56f 1794 {
1795 TVector3 test(firstPx,0.0,0.0);
1796 test.RotateY(rotAngleY);
1797 test.RotateZ(rotAngleZ);
e087ab36 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()));
695dc56f 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
e087ab36 1819 if (AliDebugLevel()>5)
5fdc214b 1820 {
e087ab36 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)));
5fdc214b 1824 }
1825
695dc56f 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
e087ab36 1834 AliInfo(Form("Calculated rot %f, modulus %f",xbis,TMath::Sqrt(x*x+y*y+z*z)));
695dc56f 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)
5fdc214b 1862 static const Double_t kFactorsqrd = 0.197*0.197;//squared conversion factor SI<->eV
695dc56f 1863
5fdc214b 1864 return 1.0 + 0.5*TMath::Exp(-qinv*qinv*fQRadius*fQRadius/kFactorsqrd);
695dc56f 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
5fdc214b 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);
695dc56f 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
02ea70d8 1880 //PH return kFALSE;
695dc56f 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}
1c56e311 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