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