]>
Commit | Line | Data |
---|---|---|
36b81802 | 1 | /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * |
2 | * See cxx source for full Copyright notice */ | |
3 | // Implementation of the interface for THBTprocessor | |
4 | // Author: Piotr Krzysztof Skowronski <Piotr.Skowronski@cern.ch> | |
5 | ////////////////////////////////////////////////////////////////////////////////// | |
6 | //Wrapper class for "hbt processor" after burner | |
7 | //The origibal code is written in fortran by Lanny Ray | |
8 | //and is put in the directory $ALICE_ROOT/HBTprocessor | |
9 | //Detailed description is on the top of the file hbt_event_processor.f | |
10 | // | |
11 | //This class can be used ONLY with AliGenCocktailAfterBurner wrapper generator | |
12 | //i.e. (in Config.C) | |
13 | // .... | |
14 | // AliGenCocktailAfterBurner = gener = new AliGenCocktailAfterBurner(); | |
15 | // gener->SetPhiRange(0, 360); //Set global parameters | |
16 | // gener->SetThetaRange(35., 145.); | |
17 | // AliGenHIJINGpara *hijing = new AliGenHIJINGpara(10000); //Add some generator | |
18 | // hijing->SetMomentumRange(0, 999); | |
19 | // gener->AddGenerator(hijing,"HIJING PARAMETRIZATION",1); | |
20 | // | |
21 | // AliGenHBTprocessor *hbtp = new AliGenHBTprocessor(); //create object | |
22 | // hbtp->SetRefControl(2); //set parameters | |
23 | // hbtp->SetSwitch1D(1); | |
24 | // hbtp->SetR0(6);//fm - radius | |
25 | // hbtp->SetLambda(0.7);//chaocity parameter | |
26 | // gener->AddAfterBurner(hbtp,"HBT PROCESSOR",1); //add to the main generator | |
27 | // | |
28 | // gener->Init(); | |
29 | // | |
30 | //CAUTIONS: | |
31 | // A) HBT PROCESSOR NEEDS MORE THAN ONE EVENT TO WORK | |
32 | // AS MORE AS IT BETTER WORKS | |
33 | // B) IT IS ABLE TO "ADD" CORRELATIONS ONLY UP TO TWO PARTICLE TYPES AT ONES | |
88cb7938 | 34 | // |
35 | // | |
36 | // Artificial particle dennsity enhancment feature | |
37 | // HBT Processor is unable to process correctly low multiplicity particles | |
38 | // even if the high statistics (more events) is supplied. | |
39 | // For that reason was implemented artificial multiplicity enhancement | |
40 | // - see also comments in HBT Processor source file (HBTP/hbt_event_processor.f) | |
41 | // or web page (http://alisoft.cern.ch/people/skowron/hbtprocessor/hbteventprocessor.html) | |
42 | // concerning trk_accep or SetTrackRejectionFactor method in this class | |
43 | // Fortran is cheated by masking several events a single one. Number is defined by fEventMerge | |
44 | // variable. In case it is equal to 1, there is no masking and the feature is off. | |
45 | // But, for example if it is 5, and we have 100 events, number of events passed to fortran is 100/5=20 | |
46 | // When fortran asks about multiplcity of event, let sey, 1 - multiplicity of events 5 to 9 is summed up | |
47 | // and returned. | |
48 | // | |
49 | // | |
36b81802 | 50 | ////////////////////////////////////////////////////////////////////////////////// |
51 | ||
52 | // 11.11.2001 Piotr Skowronski | |
53 | // Setting seed (date) in RNG in the constructor | |
54 | ||
55 | // 09.10.2001 Piotr Skowronski | |
56 | // | |
57 | // Theta - Eta cohernecy facilities added: | |
58 | // AliGenerator::SetThetaRange method overriden | |
59 | // Static methods added | |
60 | // EtaToTheta | |
61 | // ThetaToEta | |
62 | // DegreesToRadians | |
63 | // RadiansToDegrees | |
64 | // | |
65 | // Class description comments put on proper place | |
66 | ||
67 | // 27.09.2001 Piotr Skowronski | |
68 | // removing of redefinition of defaults velues | |
69 | // in method's implementation. | |
70 | // | |
71 | // | |
72 | ||
73 | #include "AliGenHBTprocessor.h" | |
88cb7938 | 74 | |
88cb7938 | 75 | #include <TParticle.h> |
76 | #include "THBTprocessor.h" | |
77 | ||
36b81802 | 78 | #include "AliStack.h" |
3e2e3ece | 79 | #include "AliMC.h" |
36b81802 | 80 | #include "AliGenCocktailAfterBurner.h" |
371660fe | 81 | #include "AliLog.h" |
36b81802 | 82 | |
83 | ||
84 | ||
85 | ClassImp(AliGenHBTprocessor) | |
86 | ||
88cb7938 | 87 | Int_t AliGenHBTprocessor::fgDebug = 0; |
b456eb2e | 88 | static TRandom* gAliRandom;//RNG to be used by the fortran code |
36b81802 | 89 | |
90 | AliGenCocktailAfterBurner* GetGenerator(); | |
91 | /*******************************************************************/ | |
36b81802 | 92 | |
b456eb2e | 93 | AliGenHBTprocessor::AliGenHBTprocessor(): |
94 | AliGenerator(), | |
95 | fHBTprocessor(0x0), | |
96 | fHbtPStatCodes(0x0), | |
b15b3ba3 | 97 | fNPDGCodes(0), |
98 | fTrackRejectionFactor(0), | |
99 | fReferenceControl(0), | |
100 | fPrintFull(0), | |
101 | fPrintSectorData(0), | |
102 | fNPidTypes(0), | |
103 | fNevents(0), | |
104 | fSwitch1d(0), | |
105 | fSwitch3d(0), | |
106 | fSwitchType(0), | |
107 | fSwitchCoherence(0), | |
108 | fSwitchCoulomb(0), | |
109 | fSwitchFermiBose(0), | |
110 | fEventLineCounter(0), | |
111 | fMaxit(0), | |
112 | fIrand(0), | |
113 | fLambda(0), | |
114 | fR1d(0), | |
115 | fRside(0), | |
116 | fRout(0), | |
117 | fRlong(0), | |
118 | fRperp(0), | |
119 | fRparallel(0), | |
120 | fR0(0), | |
121 | fQ0(0), | |
122 | fDeltap(0), | |
123 | fDelchi(0), | |
124 | fNPtBins(0), | |
125 | fNPhiBins(0), | |
126 | fNEtaBins(0), | |
127 | fN1dFine(0), | |
128 | fN1dCoarse(0), | |
129 | fN1dTotal(0), | |
130 | fN3dFine(0), | |
131 | fN3dCoarse(0), | |
132 | fN3dTotal(0), | |
133 | fN3dFineProject(0), | |
134 | fNPxBins(0), | |
135 | fNPyBins(0), | |
136 | fNPzBins(0), | |
137 | fNSectors(0), | |
138 | fPtBinSize(0), | |
139 | fPhiBinSize(0), | |
140 | fEtaBinSize(0), | |
141 | fEtaMin(0), | |
142 | fEtaMax(0), | |
143 | fBinsize1dFine(0), | |
144 | fBinsize1dCoarse(0), | |
145 | fQmid1d(0), | |
146 | fQmax1d(0), | |
147 | fBinsize3dFine(0), | |
148 | fBinsize3dCoarse(0), | |
149 | fQmid3d(0), | |
150 | fQmax3d(0), | |
151 | fPxMin(0), | |
152 | fPxMax(0), | |
153 | fDelpx(0), | |
154 | fPyMin(0), | |
155 | fPyMax(0), | |
156 | fDelpy(0), | |
157 | fPzMin(0), | |
158 | fPzMax(0), | |
159 | fDelpz(0), | |
160 | fEventMerge(1), | |
161 | fActiveStack(0) | |
36b81802 | 162 | { |
163 | // | |
164 | // Standard constructor | |
165 | // Sets default veues of all parameters | |
36b81802 | 166 | |
167 | SetName("AliGenHBTprocessor"); | |
168 | SetTitle("AliGenHBTprocessor"); | |
169 | ||
b456eb2e | 170 | gAliRandom = fRandom; |
36b81802 | 171 | fHBTprocessor = new THBTprocessor(); |
172 | ||
173 | fNPDGCodes = 0; | |
174 | DefineParticles(); | |
175 | ||
176 | SetTrackRejectionFactor(); | |
177 | SetRefControl(); | |
178 | SetPIDs(); | |
179 | SetNPIDtypes(); | |
180 | SetDeltap(); | |
181 | SetMaxIterations(); | |
182 | SetDelChi(); | |
183 | SetIRand(); | |
184 | SetLambda(); | |
185 | SetR1d() ; | |
186 | SetRSide(); | |
187 | SetROut() ; | |
188 | SetRLong() ; | |
189 | SetRPerp(); | |
190 | SetRParallel(); | |
191 | SetR0(); | |
192 | SetQ0(); | |
193 | SetSwitch1D(); | |
194 | SetSwitch3D(); | |
195 | SetSwitchType(); | |
196 | SetSwitchCoherence(); | |
197 | SetSwitchCoulomb(); | |
198 | SetSwitchFermiBose(); | |
199 | //SetMomentumRange(); | |
200 | SetPtRange(); | |
201 | SetPxRange(); | |
202 | SetPyRange(); | |
203 | SetPzRange(); | |
204 | SetPhiRange(); | |
205 | SetEtaRange(); | |
206 | SetNPtBins(); | |
207 | SetNPhiBins(); | |
208 | SetNEtaBins(); | |
209 | SetNPxBins(); | |
210 | SetNPyBins(); | |
211 | SetNPzBins(); | |
212 | SetNBins1DFineMesh(); | |
213 | SetBinSize1DFineMesh(); | |
214 | SetNBins1DCoarseMesh(); | |
215 | SetBinSize1DCoarseMesh(); | |
216 | SetNBins3DFineMesh(); | |
217 | SetBinSize3DFineMesh(); | |
218 | SetNBins3DCoarseMesh(); | |
219 | SetBinSize3DCoarseMesh(); | |
220 | SetNBins3DFineProjectMesh(); | |
221 | } | |
222 | ||
223 | /*******************************************************************/ | |
224 | ||
225 | ||
226 | /*******************************************************************/ | |
227 | ||
228 | AliGenHBTprocessor::~AliGenHBTprocessor() | |
229 | { | |
230 | //destructor | |
231 | CleanStatusCodes(); | |
232 | if (fHBTprocessor) delete fHBTprocessor; //delete generator | |
233 | ||
234 | } | |
235 | ||
236 | /*******************************************************************/ | |
237 | ||
238 | void AliGenHBTprocessor::InitStatusCodes() | |
239 | { | |
240 | //creates and inits status codes array to zero | |
241 | AliGenCocktailAfterBurner *cab = GetGenerator(); | |
242 | ||
243 | if(!cab) Fatal("InitStatusCodes()","Can not find AliGenCocktailAfterBurner generator"); | |
244 | ||
245 | Int_t nev = cab->GetNumberOfEvents(); | |
246 | ||
247 | fHbtPStatCodes = new Int_t* [nev]; | |
248 | for( Int_t i =0; i<nev;i++) | |
249 | { | |
250 | Int_t nprim = cab->GetStack(i)->GetNprimary(); | |
251 | fHbtPStatCodes[i] = new Int_t[nprim]; | |
252 | for (Int_t k =0 ;k<nprim;k++) | |
253 | fHbtPStatCodes[i][k] =0; | |
254 | ||
255 | } | |
256 | ||
257 | } | |
258 | /*******************************************************************/ | |
259 | ||
260 | void AliGenHBTprocessor::CleanStatusCodes() | |
261 | { | |
262 | //Cleans up status codes | |
263 | if (fHbtPStatCodes) | |
264 | { | |
265 | for (Int_t i =0; i<GetGenerator()->GetNumberOfEvents(); i++) | |
266 | delete [] fHbtPStatCodes[i]; | |
267 | delete fHbtPStatCodes; | |
268 | fHbtPStatCodes = 0; | |
269 | } | |
270 | ||
271 | } | |
272 | /*******************************************************************/ | |
273 | ||
274 | void AliGenHBTprocessor::Init() | |
275 | { | |
276 | //sets up parameters in generator | |
277 | ||
278 | THBTprocessor *thbtp = fHBTprocessor; | |
279 | ||
280 | ||
281 | thbtp->SetTrackRejectionFactor(fTrackRejectionFactor); | |
282 | thbtp->SetRefControl(fReferenceControl); | |
283 | ||
284 | if ((fPid[0] == fPid[1]) || (fPid[0] == 0) || (fPid[1] == 0)) | |
285 | { | |
286 | if (fPid[0] == 0) | |
287 | thbtp->SetPIDs(IdFromPDG(fPid[1]) ,0); | |
288 | else | |
289 | thbtp->SetPIDs(IdFromPDG(fPid[0]) ,0); | |
290 | thbtp->SetNPIDtypes(1); | |
291 | ||
292 | if (fSwitchType !=1) | |
293 | Warning("AliGenHBTprocessor::Init","\nThere is only one particle type set,\n\ | |
294 | and Switch_Type differnt then 1,\n which does not make sense.\n\ | |
295 | Setting it to 1.\n"); | |
296 | ||
297 | SetSwitchType(1); | |
298 | } | |
299 | else | |
300 | { | |
301 | thbtp->SetPIDs(IdFromPDG(fPid[0]) ,IdFromPDG(fPid[1])); | |
302 | SetNPIDtypes(2); | |
303 | thbtp->SetSwitchType(fSwitchType); | |
304 | } | |
305 | ||
306 | thbtp->SetMaxIterations(fMaxit); | |
307 | thbtp->SetDelChi(fDelchi); | |
308 | thbtp->SetIRand(fIrand); | |
309 | thbtp->SetLambda(fLambda); | |
310 | thbtp->SetR1d(fR1d); | |
311 | thbtp->SetRSide(fRside); | |
312 | thbtp->SetROut(fRout); | |
313 | thbtp->SetRLong(fRlong); | |
314 | thbtp->SetRPerp(fRperp); | |
315 | thbtp->SetRParallel(fRparallel); | |
316 | thbtp->SetR0(fR0); | |
317 | thbtp->SetQ0(fQ0); | |
318 | thbtp->SetSwitch1D(fSwitch1d); | |
319 | thbtp->SetSwitch3D(fSwitch3d); | |
320 | thbtp->SetSwitchType(fSwitchType); | |
321 | thbtp->SetSwitchCoherence(fSwitchCoherence); | |
322 | thbtp->SetSwitchCoulomb(fSwitchCoulomb); | |
323 | thbtp->SetSwitchFermiBose(fSwitchFermiBose); | |
324 | thbtp->SetPtRange(fPtMin,fPtMax); | |
325 | thbtp->SetPxRange(fPxMin,fPxMax); | |
326 | thbtp->SetPyRange(fPyMin,fPyMax); | |
327 | thbtp->SetPzRange(fPzMin,fPzMax); | |
88cb7938 | 328 | thbtp->SetPhiRange(fPhiMin*180./((Float_t)TMath::Pi())+180.0, //casting is because if fPhiMin = 180.0 then |
329 | fPhiMax*180./((Float_t)TMath::Pi())+180.0);//TMath::Pi() != TMath::Pi()*fPhiMin/180.0, | |
36b81802 | 330 | thbtp->SetEtaRange(fEtaMin,fEtaMax); |
331 | thbtp->SetNPtBins(fNPtBins); | |
332 | thbtp->SetNPhiBins(fNPhiBins); | |
333 | thbtp->SetNEtaBins(fNEtaBins); | |
334 | thbtp->SetNPxBins(fNPxBins); | |
335 | thbtp->SetNPyBins(fNPyBins); | |
336 | thbtp->SetNPzBins(fNPzBins); | |
337 | thbtp->SetNBins1DFineMesh(fN1dFine); | |
338 | thbtp->SetBinSize1DFineMesh(fBinsize1dFine); | |
339 | thbtp->SetNBins1DCoarseMesh(fN1dCoarse); | |
340 | thbtp->SetBinSize1DCoarseMesh(fBinsize1dCoarse); | |
341 | thbtp->SetNBins3DFineMesh(fN3dFine); | |
342 | thbtp->SetBinSize3DFineMesh(fBinsize3dFine); | |
343 | thbtp->SetNBins3DCoarseMesh(fN3dCoarse); | |
344 | thbtp->SetBinSize3DCoarseMesh(fBinsize3dCoarse); | |
345 | thbtp->SetNBins3DFineProjectMesh(fN3dFineProject); | |
88cb7938 | 346 | |
347 | thbtp->SetPrintFull(fPrintFull); | |
36b81802 | 348 | |
349 | } | |
350 | /*******************************************************************/ | |
351 | ||
352 | void AliGenHBTprocessor::Generate() | |
353 | { | |
354 | //starts processig | |
355 | AliGenCocktailAfterBurner* cab = GetGenerator(); | |
356 | if (cab == 0x0) | |
357 | { | |
358 | Fatal("Generate()","AliGenHBTprocessor needs AliGenCocktailAfterBurner to be main generator"); | |
359 | } | |
360 | if (cab->GetNumberOfEvents() <2) | |
361 | { | |
362 | Fatal("Generate()", | |
363 | "HBT Processor needs more than 1 event to work properly,\ | |
364 | but there is only %d set", cab->GetNumberOfEvents()); | |
365 | } | |
366 | ||
367 | ||
368 | fHBTprocessor->Initialize(); //reset all fields of common blocks | |
369 | //in case there are many HBT processors | |
370 | //run one after one (i.e. in AliCocktailAfterBurner) | |
371 | //between Init() called and Generate there might | |
372 | Init(); //be different instance running - be sure that we have our settings | |
373 | ||
374 | InitStatusCodes(); //Init status codes | |
375 | ||
376 | fHBTprocessor->GenerateEvent(); //Generates event | |
377 | ||
378 | CleanStatusCodes(); //Clean Status codes - thet are not needed anymore | |
379 | } | |
380 | ||
381 | /*******************************************************************/ | |
382 | ||
383 | ||
384 | /*******************************************************************/ | |
385 | void AliGenHBTprocessor::GetParticles(TClonesArray * particles) const | |
386 | { | |
387 | //practically dumm | |
388 | if (!particles) | |
389 | { | |
88cb7938 | 390 | Error("GetParticles","Particles has to be initialized"); |
36b81802 | 391 | return; |
392 | } | |
393 | fHBTprocessor->ImportParticles(particles); | |
394 | } | |
395 | ||
396 | /*******************************************************************/ | |
397 | ||
398 | Int_t AliGenHBTprocessor::GetHbtPStatusCode(Int_t part) const | |
399 | { | |
400 | //returns the status code of the given particle in the active event | |
401 | //see SetActiveEvent in the bottom of AliGenHBTprocessor.cxx | |
402 | //and in AliCocktailAfterBurner | |
88cb7938 | 403 | Int_t ev, idx; |
404 | GetTrackEventIndex(part,ev,idx); | |
405 | if ( (ev<0) || (idx<0) ) | |
406 | { | |
407 | Error("GetHbtPStatusCode","GetTrackEventIndex returned error"); | |
408 | return 0; | |
409 | } | |
0798b21e | 410 | return fHbtPStatCodes[ev][idx]; |
88cb7938 | 411 | |
36b81802 | 412 | } |
413 | ||
414 | /*******************************************************************/ | |
415 | void AliGenHBTprocessor::SetHbtPStatusCode(Int_t hbtstatcode, Int_t part) | |
416 | { | |
417 | //Sets the given status code to given particle number (part) in the active event | |
88cb7938 | 418 | Int_t ev, idx; |
419 | GetTrackEventIndex(part,ev,idx); | |
420 | if ( (ev<0) || (idx<0) ) | |
421 | { | |
422 | Error("GetHbtPStatusCode","GetTrackEventIndex returned error"); | |
423 | return; | |
424 | } | |
425 | else fHbtPStatCodes[ev][idx] = hbtstatcode; | |
36b81802 | 426 | } |
427 | ||
428 | /*******************************************************************/ | |
429 | ||
430 | void AliGenHBTprocessor::SetTrackRejectionFactor(Float_t trf) //def 1.0 | |
431 | { | |
432 | //Sets the Track Rejection Factor | |
433 | //variates in range 0.0 <-> 1.0 | |
434 | //Describes the factor of particles rejected from the output. | |
435 | //Used only in case of low muliplicity particles e.g. lambdas. | |
436 | //Processor generates addisional particles and builds the | |
437 | //correletions on such a statistics. | |
438 | //At the end these particels are left in the event according | |
439 | //to this factor: 1==all particles are left | |
440 | // 0==all are removed | |
441 | ||
442 | fTrackRejectionFactor=trf; | |
443 | fHBTprocessor->SetTrackRejectionFactor(trf); | |
444 | } | |
445 | /*******************************************************************/ | |
446 | ||
447 | void AliGenHBTprocessor::SetRefControl(Int_t rc) //default 2 | |
448 | { | |
449 | //Sets the Refernce Control Switch | |
450 | //switch wether read reference histograms from file =1 | |
451 | // compute from input events =2 - default | |
452 | fReferenceControl=rc; | |
453 | fHBTprocessor->SetRefControl(rc); | |
454 | } | |
455 | /*******************************************************************/ | |
456 | ||
457 | void AliGenHBTprocessor::SetPIDs(Int_t pid1,Int_t pid2) | |
458 | { | |
459 | //default pi+ pi- | |
460 | //Sets PDG Codes of particles to be processed, default \\Pi^{+} and \\Pi^{-} | |
461 | //This method accepts PDG codes which is | |
462 | //in opposite to THBBProcessor which accepts GEANT PID | |
463 | if ( (pid1 == 0) && (pid2 == 0) ) | |
464 | { | |
465 | Error("AliGenHBTprocessor::SetPIDs","Sensless Particle Types setting: 0 0, Ignoring\n"); | |
466 | } | |
467 | ||
468 | fPid[0]=pid1; | |
469 | fPid[1]=pid2; | |
470 | ||
471 | if(pid1 == pid2) | |
472 | { | |
473 | fHBTprocessor->SetPIDs(IdFromPDG(pid1) ,0); | |
474 | SetNPIDtypes(1); | |
475 | SetSwitchType(1); | |
476 | } | |
477 | else | |
478 | { | |
479 | fHBTprocessor->SetPIDs(IdFromPDG(pid1) ,IdFromPDG(pid2)); | |
480 | SetNPIDtypes(2); | |
481 | } | |
482 | } | |
483 | /*******************************************************************/ | |
484 | ||
485 | void AliGenHBTprocessor::SetNPIDtypes(Int_t npidt) | |
486 | { | |
487 | //Number ofparticle types to be processed - default 2 | |
488 | //see AliGenHBTprocessor::SetPIDs | |
489 | ||
490 | fNPidTypes = npidt; | |
491 | fHBTprocessor->SetNPIDtypes(npidt); | |
492 | } | |
493 | /*******************************************************************/ | |
494 | ||
495 | void AliGenHBTprocessor::SetDeltap(Float_t deltp) | |
496 | { | |
497 | //default = 0.1 GeV | |
498 | //maximum range for random momentum shifts in GeV/c; | |
499 | //px,py,pz independent; Default = 0.1 GeV/c. | |
500 | fDeltap=deltp; | |
501 | fHBTprocessor->SetDeltap(deltp); | |
502 | } | |
503 | /*******************************************************************/ | |
504 | ||
505 | void AliGenHBTprocessor::SetMaxIterations(Int_t maxiter) | |
506 | { | |
507 | //maximum # allowed iterations thru all the | |
508 | //tracks for each event. Default = 50. | |
509 | //If maxit=0, then calculate the correlations | |
510 | //for the input set of events without doing the | |
511 | //track adjustment procedure. | |
512 | ||
513 | fMaxit=maxiter; | |
514 | fHBTprocessor->SetMaxIterations(maxiter); | |
515 | } | |
516 | ||
517 | /*******************************************************************/ | |
518 | void AliGenHBTprocessor::SetDelChi(Float_t dc) | |
519 | { | |
520 | //min % change in total chi-square for which | |
521 | //the track adjustment iterations may stop, | |
522 | //Default = 0.1%. | |
523 | ||
524 | fDelchi=dc; | |
525 | fHBTprocessor->SetDelChi(dc); | |
526 | } | |
527 | /*******************************************************************/ | |
528 | ||
529 | void AliGenHBTprocessor::SetIRand(Int_t irnd) | |
530 | { | |
531 | //if fact dummy - used only for compatibility | |
532 | //we are using AliRoot built in RNG | |
533 | //dummy in fact since we are using aliroot build-in RNG | |
534 | //Sets renaodom number generator | |
535 | fIrand=irnd; | |
536 | fHBTprocessor->SetIRand(irnd); | |
537 | } | |
538 | /*******************************************************************/ | |
539 | ||
540 | void AliGenHBTprocessor::SetLambda(Float_t lam) | |
541 | { | |
542 | //default = 0.6 | |
543 | // Sets Chaoticity Parameter | |
544 | fLambda = lam; | |
545 | fHBTprocessor->SetLambda(lam); | |
546 | } | |
547 | /*******************************************************************/ | |
548 | void AliGenHBTprocessor::SetR1d(Float_t r) | |
549 | { | |
550 | //default 7.0 | |
551 | //Sets Spherical source model radius (fm) | |
552 | fR1d = r; | |
553 | fHBTprocessor->SetR1d(r); | |
554 | } | |
555 | /*******************************************************************/ | |
556 | void AliGenHBTprocessor::SetRSide(Float_t rs) | |
557 | { | |
558 | //default rs = 6.0 | |
559 | //Rside,Rout,Rlong = Non-spherical Bertsch-Pratt source model (fm) | |
560 | ||
561 | fRside = rs; | |
562 | fHBTprocessor->SetRSide(rs); | |
563 | } | |
564 | /*******************************************************************/ | |
565 | void AliGenHBTprocessor::SetROut(Float_t ro) | |
566 | { | |
567 | //default ro = 7.0 | |
568 | //Rside,Rout,Rlong = Non-spherical Bertsch-Pratt source model (fm) | |
569 | fRout = ro; | |
570 | fHBTprocessor->SetROut(ro); | |
571 | } | |
572 | /*******************************************************************/ | |
573 | void AliGenHBTprocessor::SetRLong(Float_t rl) | |
574 | { | |
575 | //default rl = 4.0 | |
576 | //Rside,Rout,Rlong = Non-spherical Bertsch-Pratt source model (fm) | |
577 | fRlong = rl; | |
578 | fHBTprocessor->SetRLong(rl); | |
579 | } | |
580 | /*******************************************************************/ | |
581 | void AliGenHBTprocessor::SetRPerp(Float_t rp) | |
582 | { | |
583 | //default rp = 6.0 | |
584 | //Rperp,Rparallel,R0= Non-spherical Yano-Koonin-Podgoretski source model (fm). | |
585 | fRperp = rp; | |
586 | fHBTprocessor->SetRPerp(rp); | |
587 | } | |
588 | /*******************************************************************/ | |
589 | void AliGenHBTprocessor::SetRParallel(Float_t rprl) | |
590 | { | |
591 | //default rprl = 4.0 | |
592 | //Rperp,Rparallel,R0= Non-spherical Yano-Koonin-Podgoretski source model (fm). | |
593 | fRparallel = rprl; | |
594 | fHBTprocessor->SetRParallel(rprl); | |
595 | } | |
596 | /*******************************************************************/ | |
597 | void AliGenHBTprocessor::SetR0(Float_t r0) | |
598 | { | |
599 | //default r0 = 4.0 | |
600 | //Rperp,Rparallel,R0= Non-spherical Yano-Koonin-Podgoretski source model (fm). | |
601 | fR0 = r0; | |
602 | fHBTprocessor->SetR0(r0); | |
603 | } | |
604 | /*******************************************************************/ | |
605 | void AliGenHBTprocessor::SetQ0(Float_t q0) | |
606 | { | |
607 | //default q0 = 9.0 | |
608 | //Sets Q0 = NA35 Coulomb parameter for finite source size in (GeV/c) | |
609 | // if fSwitchCoulomb = 2 | |
610 | // = Spherical Coulomb source radius in (fm) | |
611 | // if switchCoulomb = 3, used to interpolate the | |
612 | // input Pratt/Cramer discrete Coulomb source | |
613 | // radii tables. | |
614 | fQ0 = q0; | |
615 | fHBTprocessor->SetQ0(q0); | |
616 | } | |
617 | ||
618 | /*******************************************************************/ | |
619 | void AliGenHBTprocessor::SetSwitch1D(Int_t s1d) | |
620 | { | |
621 | //default s1d = 3 | |
622 | // Sets fSwitch1d | |
623 | // = 0 to not compute the 1D two-body //orrelations. | |
624 | // = 1 to compute this using Q-invariant | |
625 | // = 2 to compute this using Q-total | |
626 | // = 3 to compute this using Q-3-ve//tor | |
627 | ||
628 | fSwitch1d = s1d; | |
629 | fHBTprocessor->SetSwitch1D(s1d); | |
630 | } | |
631 | /*******************************************************************/ | |
632 | void AliGenHBTprocessor::SetSwitch3D(Int_t s3d) | |
633 | { | |
634 | //default s3d = 0 | |
635 | // Sets fSwitch3d | |
636 | // = 0 to not compute the 3D two-body correlations. | |
637 | // = 1 to compute this using the side-out-long form | |
638 | // = 2 to compute this using the Yanno-Koonin-Pogoredskij form. | |
639 | ||
640 | fSwitch3d = s3d; | |
641 | fHBTprocessor->SetSwitch3D(s3d); | |
642 | } | |
643 | /*******************************************************************/ | |
644 | void AliGenHBTprocessor::SetSwitchType(Int_t st) | |
645 | { | |
646 | //default st = 3 | |
647 | // Sets switch_type = 1 to fit only the like pair correlations | |
648 | // = 2 to fit only the unlike pair correlations | |
649 | // = 3 to fit both the like and unlike pair correl. | |
650 | //See SetPIDs and Init | |
651 | //If only one particle type is set, unly==1 makes sens | |
652 | ||
653 | fSwitchType = st; | |
654 | fHBTprocessor->SetSwitchType(st); | |
655 | } | |
656 | /*******************************************************************/ | |
657 | void AliGenHBTprocessor::SetSwitchCoherence(Int_t sc) | |
658 | { | |
659 | // default sc = 0 | |
660 | // switchCoherence = 0 for purely incoherent source (but can have | |
661 | // lambda < 1.0) | |
662 | // = 1 for mixed incoherent and coherent source | |
663 | ||
664 | fSwitchCoherence = sc; | |
665 | fHBTprocessor->SetSwitchCoherence(sc); | |
666 | } | |
667 | /*******************************************************************/ | |
668 | void AliGenHBTprocessor::SetSwitchCoulomb(Int_t scol) | |
669 | { | |
670 | //default scol = 2 | |
671 | // switchCoulomb = 0 no Coulomb correction | |
672 | // = 1 Point source Gamow correction only | |
673 | // = 2 NA35 finite source size correction | |
674 | // = 3 Pratt/Cramer finite source size correction; | |
675 | // interpolated from input tables. | |
676 | fSwitchCoulomb =scol; | |
677 | fHBTprocessor->SetSwitchCoulomb(scol); | |
678 | } | |
679 | /*******************************************************************/ | |
680 | void AliGenHBTprocessor::SetSwitchFermiBose(Int_t sfb) | |
681 | { | |
682 | //default sfb = 1 | |
683 | // switchFermiBose = 1 Boson pairs | |
684 | // = -1 Fermion pairs | |
685 | ||
686 | fSwitchFermiBose = sfb; | |
687 | fHBTprocessor->SetSwitchFermiBose(sfb); | |
688 | } | |
689 | /*******************************************************************/ | |
690 | void AliGenHBTprocessor::SetPtRange(Float_t ptmin, Float_t ptmax) | |
691 | { | |
692 | // default ptmin = 0.1, ptmax = 0.98 | |
693 | //Sets Pt range (GeV) | |
694 | AliGenerator::SetPtRange(ptmin,ptmax); | |
695 | fHBTprocessor->SetPtRange(ptmin,ptmax); | |
696 | } | |
697 | ||
698 | /*******************************************************************/ | |
699 | void AliGenHBTprocessor::SetPxRange(Float_t pxmin, Float_t pxmax) | |
700 | { | |
701 | //default pxmin = -1.0, pxmax = 1.0 | |
702 | //Sets Px range | |
703 | fPxMin =pxmin; | |
704 | fPxMax =pxmax; | |
705 | fHBTprocessor->SetPxRange(pxmin,pxmax); | |
706 | } | |
707 | /*******************************************************************/ | |
708 | void AliGenHBTprocessor::SetPyRange(Float_t pymin, Float_t pymax) | |
709 | { | |
710 | //default pymin = -1.0, pymax = 1.0 | |
711 | //Sets Py range | |
712 | fPyMin =pymin; | |
713 | fPyMax =pymax; | |
714 | fHBTprocessor->SetPyRange(pymin,pymax); | |
715 | } | |
716 | /*******************************************************************/ | |
717 | void AliGenHBTprocessor::SetPzRange(Float_t pzmin, Float_t pzmax) | |
718 | { | |
719 | //default pzmin = -3.6, pzmax = 3.6 | |
720 | //Sets Py range | |
721 | fPzMin =pzmin; | |
722 | fPzMax =pzmax; | |
723 | fHBTprocessor->SetPzRange(pzmin,pzmax); | |
724 | } | |
371660fe | 725 | void AliGenHBTprocessor::SetMomentumRange(Float_t /*pmin*/, Float_t /*pmax*/) |
36b81802 | 726 | { |
727 | //default pmin=0, pmax=0 | |
728 | //Do not use this method! | |
729 | MayNotUse("AliGenHBTprocessor::SetMomentumRange Method is Dummy"); | |
730 | } | |
731 | ||
732 | /*******************************************************************/ | |
733 | void AliGenHBTprocessor::SetPhiRange(Float_t phimin, Float_t phimax) | |
734 | { | |
735 | // | |
736 | //Sets \\Phi range | |
737 | AliGenerator::SetPhiRange(phimin,phimax); | |
738 | ||
88cb7938 | 739 | fHBTprocessor->SetPhiRange(phimin+180.0,phimax+180.0); |
36b81802 | 740 | } |
741 | /*******************************************************************/ | |
742 | void AliGenHBTprocessor::SetEtaRange(Float_t etamin, Float_t etamax) | |
743 | { | |
744 | //default etamin = -1.5, etamax = 1.5 | |
745 | //Sets \\Eta range | |
746 | fEtaMin= etamin; | |
747 | fEtaMax =etamax; | |
748 | fHBTprocessor->SetEtaRange(etamin,etamax); | |
749 | ||
750 | //set the azimothal angle range in the AliGeneraor - | |
751 | //to keep coherency between azimuthal angle and pseudorapidity | |
752 | //DO NOT CALL this->SetThetaRange, because it calls this method (where we are) | |
753 | //which must cause INFINITE LOOP | |
754 | AliGenerator::SetThetaRange(RadiansToDegrees(EtaToTheta(fEtaMin)), | |
755 | RadiansToDegrees(EtaToTheta(fEtaMax))); | |
756 | ||
757 | } | |
758 | /*******************************************************************/ | |
759 | void AliGenHBTprocessor::SetThetaRange(Float_t thetamin, Float_t thetamax) | |
760 | { | |
761 | //default thetamin = 0, thetamax = 180 | |
762 | //Azimuthal angle, override AliGenerator method which uses widely (i.e. wrapper generators) | |
763 | //core fortran HBTProcessor uses Eta (pseudorapidity) | |
764 | //so these methods has to be synchronized | |
765 | ||
766 | AliGenerator::SetThetaRange(thetamin,thetamax); | |
767 | ||
768 | SetEtaRange( ThetaToEta(fThetaMin) , ThetaToEta(fThetaMax) ); | |
769 | ||
770 | } | |
771 | ||
772 | /*******************************************************************/ | |
773 | void AliGenHBTprocessor::SetNPtBins(Int_t nptbin) | |
774 | { | |
775 | //default nptbin = 50 | |
776 | //set number of Pt bins | |
777 | fNPtBins= nptbin; | |
778 | fHBTprocessor->SetNPtBins(nptbin); | |
779 | } | |
780 | /*******************************************************************/ | |
781 | void AliGenHBTprocessor::SetNPhiBins(Int_t nphibin) | |
782 | { | |
783 | //default nphibin = 50 | |
784 | //set number of Phi bins | |
785 | fNPhiBins=nphibin; | |
786 | fHBTprocessor->SetNPhiBins(nphibin); | |
787 | } | |
788 | /*******************************************************************/ | |
789 | void AliGenHBTprocessor::SetNEtaBins(Int_t netabin) | |
790 | { | |
791 | //default netabin = 50 | |
792 | //set number of Eta bins | |
793 | fNEtaBins = netabin; | |
794 | fHBTprocessor->SetNEtaBins(netabin); | |
795 | } | |
796 | /*******************************************************************/ | |
797 | void AliGenHBTprocessor::SetNPxBins(Int_t npxbin) | |
798 | { | |
799 | //default npxbin = 20 | |
800 | //set number of Px bins | |
801 | fNPxBins = npxbin; | |
802 | fHBTprocessor->SetNPxBins(npxbin); | |
803 | } | |
804 | /*******************************************************************/ | |
805 | void AliGenHBTprocessor::SetNPyBins(Int_t npybin) | |
806 | { | |
807 | //default npybin = 20 | |
808 | //set number of Py bins | |
809 | fNPyBins = npybin; | |
810 | fHBTprocessor->SetNPyBins(npybin); | |
811 | } | |
812 | /*******************************************************************/ | |
813 | void AliGenHBTprocessor::SetNPzBins(Int_t npzbin) | |
814 | { | |
815 | //default npzbin = 70 | |
816 | //set number of Pz bins | |
817 | fNPzBins = npzbin; | |
818 | fHBTprocessor->SetNPzBins(npzbin); | |
819 | } | |
820 | /*******************************************************************/ | |
821 | void AliGenHBTprocessor::SetNBins1DFineMesh(Int_t n) | |
822 | { | |
823 | //default n = 10 | |
824 | //Sets the number of bins in the 1D mesh | |
825 | fN1dFine =n; | |
826 | fHBTprocessor->SetNBins1DFineMesh(n); | |
827 | ||
828 | } | |
829 | /*******************************************************************/ | |
830 | void AliGenHBTprocessor::SetBinSize1DFineMesh(Float_t x) | |
831 | { | |
832 | //default x=0.01 | |
833 | //Sets the bin size in the 1D mesh | |
834 | fBinsize1dFine = x; | |
835 | fHBTprocessor->SetBinSize1DFineMesh(x); | |
836 | } | |
837 | /*******************************************************************/ | |
838 | ||
839 | void AliGenHBTprocessor::SetNBins1DCoarseMesh(Int_t n) | |
840 | { | |
841 | //default n =2 | |
842 | //Sets the number of bins in the coarse 1D mesh | |
843 | fN1dCoarse =n; | |
844 | fHBTprocessor->SetNBins1DCoarseMesh(n); | |
845 | } | |
846 | /*******************************************************************/ | |
847 | void AliGenHBTprocessor::SetBinSize1DCoarseMesh(Float_t x) | |
848 | { | |
849 | //default x=0.05 | |
850 | //Sets the bin size in the coarse 1D mesh | |
851 | fBinsize1dCoarse =x; | |
852 | fHBTprocessor->SetBinSize1DCoarseMesh(x); | |
853 | } | |
854 | /*******************************************************************/ | |
855 | ||
856 | void AliGenHBTprocessor::SetNBins3DFineMesh(Int_t n) | |
857 | { | |
858 | //default n = 8 | |
859 | //Sets the number of bins in the 3D mesh | |
860 | fN3dFine =n; | |
861 | fHBTprocessor->SetNBins3DFineMesh(n); | |
862 | } | |
863 | /*******************************************************************/ | |
864 | void AliGenHBTprocessor::SetBinSize3DFineMesh(Float_t x) | |
865 | { | |
866 | //default x=0.01 | |
867 | //Sets the bin size in the 3D mesh | |
868 | fBinsize3dFine =x; | |
869 | fHBTprocessor->SetBinSize3DFineMesh(x); | |
870 | } | |
871 | /*******************************************************************/ | |
872 | ||
873 | void AliGenHBTprocessor::SetNBins3DCoarseMesh(Int_t n) | |
874 | { | |
875 | //default n = 2 | |
876 | //Sets the number of bins in the coarse 3D mesh | |
877 | ||
878 | fN3dCoarse = n; | |
879 | fHBTprocessor->SetNBins3DCoarseMesh(n); | |
880 | } | |
881 | /*******************************************************************/ | |
882 | void AliGenHBTprocessor::SetBinSize3DCoarseMesh(Float_t x) | |
883 | { | |
884 | //default x=0.08 | |
885 | //Sets the bin size in the coarse 3D mesh | |
886 | fBinsize3dCoarse = x; | |
887 | fHBTprocessor->SetBinSize3DCoarseMesh(x); | |
888 | } | |
889 | /*******************************************************************/ | |
890 | ||
891 | void AliGenHBTprocessor::SetNBins3DFineProjectMesh(Int_t n ) | |
892 | { | |
893 | //default n =3 | |
894 | //Sets the number of bins in the fine project mesh | |
895 | fN3dFineProject = n; | |
896 | fHBTprocessor->SetNBins3DFineProjectMesh(n); | |
897 | } | |
88cb7938 | 898 | /*******************************************************************/ |
899 | void AliGenHBTprocessor::SetPrintFull(Int_t flag) | |
900 | { | |
b456eb2e | 901 | //sets the print full flag |
88cb7938 | 902 | fPrintFull = flag; |
903 | fHBTprocessor->SetPrintFull(flag); | |
904 | } | |
905 | ||
906 | ||
907 | /*******************************************************************/ | |
908 | ||
909 | Int_t AliGenHBTprocessor::GetNumberOfEvents() | |
910 | { | |
b456eb2e | 911 | //returns number of available events |
3e2e3ece | 912 | AliGenerator* g = gAlice->GetMCApp()->Generator(); |
88cb7938 | 913 | AliGenCocktailAfterBurner* cab = (g)?dynamic_cast<AliGenCocktailAfterBurner*>(g):0x0; |
914 | if (cab == 0x0) | |
915 | { | |
916 | Fatal("GetNumberOfEvents","Master Generator is not an AliGenCocktailAfterBurner"); | |
917 | return 0; | |
918 | } | |
919 | return (Int_t)TMath::Ceil(cab->GetNumberOfEvents()/((Float_t)fEventMerge)); | |
920 | } | |
921 | ||
36b81802 | 922 | /*******************************************************************/ |
923 | ||
88cb7938 | 924 | void AliGenHBTprocessor::SetActiveEventNumber(Int_t n) |
925 | { | |
b456eb2e | 926 | //sets the active event |
88cb7938 | 927 | fActiveStack = n*fEventMerge; |
371660fe | 928 | AliDebug(1,Form("Settimg active event %d passed %d",fActiveStack,n)); |
88cb7938 | 929 | } |
930 | /*******************************************************************/ | |
36b81802 | 931 | |
88cb7938 | 932 | Int_t AliGenHBTprocessor::GetNumberOfTracks() |
933 | { | |
b456eb2e | 934 | //returns number of tracks in active event |
3e2e3ece | 935 | AliGenerator* g = gAlice->GetMCApp()->Generator(); |
88cb7938 | 936 | AliGenCocktailAfterBurner* cab = (g)?dynamic_cast<AliGenCocktailAfterBurner*>(g):0x0; |
937 | if (cab == 0x0) | |
938 | { | |
939 | Fatal("GetNumberOfEvents","Master Generator is not an AliGenCocktailAfterBurner"); | |
940 | return 0; | |
941 | } | |
942 | Int_t n = 0; | |
943 | for (Int_t i = fActiveStack;i < fActiveStack+fEventMerge; i++) | |
944 | { | |
945 | if (i >= GetNumberOfEvents()) break; //protection not to overshoot nb of events | |
946 | AliStack* stack = cab->GetStack(i); | |
947 | if (stack == 0x0) | |
948 | Error("GetNumberOfTracks","There is no stack %d",i); | |
949 | ||
950 | n+=stack->GetNprimary(); | |
951 | } | |
952 | return n; | |
953 | } | |
954 | /*******************************************************************/ | |
955 | ||
956 | void AliGenHBTprocessor::SetNEventsToMerge(Int_t nev) | |
957 | { | |
b456eb2e | 958 | //sets number of events to merge |
88cb7938 | 959 | if (nev > 0 ) fEventMerge = nev; |
960 | } | |
961 | /*******************************************************************/ | |
962 | ||
963 | TParticle* AliGenHBTprocessor::GetTrack(Int_t n) | |
964 | { | |
965 | //returns track that hbtp thinks is n in active event | |
371660fe | 966 | AliDebug(5,Form("n = %d",n)); |
3e2e3ece | 967 | AliGenerator* g = gAlice->GetMCApp()->Generator(); |
88cb7938 | 968 | AliGenCocktailAfterBurner* cab = (g)?dynamic_cast<AliGenCocktailAfterBurner*>(g):0x0; |
969 | if (cab == 0x0) | |
970 | { | |
971 | Fatal("GetTrackEventIndex","Master Generator is not an AliGenCocktailAfterBurner"); | |
972 | return 0; | |
973 | } | |
974 | ||
975 | Int_t ev, idx; | |
976 | GetTrackEventIndex(n,ev,idx); | |
371660fe | 977 | AliDebug(5,Form("Event = %d Particle = %d",ev,idx)); |
88cb7938 | 978 | if ( (ev<0) || (idx<0) ) |
979 | { | |
980 | Error("GetTrack","GetTrackEventIndex returned error"); | |
981 | return 0x0; | |
982 | } | |
371660fe | 983 | AliDebug(5,Form("Number of Tracks in Event(%d) = %d",ev,cab->GetStack(ev)->GetNprimary())); |
0798b21e | 984 | return cab->GetStack(ev)->Particle(idx);//safe - in case stack does not exist |
88cb7938 | 985 | } |
36b81802 | 986 | /*******************************************************************/ |
987 | ||
88cb7938 | 988 | void AliGenHBTprocessor::GetTrackEventIndex(Int_t n, Int_t &evno, Int_t &index) const |
989 | { | |
990 | //returns event(stack) number and particle index | |
3e2e3ece | 991 | AliGenerator* g = gAlice->GetMCApp()->Generator(); |
88cb7938 | 992 | AliGenCocktailAfterBurner* cab = (g)?dynamic_cast<AliGenCocktailAfterBurner*>(g):0x0; |
993 | if (cab == 0x0) | |
994 | { | |
995 | Fatal("GetTrackEventIndex","Master Generator is not an AliGenCocktailAfterBurner"); | |
996 | return; | |
997 | } | |
998 | ||
999 | evno = -1; | |
1000 | index = -1; | |
1001 | for (Int_t i = fActiveStack;i < fActiveStack+fEventMerge; i++) | |
1002 | { | |
1003 | AliStack* stack = cab->GetStack(i); | |
1004 | if (stack == 0x0) | |
1005 | { | |
1006 | Error("GetTrackEventIndex","There is no stack %d",i); | |
1007 | return; | |
1008 | } | |
1009 | ||
1010 | Int_t ntracks = stack->GetNprimary(); | |
371660fe | 1011 | AliDebug(10,Form("Event %d has %d tracks. n = %d",i,ntracks,n)); |
88cb7938 | 1012 | |
1013 | if ( ntracks > n) | |
1014 | { | |
1015 | evno = i; | |
1016 | index = n; | |
1017 | return ; | |
1018 | } | |
1019 | else | |
1020 | { | |
1021 | n-=ntracks; | |
1022 | continue; | |
1023 | } | |
1024 | } | |
1025 | Error("GetTrackEventIndex","Could not find given track"); | |
1026 | } | |
1027 | ||
1028 | /*******************************************************************/ | |
1029 | /*******************************************************************/ | |
1030 | /*******************************************************************/ | |
36b81802 | 1031 | |
1032 | ||
1033 | ||
1034 | ||
1035 | ||
1036 | void AliGenHBTprocessor::DefineParticles() | |
1037 | { | |
1038 | // | |
1039 | // Load standard numbers for GEANT particles and PDG conversion | |
1040 | fNPDGCodes = 0; //this is done in the constructor - but in any case ... | |
1041 | ||
1042 | fPDGCode[fNPDGCodes++]=-99; // 0 = unused location | |
1043 | fPDGCode[fNPDGCodes++]=22; // 1 = photon | |
1044 | fPDGCode[fNPDGCodes++]=-11; // 2 = positron | |
1045 | fPDGCode[fNPDGCodes++]=11; // 3 = electron | |
1046 | fPDGCode[fNPDGCodes++]=12; // 4 = neutrino e | |
1047 | fPDGCode[fNPDGCodes++]=-13; // 5 = muon + | |
1048 | fPDGCode[fNPDGCodes++]=13; // 6 = muon - | |
1049 | fPDGCode[fNPDGCodes++]=111; // 7 = pi0 | |
1050 | fPDGCode[fNPDGCodes++]=211; // 8 = pi+ | |
1051 | fPDGCode[fNPDGCodes++]=-211; // 9 = pi- | |
1052 | fPDGCode[fNPDGCodes++]=130; // 10 = Kaon Long | |
1053 | fPDGCode[fNPDGCodes++]=321; // 11 = Kaon + | |
1054 | fPDGCode[fNPDGCodes++]=-321; // 12 = Kaon - | |
1055 | fPDGCode[fNPDGCodes++]=2112; // 13 = Neutron | |
1056 | fPDGCode[fNPDGCodes++]=2212; // 14 = Proton | |
1057 | fPDGCode[fNPDGCodes++]=-2212; // 15 = Anti Proton | |
1058 | fPDGCode[fNPDGCodes++]=310; // 16 = Kaon Short | |
1059 | fPDGCode[fNPDGCodes++]=221; // 17 = Eta | |
1060 | fPDGCode[fNPDGCodes++]=3122; // 18 = Lambda | |
1061 | fPDGCode[fNPDGCodes++]=3222; // 19 = Sigma + | |
1062 | fPDGCode[fNPDGCodes++]=3212; // 20 = Sigma 0 | |
1063 | fPDGCode[fNPDGCodes++]=3112; // 21 = Sigma - | |
1064 | fPDGCode[fNPDGCodes++]=3322; // 22 = Xi0 | |
1065 | fPDGCode[fNPDGCodes++]=3312; // 23 = Xi- | |
1066 | fPDGCode[fNPDGCodes++]=3334; // 24 = Omega- | |
1067 | fPDGCode[fNPDGCodes++]=-2112; // 25 = Anti Proton | |
1068 | fPDGCode[fNPDGCodes++]=-3122; // 26 = Anti Proton | |
1069 | fPDGCode[fNPDGCodes++]=-3222; // 27 = Anti Sigma - | |
1070 | fPDGCode[fNPDGCodes++]=-3212; // 28 = Anti Sigma 0 | |
1071 | fPDGCode[fNPDGCodes++]=-3112; // 29 = Anti Sigma 0 | |
1072 | fPDGCode[fNPDGCodes++]=-3322; // 30 = Anti Xi 0 | |
1073 | fPDGCode[fNPDGCodes++]=-3312; // 31 = Anti Xi + | |
1074 | fPDGCode[fNPDGCodes++]=-3334; // 32 = Anti Omega + | |
1075 | } | |
1076 | ||
1077 | /*******************************************************************/ | |
1078 | Int_t AliGenHBTprocessor::IdFromPDG(Int_t pdg) const | |
1079 | { | |
1080 | // | |
1081 | // Return Geant3 code from PDG and pseudo ENDF code | |
1082 | // | |
1083 | for(Int_t i=0;i<fNPDGCodes;++i) | |
1084 | if(pdg==fPDGCode[i]) return i; | |
1085 | return -1; | |
1086 | } | |
1087 | Int_t AliGenHBTprocessor::PDGFromId(Int_t id) const | |
1088 | { | |
1089 | // | |
1090 | // Return PDG code and pseudo ENDF code from Geant3 code | |
1091 | // | |
1092 | if(id>0 && id<fNPDGCodes) return fPDGCode[id]; | |
1093 | else return -1; | |
1094 | } | |
1095 | Double_t AliGenHBTprocessor::ThetaToEta(Double_t arg) | |
1096 | { | |
1097 | //converts etha(azimuthal angle) to Eta (pseudorapidity). Argument in radians | |
88cb7938 | 1098 | |
1099 | if(arg>= TMath::Pi()) return 708.0;//This number is the biggest wich not crashes exp(x) | |
1100 | if(arg<= 0.0) return -708.0;// | |
1101 | if(arg == TMath::Pi()/2.) return 0.0;// | |
36b81802 | 1102 | |
36b81802 | 1103 | if (arg > 0.0) |
1104 | { | |
88cb7938 | 1105 | return TMath::Log( TMath::Tan(arg/2.)) ; |
36b81802 | 1106 | } |
1107 | else | |
88cb7938 | 1108 | { |
1109 | return -TMath::Log( TMath::Tan(-arg/2.)) ; | |
36b81802 | 1110 | } |
1111 | } | |
1112 | ||
1113 | /*******************************************************************/ | |
1114 | /****** ROUTINES USED FOR COMMUNUCATION ********/ | |
1115 | /******************** WITH FORTRAN ********************/ | |
1116 | /*******************************************************************/ | |
1117 | ||
1118 | #ifndef WIN32 | |
1119 | # define hbtpran hbtpran_ | |
1120 | # define alihbtp_puttrack alihbtp_puttrack_ | |
1121 | # define alihbtp_gettrack alihbtp_gettrack_ | |
1122 | # define alihbtp_getnumberevents alihbtp_getnumberevents_ | |
1123 | # define alihbtp_getnumbertracks alihbtp_getnumbertracks_ | |
1124 | # define alihbtp_initialize alihbtp_initialize_ | |
1125 | # define alihbtp_setactiveeventnumber alihbtp_setactiveeventnumber_ | |
1126 | # define alihbtp_setparameters alihbtp_setparameters_ | |
1127 | # define type_ofCall | |
1128 | ||
1129 | #else | |
1130 | # define hbtpran HBTPRAN | |
1131 | # define alihbtp_puttrack ALIHBTP_PUTTRACK | |
1132 | # define alihbtp_gettrack ALIHBTP_GETTRACK | |
1133 | # define alihbtp_getnumberevents ALIHBTP_GETNUMBEREVENTS | |
1134 | # define alihbtp_getnumbertracks ALIHBTP_GETNUMBERTRACKS | |
1135 | # define alihbtp_initialize ALIHBTP_INITIALIZE | |
1136 | # define alihbtp_setactiveeventnumber ALIHBTP_SETACTIVEEVENTNUMBER | |
1137 | # define alihbtp_setparameters ALIHBTP_SETPARAMETERS | |
1138 | # define type_ofCall _stdcall | |
1139 | #endif | |
1140 | ||
36b81802 | 1141 | /*******************************************************************/ |
1142 | ||
1143 | AliGenCocktailAfterBurner* GetGenerator() | |
1144 | { | |
1145 | // This function has two tasks: | |
1146 | // Check if environment is OK (exist gAlice and generator) | |
1147 | // Returns pointer to genrator | |
1148 | //to be changed with TFolders | |
1149 | ||
1150 | if(!gAlice) | |
1151 | { | |
88cb7938 | 1152 | ::Error("AliGenHBTprocessor.cxx: GetGenerator()", |
1153 | "There is NO gALICE! Check what you are doing!"); | |
1154 | ::Fatal("AliGenHBTprocessor.cxx: GetGenerator()", | |
1155 | "Running HBT Processor without gAlice... Exiting \n"); | |
1156 | return 0x0;//pro forma | |
36b81802 | 1157 | } |
3e2e3ece | 1158 | AliGenerator * gen = gAlice->GetMCApp()->Generator(); |
36b81802 | 1159 | |
1160 | if (!gen) | |
1161 | { | |
88cb7938 | 1162 | ::Fatal("AliGenHBTprocessor.cxx: GetGenerator()", |
1163 | "There is no generator in gAlice, exiting\n"); | |
36b81802 | 1164 | return 0x0; |
1165 | } | |
1166 | ||
1167 | //we do not sure actual type of the genetator | |
1168 | //and simple casting is risky - we use ROOT machinery to do safe cast | |
1169 | ||
1170 | TClass* cabclass = AliGenCocktailAfterBurner::Class(); //get AliGenCocktailAfterBurner TClass | |
1171 | TClass* genclass = gen->IsA();//get TClass of the generator we got from galice | |
1172 | //use casting implemented in TClass | |
1173 | //cast gen to cabclass | |
1174 | AliGenCocktailAfterBurner* cab=(AliGenCocktailAfterBurner*)genclass->DynamicCast(cabclass,gen); | |
1175 | ||
1176 | if (cab == 0x0)//if generator that we got is not AliGenCocktailAfterBurner or its descendant we get null | |
1177 | { //then quit with error | |
88cb7938 | 1178 | ::Fatal("AliGenHBTprocessor.cxx: GetGenerator()", |
1179 | "The main Generator is not a AliGenCocktailAfterBurner, exiting\n"); | |
36b81802 | 1180 | return 0x0; |
1181 | } | |
36b81802 | 1182 | return cab; |
36b81802 | 1183 | } |
1184 | /*******************************************************************/ | |
1185 | ||
1186 | AliGenHBTprocessor* GetAliGenHBTprocessor() | |
1187 | { | |
1188 | //returns pointer to the current instance of AliGenHBTprocessor in | |
1189 | //AliGenCocktailAfterBurner (can be more than one) | |
1190 | // | |
1191 | AliGenCocktailAfterBurner* gen = GetGenerator(); | |
1192 | AliGenerator* g = gen->GetCurrentGenerator(); | |
1193 | if(g == 0x0) | |
1194 | { | |
88cb7938 | 1195 | ::Fatal("AliGenHBTprocessor.cxx: GetAliGenHBTprocessor()", |
36b81802 | 1196 | "Can not get the current generator. Exiting"); |
88cb7938 | 1197 | return 0x0;//pro forma |
36b81802 | 1198 | } |
1199 | ||
1200 | TClass* hbtpclass = AliGenHBTprocessor::Class(); //get AliGenCocktailAfterBurner TClass | |
1201 | TClass* gclass = g->IsA();//get TClass of the current generator we got from CAB | |
1202 | AliGenHBTprocessor* hbtp = (AliGenHBTprocessor*)gclass->DynamicCast(hbtpclass,g);//try to cast | |
1203 | if (hbtp == 0x0) | |
1204 | { | |
88cb7938 | 1205 | ::Fatal("AliGenHBTprocessor.cxx: GetAliGenHBTprocessor()", |
1206 | "\nCurrernt generator in AliGenCocktailAfterBurner is not a AliGenHBTprocessor, exiting\n"); | |
36b81802 | 1207 | return 0x0; |
1208 | } | |
1209 | return hbtp; | |
1210 | } | |
1211 | ||
1212 | /*******************************************************************/ | |
1213 | extern "C" void type_ofCall alihbtp_setparameters() | |
1214 | { | |
1215 | //dummy | |
1216 | } | |
1217 | ||
1218 | extern "C" void type_ofCall alihbtp_initialize() | |
1219 | { | |
1220 | //dummy | |
1221 | } | |
1222 | ||
1223 | /*******************************************************************/ | |
1224 | ||
1225 | extern "C" void type_ofCall alihbtp_getnumberevents(Int_t &nev) | |
1226 | { | |
371660fe | 1227 | //passes number of events to the fortran |
b456eb2e | 1228 | if(AliGenHBTprocessor::GetDebug()) |
371660fe | 1229 | AliInfoGeneral("alihbtp_getnumberevents",Form("(%d) ....",nev)); |
88cb7938 | 1230 | |
1231 | AliGenHBTprocessor* gen = GetAliGenHBTprocessor();//we dont check because it is done in function | |
1232 | nev = gen->GetNumberOfEvents(); | |
36b81802 | 1233 | |
b456eb2e | 1234 | if(AliGenHBTprocessor::GetDebug()>5) |
371660fe | 1235 | AliInfoGeneral("alihbtp_getnumberevents",Form("EXITED N Ev = %d",nev)); |
36b81802 | 1236 | } |
1237 | ||
1238 | /*******************************************************************/ | |
1239 | ||
1240 | extern "C" void type_ofCall alihbtp_setactiveeventnumber(Int_t & nev) | |
1241 | { | |
1242 | //sets active event in generator (AliGenCocktailAfterBurner) | |
1243 | ||
b456eb2e | 1244 | if(AliGenHBTprocessor::GetDebug()>5) |
88cb7938 | 1245 | ::Info("AliGenHBTpocessor.cxx: alihbtp_setactiveeventnumber","(%d)",nev); |
b456eb2e | 1246 | if(AliGenHBTprocessor::GetDebug()>0) |
88cb7938 | 1247 | ::Info("AliGenHBTpocessor.cxx: alihbtp_setactiveeventnumber","Asked for event %d",nev-1); |
1248 | ||
1249 | AliGenHBTprocessor* gen = GetAliGenHBTprocessor();//we dont check because it is done in function | |
36b81802 | 1250 | |
88cb7938 | 1251 | gen->SetActiveEventNumber(nev - 1); //fortran numerates events from 1 to N |
36b81802 | 1252 | |
b456eb2e | 1253 | if(AliGenHBTprocessor::GetDebug()>5) |
88cb7938 | 1254 | ::Info("AliGenHBTpocessor.cxx: alihbtp_setactiveeventnumber","EXITED returned %d",nev); |
36b81802 | 1255 | } |
1256 | /*******************************************************************/ | |
1257 | ||
1258 | extern "C" void type_ofCall alihbtp_getnumbertracks(Int_t &ntracks) | |
1259 | { | |
1260 | //passes number of particles in active event to the fortran | |
b456eb2e | 1261 | if(AliGenHBTprocessor::GetDebug()>5) |
88cb7938 | 1262 | ::Info("AliGenHBTpocessor.cxx: alihbtp_getnumbertracks","(%d)",ntracks); |
36b81802 | 1263 | |
88cb7938 | 1264 | AliGenHBTprocessor* gen = GetAliGenHBTprocessor();//we dont check because it is done in function |
36b81802 | 1265 | |
88cb7938 | 1266 | ntracks = gen->GetNumberOfTracks(); |
b456eb2e | 1267 | if(AliGenHBTprocessor::GetDebug()>5) |
88cb7938 | 1268 | ::Info("AliGenHBTpocessor.cxx: alihbtp_getnumbertracks","EXITED Ntracks = %d",ntracks); |
36b81802 | 1269 | } |
1270 | ||
1271 | /*******************************************************************/ | |
1272 | ||
1273 | extern "C" void type_ofCall | |
1274 | alihbtp_puttrack(Int_t & n,Int_t& flag, Float_t& px, | |
1275 | Float_t& py, Float_t& pz, Int_t& geantpid) | |
1276 | { | |
1277 | //sets new parameters (momenta) in track number n | |
1278 | // in the active event | |
1279 | // n - number of the track in active event | |
1280 | // flag - flag of the track | |
1281 | // px,py,pz - momenta | |
1282 | // geantpid - type of the particle - Geant Particle ID | |
1283 | ||
b456eb2e | 1284 | if(AliGenHBTprocessor::GetDebug()>5) |
88cb7938 | 1285 | ::Info("AliGenHBTpocessor.cxx: alihbtp_puttrack","(%d)",n); |
36b81802 | 1286 | |
88cb7938 | 1287 | AliGenHBTprocessor* gen = GetAliGenHBTprocessor();//we dont check because it is done in function |
36b81802 | 1288 | |
88cb7938 | 1289 | TParticle * track = gen->GetTrack(n-1); |
1290 | if (track == 0x0) | |
1291 | { | |
1292 | ::Error("AliGenHBTprocessor.cxx","Can not get track from AliGenHBTprocessor"); | |
1293 | return; | |
1294 | } | |
36b81802 | 1295 | |
1296 | //check to be deleted | |
88cb7938 | 1297 | if (geantpid != (gen->IdFromPDG( track->GetPdgCode() ))) |
36b81802 | 1298 | { |
88cb7938 | 1299 | ::Error("AliGenHBTprocessor.cxx: alihbtp_puttrack", |
1300 | "SOMETHING IS GOING BAD:\n GEANTPIDS ARE NOT THE SAME"); | |
36b81802 | 1301 | } |
1302 | ||
b456eb2e | 1303 | if(AliGenHBTprocessor::GetDebug()>0) |
36b81802 | 1304 | if (px != track->Px()) |
88cb7938 | 1305 | ::Info("AliGenHBTprocessor.cxx: alihbtp_puttrack", |
1306 | "Px diff. = %f", px - track->Px()); | |
36b81802 | 1307 | |
b456eb2e | 1308 | if(AliGenHBTprocessor::GetDebug()>3) |
88cb7938 | 1309 | ::Info("AliGenHBTprocessor.cxx: alihbtp_puttrack", |
1310 | "track->GetPdgCode() --> %d",track->GetPdgCode()); | |
36b81802 | 1311 | |
1312 | ||
1313 | ||
1314 | Float_t m = track->GetMass(); | |
1315 | Float_t e = TMath::Sqrt(m*m+px*px+py*py+pz*pz); | |
1316 | track->SetMomentum(px,py,pz,e); | |
1317 | ||
88cb7938 | 1318 | gen->SetHbtPStatusCode(flag,n-1); |
36b81802 | 1319 | |
b456eb2e | 1320 | if(AliGenHBTprocessor::GetDebug()>5) ::Info("AliGenHBTprocessor.cxx: alihbtp_puttrack","EXITED "); |
36b81802 | 1321 | } |
1322 | ||
1323 | /*******************************************************************/ | |
1324 | ||
1325 | extern "C" void type_ofCall | |
1326 | alihbtp_gettrack(Int_t & n,Int_t & flag, Float_t & px, | |
1327 | Float_t & py, Float_t & pz, Int_t & geantpid) | |
1328 | ||
1329 | { | |
1330 | //passes track parameters to the fortran | |
1331 | // n - number of the track in active event | |
1332 | // flag - flag of the track | |
1333 | // px,py,pz - momenta | |
1334 | // geantpid - type of the particle - Geant Particle ID | |
1335 | ||
b456eb2e | 1336 | if(AliGenHBTprocessor::GetDebug()>5) ::Info("AliGenHBTprocessor.cxx: alihbtp_gettrack","(%d)",n); |
36b81802 | 1337 | AliGenHBTprocessor* g = GetAliGenHBTprocessor(); |
88cb7938 | 1338 | |
1339 | TParticle * track = g->GetTrack(n-1); | |
36b81802 | 1340 | |
1341 | flag = g->GetHbtPStatusCode(n-1); | |
1342 | ||
1343 | px = (Float_t)track->Px(); | |
1344 | py = (Float_t)track->Py(); | |
1345 | pz = (Float_t)track->Pz(); | |
1346 | ||
1347 | geantpid = g->IdFromPDG( track->GetPdgCode() ); | |
1348 | ||
b456eb2e | 1349 | if(AliGenHBTprocessor::GetDebug()>5) ::Info("AliGenHBTprocessor.cxx: alihbtp_gettrack","EXITED"); |
36b81802 | 1350 | } |
1351 | ||
1352 | /*******************************************************************/ | |
1353 | extern "C" Float_t type_ofCall hbtpran(Int_t &) | |
1354 | { | |
1355 | //interface to the random number generator | |
b456eb2e | 1356 | return gAliRandom->Rndm(); |
36b81802 | 1357 | } |
1358 | ||
1359 | /*******************************************************************/ | |
1360 | ||
1361 | ||
1362 | /*******************************************************************/ |