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