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