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