Temporari implementation of LocaltoGlobal and GlobaltoLocal
[u/mrichter/AliRoot.git] / HBTAN / AliHBTAnalysis.cxx
CommitLineData
bfb09ece 1//_________________________________________________________
2////////////////////////////////////////////////////////////////////////////
3//
4// class AliHBTAnalysis
5//
6// Central Object Of HBTAnalyser:
7// This class performs main looping within HBT Analysis
8// User must plug a reader of Type AliHBTReader
9// User plugs in coorelation and monitor functions
10// as well as monitor functions
11//
12// HBT Analysis Tool, which is integral part of AliRoot,
13// ALICE Off-Line framework:
14//
15// Piotr.Skowronski@cern.ch
16// more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
17//
18////////////////////////////////////////////////////////////////////////////
19//_________________________________________________________
20
1b446896 21#include "AliHBTAnalysis.h"
1b446896 22#include "AliHBTRun.h"
dc2c3f36 23#include "AliHBTEvent.h"
1b446896 24#include "AliHBTReader.h"
25#include "AliHBTParticle.h"
26#include "AliHBTParticleCut.h"
27#include "AliHBTPair.h"
28#include "AliHBTPairCut.h"
29#include "AliHBTFunction.h"
5c58441a 30#include "AliHBTMonitorFunction.h"
bed069a4 31#include "AliHBTEventBuffer.h"
9616170a 32#include "AliHBTPairCut.h"
81b7b887 33
9616170a 34#include <TSystem.h>
dc2c3f36 35#include <TBenchmark.h>
1b446896 36#include <TList.h>
37
9616170a 38#include <sys/types.h>
39#include <sys/stat.h>
40#include <unistd.h>
41#include <fcntl.h>
e4f2b1da 42
1b446896 43ClassImp(AliHBTAnalysis)
44
45const UInt_t AliHBTAnalysis::fgkFctnArraySize = 100;
81b7b887 46const UInt_t AliHBTAnalysis::fgkDefaultMixingInfo = 1000;
47const Int_t AliHBTAnalysis::fgkDefaultBufferSize = 5;
48
49AliHBTAnalysis::AliHBTAnalysis():
2dc7203b 50 fReader(0x0),
51 fNTrackFunctions(0),
52 fNParticleFunctions(0),
53 fNParticleAndTrackFunctions(0),
54 fNTrackMonitorFunctions(0),
55 fNParticleMonitorFunctions(0),
56 fNParticleAndTrackMonitorFunctions(0),
9616170a 57 fTrackFunctions ( new AliHBTOnePairFctn* [fgkFctnArraySize]),
58 fParticleFunctions ( new AliHBTOnePairFctn* [fgkFctnArraySize]),
59 fParticleAndTrackFunctions ( new AliHBTTwoPairFctn* [fgkFctnArraySize]),
60 fParticleMonitorFunctions ( new AliHBTMonOneParticleFctn* [fgkFctnArraySize]),
61 fTrackMonitorFunctions ( new AliHBTMonOneParticleFctn* [fgkFctnArraySize]),
62 fParticleAndTrackMonitorFunctions ( new AliHBTMonTwoParticleFctn* [fgkFctnArraySize]),
63 fPairCut(new AliHBTEmptyPairCut()),//empty cut - accepts all particles
64 fAntiMergingCut(0x0),
2dc7203b 65 fBufferSize(2),
66 fDisplayMixingInfo(fgkDefaultMixingInfo),
67 fIsOwner(kFALSE)
1b446896 68 {
9616170a 69 //default constructor
1b446896 70 }
491d1b5d 71/*************************************************************************************/
1b446896 72
81b7b887 73AliHBTAnalysis::AliHBTAnalysis(const AliHBTAnalysis& in):
bdec021f 74 TObject(in),
2dc7203b 75 fReader(0x0),
76 fNTrackFunctions(0),
77 fNParticleFunctions(0),
78 fNParticleAndTrackFunctions(0),
79 fNTrackMonitorFunctions(0),
80 fNParticleMonitorFunctions(0),
81 fNParticleAndTrackMonitorFunctions(0),
82 fTrackFunctions(0x0),
83 fParticleFunctions(0x0),
84 fParticleAndTrackFunctions(0x0),
85 fParticleMonitorFunctions(0x0),
86 fTrackMonitorFunctions(0x0),
87 fParticleAndTrackMonitorFunctions(0x0),
88 fPairCut(0x0),
9616170a 89 fAntiMergingCut(0x0),
2dc7203b 90 fBufferSize(fgkDefaultBufferSize),
91 fDisplayMixingInfo(fgkDefaultMixingInfo),
92 fIsOwner(kFALSE)
81b7b887 93 {
2dc7203b 94//copy constructor
81b7b887 95 Fatal("AliHBTAnalysis(const AliHBTAnalysis&)","Sensless");
96 }
97/*************************************************************************************/
bdec021f 98const AliHBTAnalysis& AliHBTAnalysis::operator=(const AliHBTAnalysis& /*right*/)
81b7b887 99 {
2dc7203b 100//operator =
81b7b887 101 Fatal("AliHBTAnalysis(const AliHBTAnalysis&)","Sensless");
102 return *this;
103 }
104/*************************************************************************************/
1b446896 105AliHBTAnalysis::~AliHBTAnalysis()
106 {
107 //destructor
108 //note that we do not delete functions itself
109 // they should be deleted by whom where created
110 //we only store pointers, and we use "new" only for pointers array
b70eb506 111 if (fIsOwner)
112 DeleteFunctions();
1b446896 113 delete [] fTrackFunctions;
114 delete [] fParticleFunctions;
115 delete [] fParticleAndTrackFunctions;
116
5c58441a 117 delete [] fParticleMonitorFunctions;
118 delete [] fTrackMonitorFunctions;
119 delete [] fParticleAndTrackMonitorFunctions;
120
1b446896 121 delete fPairCut; // always have an copy of an object - we create we dstroy
9616170a 122 delete fAntiMergingCut;
1b446896 123 }
124
125/*************************************************************************************/
e4f2b1da 126
81b7b887 127void AliHBTAnalysis::DeleteFunctions()
128{
e4f2b1da 129 //Deletes all functions added to analysis
81b7b887 130 UInt_t ii;
131 for(ii = 0;ii<fNParticleFunctions;ii++)
132 delete fParticleFunctions[ii];
e4f2b1da 133 fNParticleFunctions = 0;
81b7b887 134
135 for(ii = 0;ii<fNTrackFunctions;ii++)
136 delete fTrackFunctions[ii];
e4f2b1da 137 fNTrackFunctions = 0;
138
81b7b887 139 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
140 delete fParticleAndTrackFunctions[ii];
e4f2b1da 141 fNParticleAndTrackFunctions = 0;
81b7b887 142
143 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
144 delete fParticleMonitorFunctions[ii];
e4f2b1da 145 fNParticleMonitorFunctions = 0;
81b7b887 146
147 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
148 delete fTrackMonitorFunctions[ii];
e4f2b1da 149 fNTrackMonitorFunctions = 0;
81b7b887 150
151 for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
152 delete fParticleAndTrackMonitorFunctions[ii];
e4f2b1da 153 fNParticleAndTrackMonitorFunctions = 0;
bfb09ece 154
81b7b887 155}
e4f2b1da 156/*************************************************************************************/
157
158void AliHBTAnalysis::Init()
159{
2dc7203b 160//Initializeation method
161//calls Init for all functions
e4f2b1da 162 UInt_t ii;
163 for(ii = 0;ii<fNParticleFunctions;ii++)
164 fParticleFunctions[ii]->Init();
165
166 for(ii = 0;ii<fNTrackFunctions;ii++)
167 fTrackFunctions[ii]->Init();
168
169 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
170 fParticleAndTrackFunctions[ii]->Init();
171
172 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
173 fParticleMonitorFunctions[ii]->Init();
174
175 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
176 fTrackMonitorFunctions[ii]->Init();
177
178 for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
179 fParticleAndTrackMonitorFunctions[ii]->Init();
bfb09ece 180
181
e4f2b1da 182}
183/*************************************************************************************/
184
185void AliHBTAnalysis::ResetFunctions()
186{
187//In case fOwner is true, deletes all functions
188//in other case, just set number of analysis to 0
189 if (fIsOwner) DeleteFunctions();
190 else
191 {
192 fNParticleFunctions = 0;
193 fNTrackFunctions = 0;
194 fNParticleAndTrackFunctions = 0;
195 fNParticleMonitorFunctions = 0;
196 fNTrackMonitorFunctions = 0;
197 fNParticleAndTrackMonitorFunctions = 0;
198 }
199}
200/*************************************************************************************/
201
1b446896 202void AliHBTAnalysis::Process(Option_t* option)
203{
204 //default option = "TracksAndParticles"
205 //Main method of the HBT Analysis Package
206 //It triggers reading with the global cut (default is an empty cut)
207 //Than it checks options and data which are read
208 //if everything is OK, then it calls one of the looping methods
209 //depending on tfReaderhe option
210 //These methods differs on what thay are looping on
211 //
212 // METHOD OPTION
213 //--------------------------------------------------------------------
214 //ProcessTracksAndParticles - "TracksAndParticles"
215 // DEFAULT
216 // it loops over both, tracks(reconstructed) and particles(simulated)
217 // all function gethered in all 3 lists are called for each (double)pair
218 //
219 //ProcessTracks - "Tracks"
220 // it loops only on tracks(reconstructed),
221 // functions ONLY from fTrackFunctions list are called
222 //
223 //ProcessParticles - "Particles"
224 // it loops only on particles(simulated),
225 // functions ONLY from fParticleAndTrackFunctions list are called
226 //
227 //
228 if (!fReader)
229 {
01725374 230 Error("Process","The reader is not set");
1b446896 231 return;
232 }
233
1b446896 234 const char *oT = strstr(option,"Tracks");
235 const char *oP = strstr(option,"Particles");
236
dc2c3f36 237 Bool_t nonid = IsNonIdentAnalysis();
238
bed069a4 239 Init();
240
1b446896 241 if(oT && oP)
242 {
dc2c3f36 243 if (nonid) ProcessTracksAndParticlesNonIdentAnal();
244 else ProcessTracksAndParticles();
1b446896 245 return;
246 }
247
248 if(oT)
249 {
dc2c3f36 250 if (nonid) ProcessTracksNonIdentAnal();
251 else ProcessTracks();
1b446896 252 return;
253 }
254
255 if(oP)
256 {
dc2c3f36 257 if (nonid) ProcessParticlesNonIdentAnal();
258 else ProcessParticles();
1b446896 259 return;
260 }
261
262}
1b446896 263/*************************************************************************************/
491d1b5d 264
1b446896 265void AliHBTAnalysis::ProcessTracksAndParticles()
266{
9616170a 267//Makes analysis for both tracks and particles
268//mainly for resolution study and analysies with weighting algirithms
1b446896 269//In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
270//the loops are splited
271
9616170a 272// cuta on particles only
1b446896 273
274 AliHBTParticle * part1, * part2;
275 AliHBTParticle * track1, * track2;
276
277 AliHBTEvent * trackEvent, *partEvent;
bed069a4 278 AliHBTEvent * trackEvent1 = 0x0,*partEvent1 = 0x0;
1b446896 279 AliHBTEvent * trackEvent2,*partEvent2;
280
281// Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
282
bed069a4 283// Int_t nev = fReader->GetNumberOfTrackEvents();
1b446896 284 AliHBTPair * trackpair = new AliHBTPair();
285 AliHBTPair * partpair = new AliHBTPair();
491d1b5d 286
287 AliHBTPair * tmptrackpair;//temprary pointers to pairs
288 AliHBTPair * tmppartpair;
81b7b887 289
bed069a4 290 AliHBTEventBuffer partbuffer(fBufferSize);
291 AliHBTEventBuffer trackbuffer(fBufferSize);
292
81b7b887 293 register UInt_t ii;
7a2c8238 294
bed069a4 295 Bool_t nocorrfctns = (fNParticleFunctions == 0) && (fNTrackFunctions == 0) && (fNParticleAndTrackFunctions == 0);
296
297 fReader->Rewind();
298 Int_t i = -1;
299 while (fReader->Next() == kFALSE)
1b446896 300 {
bed069a4 301 i++;
302 partEvent= fReader->GetParticleEvent();
303 trackEvent = fReader->GetTrackEvent();
1b446896 304
bed069a4 305 if ( !partEvent || !trackEvent )
306 {
307 Error("ProcessTracksAndParticles","Can not get event");
308 return;
309 }
310
311 if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
312 {
313 Fatal("ProcessTracksAndParticles",
314 "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
315 i,partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
316 return;
317 }
1b446896 318
bed069a4 319 if(partEvent1 == 0x0)
1b446896 320 {
bed069a4 321 partEvent1 = new AliHBTEvent();
322 partEvent1->SetOwner(kTRUE);
81b7b887 323
bed069a4 324 trackEvent1 = new AliHBTEvent();
325 trackEvent1->SetOwner(kTRUE);
326 }
327 else
328 {
329 partEvent1->Reset();
330 trackEvent1->Reset();
331 }
332
333 for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
334 {
335 /***************************************/
336 /****** Looping same events ********/
337 /****** filling numerators ********/
338 /***************************************/
339 if ( (j%fDisplayMixingInfo) == 0)
81b7b887 340 Info("ProcessTracksAndParticles",
341 "Mixing particle %d from event %d with particles from event %d",j,i,i);
1b446896 342
343 part1= partEvent->GetParticle(j);
bed069a4 344 track1= trackEvent->GetParticle(j);
4d91c73a 345
346//PID imperfections ???
347// if( part1->GetPdgCode() != track1->GetPdgCode() )
348// {
349// Fatal("ProcessTracksAndParticles",
350// "Event %d: Particle %d: PID of simulated particle (%d) not the same of reconstructed track (%d)",
351// i,j, part1->GetPdgCode(),track1->GetPdgCode() );
352// return;
353// }
bed069a4 354
355 Bool_t firstcut = fPairCut->GetFirstPartCut()->Pass(part1);
9277d10a 356 if (fBufferSize != 0)
47d9a058 357 if ( (firstcut == kFALSE) || (fPairCut->GetSecondPartCut()->Pass(part1) == kFALSE) )
358 {
359 //accepted by any cut
360 // we have to copy because reader keeps only one event
361
362 partEvent1->AddParticle(new AliHBTParticle(*part1));
363 trackEvent1->AddParticle(new AliHBTParticle(*track1));
364 }
81b7b887 365
bed069a4 366 if (firstcut) continue;
367
81b7b887 368 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
e4f2b1da 369 fParticleMonitorFunctions[ii]->Process(part1);
81b7b887 370 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
e4f2b1da 371 fTrackMonitorFunctions[ii]->Process(track1);
81b7b887 372 for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
e4f2b1da 373 fParticleAndTrackMonitorFunctions[ii]->Process(track1,part1);
5c58441a 374
bed069a4 375 if (nocorrfctns) continue;
5c58441a 376
1b446896 377 for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
378 {
379 part2= partEvent->GetParticle(k);
d74e6a27 380 if (part1->GetUID() == part2->GetUID()) continue;
1b446896 381 partpair->SetParticles(part1,part2);
382
4d91c73a 383 track2= trackEvent->GetParticle(k);
1b446896 384 trackpair->SetParticles(track1,track2);
385
9616170a 386 if(fPairCut->Pass(partpair) ) //check pair cut
bed069a4 387 { //do not meets crietria of the pair cut, try with swapped pairs
9616170a 388 if( fPairCut->Pass(partpair->GetSwapedPair()) )
1b446896 389 continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
390 else
bed069a4 391 { //swaped pair meets all the criteria
491d1b5d 392 tmppartpair = partpair->GetSwapedPair();
393 tmptrackpair = trackpair->GetSwapedPair();
1b446896 394 }
395 }
491d1b5d 396 else
397 {//meets criteria of the pair cut
398 tmptrackpair = trackpair;
399 tmppartpair = partpair;
400 }
9616170a 401
1b446896 402 for(ii = 0;ii<fNParticleFunctions;ii++)
491d1b5d 403 fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
1b446896 404
405 for(ii = 0;ii<fNTrackFunctions;ii++)
491d1b5d 406 fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
1b446896 407
408 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
491d1b5d 409 fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair,tmppartpair);
bed069a4 410 //end of 2nd loop over particles from the same event
411 }//for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
412
413 /***************************************/
414 /***** Filling denominators *********/
415 /***************************************/
416 if (fBufferSize == 0) continue;
417
418 partbuffer.ResetIter();
419 trackbuffer.ResetIter();
420 Int_t m = 0;
421 while (( partEvent2 = partbuffer.Next() ))
422 {
423 trackEvent2 = trackbuffer.Next();
5c58441a 424
bed069a4 425 m++;
426 if ( (j%fDisplayMixingInfo) == 0)
427 Info("ProcessTracksAndParticles",
428 "Mixing particle %d from event %d with particles from event %d",j,i,i-m);
429
430 for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++) // ... on all particles
1b446896 431 {
1b446896 432 part2= partEvent2->GetParticle(l);
433 partpair->SetParticles(part1,part2);
434
435 track2= trackEvent2->GetParticle(l);
436 trackpair->SetParticles(track1,track2);
437
81b7b887 438 if( fPairCut->Pass(partpair) ) //check pair cut
1b446896 439 { //do not meets crietria of the
81b7b887 440 if( fPairCut->Pass(partpair->GetSwapedPair()) )
441 continue;
442 else
443 {
444 tmppartpair = partpair->GetSwapedPair();
445 tmptrackpair = trackpair->GetSwapedPair();
446 }
1b446896 447 }
491d1b5d 448 else
449 {//meets criteria of the pair cut
450 tmptrackpair = trackpair;
451 tmppartpair = partpair;
452 }
9616170a 453
454 //Anti merging cut is only on tracks makes the background
455 if (fAntiMergingCut)
456 if (fAntiMergingCut->Pass(trackpair)) continue;
457
1b446896 458 for(ii = 0;ii<fNParticleFunctions;ii++)
491d1b5d 459 fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
1b446896 460
461 for(ii = 0;ii<fNTrackFunctions;ii++)
491d1b5d 462 fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
1b446896 463
464 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
491d1b5d 465 fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair,tmppartpair);
1b446896 466 }//for(Int_t l = 0; l<N2;l++) // ... on all particles
bed069a4 467
468 }
469 //end of loop over particles from first event
470 }//for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
471 partEvent1 = partbuffer.Push(partEvent1);
472 trackEvent1 = trackbuffer.Push(trackEvent1);
473 //end of loop over events
474 }//while (fReader->Next() == kFALSE)
1b446896 475}
476/*************************************************************************************/
477
478void AliHBTAnalysis::ProcessTracks()
479{
2dc7203b 480//In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
1b446896 481//the loops are splited
482 AliHBTParticle * track1, * track2;
1b446896 483 AliHBTEvent * trackEvent;
bed069a4 484 AliHBTEvent * trackEvent1 = 0x0;
1b446896 485 AliHBTEvent * trackEvent2;
81b7b887 486
bed069a4 487 register UInt_t ii;
81b7b887 488
489 AliHBTPair * trackpair = new AliHBTPair();
490 AliHBTPair * tmptrackpair; //temporary pointer
1b446896 491
bed069a4 492 AliHBTEventBuffer trackbuffer(fBufferSize);
493
494 fReader->Rewind();
495 Int_t i = -1;
496 while (fReader->Next() == kFALSE)
1b446896 497 {
bed069a4 498 i++;
499 trackEvent = fReader->GetTrackEvent();
1b446896 500 if (!trackEvent) continue;
1b446896 501
bed069a4 502 if(trackEvent1 == 0x0)
503 {
504 trackEvent1 = new AliHBTEvent();
505 trackEvent1->SetOwner(kTRUE);
506 }
507 else
508 {
509 trackEvent1->Reset();
510 }
511
1b446896 512 for (Int_t j = 0; j<trackEvent->GetNumberOfParticles() ; j++)
513 {
bed069a4 514 /***************************************/
515 /****** Looping same events ********/
516 /****** filling numerators ********/
517 /***************************************/
81b7b887 518 if ( (j%fDisplayMixingInfo) == 0)
519 Info("ProcessTracks",
520 "Mixing particle %d from event %d with particles from event %d",j,i,i);
521
bed069a4 522 track1= trackEvent->GetParticle(j);
523 Bool_t firstcut = fPairCut->GetFirstPartCut()->Pass(track1);
524
9277d10a 525 if (fBufferSize != 0)
47d9a058 526 if ( (firstcut == kFALSE) || (fPairCut->GetSecondPartCut()->Pass(track1) == kFALSE) )
527 {
528 //accepted by any cut
529 // we have to copy because reader keeps only one event
530 trackEvent1->AddParticle(new AliHBTParticle(*track1));
531 }
bed069a4 532
533 if (firstcut) continue;
81b7b887 534
535 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
e4f2b1da 536 fTrackMonitorFunctions[ii]->Process(track1);
5c58441a 537
bed069a4 538 if ( fNTrackFunctions ==0 ) continue;
5c58441a 539
1b446896 540 for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
541 {
d74e6a27 542 track2= trackEvent->GetParticle(k);
543 if (track1->GetUID() == track2->GetUID()) continue;
544
5c58441a 545 trackpair->SetParticles(track1,track2);
546 if(fPairCut->Pass(trackpair)) //check pair cut
81b7b887 547 { //do not meets crietria of the
5c58441a 548 if( fPairCut->Pass(trackpair->GetSwapedPair()) ) continue;
549 else tmptrackpair = trackpair->GetSwapedPair();
81b7b887 550 }
551 else
552 {
553 tmptrackpair = trackpair;
554 }
555 for(ii = 0;ii<fNTrackFunctions;ii++)
556 fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
557 }
bed069a4 558 /***************************************/
559 /***** Filling denominators *********/
560 /***************************************/
561
562 if (fBufferSize == 0) continue;
563
564 trackbuffer.ResetIter();
565
40d40d72 566 Int_t m = 0;
bed069a4 567 while (( trackEvent2 = trackbuffer.Next() ))
81b7b887 568 {
40d40d72 569 m++;
570 if ( (j%fDisplayMixingInfo) == 0)
571 Info("ProcessTracks",
572 "Mixing track %d from event %d with tracks from event %d",j,i,i-m);
81b7b887 573 for(Int_t l = 0; l<trackEvent2->GetNumberOfParticles();l++) // ... on all particles
bed069a4 574 {
575
576 track2= trackEvent2->GetParticle(l);
577 trackpair->SetParticles(track1,track2);
578
579 if( fPairCut->Pass(trackpair) ) //check pair cut
580 { //do not meets crietria of the
581 if( fPairCut->Pass(trackpair->GetSwapedPair()) )
582 continue;
583 else
584 {
585 tmptrackpair = trackpair->GetSwapedPair();
586 }
587 }
588 else
589 {//meets criteria of the pair cut
81b7b887 590 tmptrackpair = trackpair;
bed069a4 591 }
592
9616170a 593 //Anti merging cut is only on tracks makes the background
594 if (fAntiMergingCut)
595 if (fAntiMergingCut->Pass(trackpair)) continue;
596
bed069a4 597 for(ii = 0;ii<fNTrackFunctions;ii++)
598 fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
599
81b7b887 600 }//for(Int_t l = 0; l<N2;l++) // ... on all particles
bed069a4 601 }
1b446896 602 }
bed069a4 603 trackEvent1 = trackbuffer.Push(trackEvent1);
604 }//while (fReader->Next() == kFALSE)
1b446896 605}
606
607/*************************************************************************************/
491d1b5d 608
1b446896 609void AliHBTAnalysis::ProcessParticles()
610{
81b7b887 611//In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
1b446896 612//the loops are splited
613 AliHBTParticle * part1, * part2;
1b446896 614 AliHBTEvent * partEvent;
bed069a4 615 AliHBTEvent * partEvent1 = 0x0;
1b446896 616 AliHBTEvent * partEvent2;
491d1b5d 617
bed069a4 618 register UInt_t ii;
619
491d1b5d 620 AliHBTPair * partpair = new AliHBTPair();
bed069a4 621 AliHBTPair * tmppartpair; //temporary pointer
1b446896 622
bed069a4 623 AliHBTEventBuffer partbuffer(fBufferSize);
1b446896 624
bed069a4 625 fReader->Rewind();
626 Int_t i = -1;
627 while (fReader->Next() == kFALSE)
1b446896 628 {
bed069a4 629 i++;
630 partEvent = fReader->GetParticleEvent();
1b446896 631 if (!partEvent) continue;
bed069a4 632
633 if(partEvent1 == 0x0)
634 {
635 partEvent1 = new AliHBTEvent();
636 partEvent1->SetOwner(kTRUE);
637 }
638 else
639 {
640 partEvent1->Reset();
641 }
1b446896 642
643 for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
644 {
bed069a4 645 /***************************************/
646 /****** Looping same events ********/
647 /****** filling numerators ********/
648 /***************************************/
81b7b887 649 if ( (j%fDisplayMixingInfo) == 0)
40d40d72 650 Info("ProcessParticles",
bed069a4 651 "Mixing particle %d from event %d with particles from event %d",j,i,i);
652
40d40d72 653 part1 = partEvent->GetParticle(j);
bed069a4 654 Bool_t firstcut = fPairCut->GetFirstPartCut()->Pass(part1);
655
9277d10a 656 if (fBufferSize != 0) //useless in case
47d9a058 657 if ( (firstcut == kFALSE) || (fPairCut->GetSecondPartCut()->Pass(part1) == kFALSE) )
658 {
659 //accepted by any cut
660 // we have to copy because reader keeps only one event
661 partEvent1->AddParticle(new AliHBTParticle(*part1));
662 }
1b446896 663
bed069a4 664 if (firstcut) continue;
5c58441a 665
bed069a4 666 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
667 fParticleMonitorFunctions[ii]->Process(part1);
1b446896 668
bed069a4 669 if ( fNParticleFunctions == 0 ) continue;
670
671 for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
672 {
673 part2= partEvent->GetParticle(k);
674 if (part1->GetUID() == part2->GetUID()) continue;
81b7b887 675
bed069a4 676 partpair->SetParticles(part1,part2);
677 if(fPairCut->Pass(partpair)) //check pair cut
678 { //do not meets crietria of the
679 if( fPairCut->Pass(partpair->GetSwapedPair()) ) continue;
680 else tmppartpair = partpair->GetSwapedPair();
7a2c8238 681 }
bed069a4 682 else
7a2c8238 683 {
bed069a4 684 tmppartpair = partpair;
7a2c8238 685 }
bed069a4 686 for(ii = 0;ii<fNParticleFunctions;ii++)
687 fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
688 }
689 /***************************************/
690 /***** Filling denominators *********/
691 /***************************************/
692
693 if (fBufferSize == 0) continue;
694
695 partbuffer.ResetIter();
696
40d40d72 697 Int_t m = 0;
bed069a4 698 while (( partEvent2 = partbuffer.Next() ))
699 {
40d40d72 700 m++;
701 if ( (j%fDisplayMixingInfo) == 0)
702 Info("ProcessParticles",
703 "Mixing particle %d from event %d with particles from event %d",j,i,i-m);
bed069a4 704 for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++) // ... on all particles
1b446896 705 {
bed069a4 706
1b446896 707 part2= partEvent2->GetParticle(l);
708 partpair->SetParticles(part1,part2);
bed069a4 709
710 if( fPairCut->Pass(partpair) ) //check pair cut
1b446896 711 { //do not meets crietria of the
bed069a4 712 if( fPairCut->Pass(partpair->GetSwapedPair()) )
713 continue;
714 else
715 {
716 tmppartpair = partpair->GetSwapedPair();
717 }
491d1b5d 718 }
719 else
bed069a4 720 {//meets criteria of the pair cut
721 tmppartpair = partpair;
722 }
723
1b446896 724 for(ii = 0;ii<fNParticleFunctions;ii++)
491d1b5d 725 fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
bed069a4 726
727 }//for(Int_t l = 0; l<N2;l++) // ... on all particles
728 }
1b446896 729 }
bed069a4 730 partEvent1 = partbuffer.Push(partEvent1);
731 }//while (fReader->Next() == kFALSE)
1b446896 732}
9616170a 733/*************************************************************************************/
734
735void AliHBTAnalysis::SetAntiMergingCut(Float_t x)
736{
737 if (x < 0)
738 {
739 Error("SetAntiMergingCut","Value less then 0");
740 return;
741 }
742 if (fAntiMergingCut == 0x0) fAntiMergingCut = new AliHBTAvSeparationCut(x);
743 else fAntiMergingCut->SetMinimum(x);
744}
1b446896 745/*************************************************************************************/
746
491d1b5d 747void AliHBTAnalysis::WriteFunctions()
1b446896 748{
81b7b887 749//Calls Write for all defined functions in analysis
750//== writes all results
1b446896 751 UInt_t ii;
752 for(ii = 0;ii<fNParticleFunctions;ii++)
753 fParticleFunctions[ii]->Write();
754
755 for(ii = 0;ii<fNTrackFunctions;ii++)
756 fTrackFunctions[ii]->Write();
757
758 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
759 fParticleAndTrackFunctions[ii]->Write();
5c58441a 760
761 for(ii = 0;ii<fNParticleMonitorFunctions;ii++)
762 fParticleMonitorFunctions[ii]->Write();
763
764 for(ii = 0;ii<fNTrackMonitorFunctions;ii++)
765 fTrackMonitorFunctions[ii]->Write();
766
767 for(ii = 0;ii<fNParticleAndTrackMonitorFunctions;ii++)
768 fParticleAndTrackMonitorFunctions[ii]->Write();
1b446896 769}
770/*************************************************************************************/
491d1b5d 771
1b446896 772void AliHBTAnalysis::SetGlobalPairCut(AliHBTPairCut* cut)
773{
81b7b887 774//Sets the global cut
1b446896 775 if (cut == 0x0)
776 {
777 Error("AliHBTAnalysis::SetGlobalPairCut","Pointer is NULL. Ignoring");
778 }
779 delete fPairCut;
780 fPairCut = (AliHBTPairCut*)cut->Clone();
781}
782
783/*************************************************************************************/
491d1b5d 784
27b3fe5d 785void AliHBTAnalysis::AddTrackFunction(AliHBTOnePairFctn* f)
1b446896 786{
81b7b887 787//Adds track function
1b446896 788 if (f == 0x0) return;
789 if (fNTrackFunctions == fgkFctnArraySize)
790 {
791 Error("AliHBTAnalysis::AddTrackFunction","Can not add this function, not enough place in the array.");
792 }
793 fTrackFunctions[fNTrackFunctions] = f;
794 fNTrackFunctions++;
795}
491d1b5d 796/*************************************************************************************/
797
27b3fe5d 798void AliHBTAnalysis::AddParticleFunction(AliHBTOnePairFctn* f)
1b446896 799{
81b7b887 800//adds particle function
1b446896 801 if (f == 0x0) return;
802
803 if (fNParticleFunctions == fgkFctnArraySize)
804 {
805 Error("AliHBTAnalysis::AddParticleFunction","Can not add this function, not enough place in the array.");
806 }
807 fParticleFunctions[fNParticleFunctions] = f;
808 fNParticleFunctions++;
1b446896 809}
5c58441a 810/*************************************************************************************/
811
27b3fe5d 812void AliHBTAnalysis::AddParticleAndTrackFunction(AliHBTTwoPairFctn* f)
1b446896 813{
81b7b887 814//add resolution function
1b446896 815 if (f == 0x0) return;
816 if (fNParticleAndTrackFunctions == fgkFctnArraySize)
817 {
818 Error("AliHBTAnalysis::AddParticleAndTrackFunction","Can not add this function, not enough place in the array.");
819 }
820 fParticleAndTrackFunctions[fNParticleAndTrackFunctions] = f;
821 fNParticleAndTrackFunctions++;
822}
5c58441a 823/*************************************************************************************/
824
825void AliHBTAnalysis::AddParticleMonitorFunction(AliHBTMonOneParticleFctn* f)
826{
81b7b887 827//add particle monitoring function
5c58441a 828 if (f == 0x0) return;
829
830 if (fNParticleMonitorFunctions == fgkFctnArraySize)
831 {
832 Error("AliHBTAnalysis::AddParticleMonitorFunction","Can not add this function, not enough place in the array.");
833 }
834 fParticleMonitorFunctions[fNParticleMonitorFunctions] = f;
835 fNParticleMonitorFunctions++;
836}
837/*************************************************************************************/
1b446896 838
5c58441a 839void AliHBTAnalysis::AddTrackMonitorFunction(AliHBTMonOneParticleFctn* f)
840{
81b7b887 841//add track monitoring function
5c58441a 842 if (f == 0x0) return;
1b446896 843
5c58441a 844 if (fNTrackMonitorFunctions == fgkFctnArraySize)
845 {
846 Error("AliHBTAnalysis::AddTrackMonitorFunction","Can not add this function, not enough place in the array.");
847 }
848 fTrackMonitorFunctions[fNTrackMonitorFunctions] = f;
849 fNTrackMonitorFunctions++;
850}
1b446896 851/*************************************************************************************/
852
5c58441a 853void AliHBTAnalysis::AddParticleAndTrackMonitorFunction(AliHBTMonTwoParticleFctn* f)
854{
81b7b887 855//add resolution monitoring function
5c58441a 856 if (f == 0x0) return;
857 if (fNParticleAndTrackMonitorFunctions == fgkFctnArraySize)
858 {
859 Error("AliHBTAnalysis::AddParticleAndTrackMonitorFunction","Can not add this function, not enough place in the array.");
860 }
861 fParticleAndTrackMonitorFunctions[fNParticleAndTrackMonitorFunctions] = f;
862 fNParticleAndTrackMonitorFunctions++;
863}
864
1b446896 865
5c58441a 866/*************************************************************************************/
1b446896 867/*************************************************************************************/
491d1b5d 868
1b446896 869Bool_t AliHBTAnalysis::RunCoherencyCheck()
870{
871 //Checks if both HBTRuns are similar
872 //return true if error found
873 //if they seem to be OK return false
874 Int_t i;
81b7b887 875 Info("RunCoherencyCheck","Checking HBT Runs Coherency");
876
877 Info("RunCoherencyCheck","Number of events ...");
1b446896 878 if (fReader->GetNumberOfPartEvents() == fReader->GetNumberOfTrackEvents() ) //check whether there is the same number of events
879 {
81b7b887 880 Info("RunCoherencyCheck","OK. %d found\n",fReader->GetNumberOfTrackEvents());
1b446896 881 }
882 else
883 { //if not the same - ERROR
884 Error("AliHBTAnalysis::RunCoherencyCheck()",
885 "Number of simulated events (%d) is not equal to number of reconstructed events(%d)",
886 fReader->GetNumberOfPartEvents(),fReader->GetNumberOfTrackEvents());
887 return kTRUE;
888 }
889
81b7b887 890 Info("RunCoherencyCheck","Checking number of Particles AND Particles Types in each event ...");
1b446896 891
892 AliHBTEvent *partEvent;
893 AliHBTEvent *trackEvent;
894 for( i = 0; i<fReader->GetNumberOfTrackEvents();i++)
895 {
896 partEvent= fReader->GetParticleEvent(i); //gets the "ith" event
897 trackEvent = fReader->GetTrackEvent(i);
898
899 if ( (partEvent == 0x0) && (partEvent == 0x0) ) continue;
900 if ( (partEvent == 0x0) || (partEvent == 0x0) )
901 {
902 Error("AliHBTAnalysis::RunCoherencyCheck()",
903 "One event is NULL and the other one not. Event Number %d",i);
904 return kTRUE;
905 }
906
907 if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
908 {
909 Error("AliHBTAnalysis::RunCoherencyCheck()",
910 "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
911 i,partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
912 return kTRUE;
913 }
914 else
915 for (Int_t j = 0; j<partEvent->GetNumberOfParticles(); j++)
916 {
917 if( partEvent->GetParticle(j)->GetPdgCode() != trackEvent->GetParticle(j)->GetPdgCode() )
918 {
919 Error("AliHBTAnalysis::RunCoherencyCheck()",
920 "Event %d: Particle %d: PID of simulated particle (%d) not the same of reconstructed track (%d)",
921 i,j, partEvent->GetParticle(j)->GetPdgCode(),trackEvent->GetParticle(j)->GetPdgCode() );
922 return kTRUE;
923
924 }
925 }
926 }
81b7b887 927 Info("RunCoherencyCheck"," Done");
928 Info("RunCoherencyCheck"," Everything looks OK");
1b446896 929 return kFALSE;
930}
931
dc2c3f36 932/*************************************************************************************/
933
934void AliHBTAnalysis::ProcessTracksAndParticlesNonIdentAnal()
935{
81b7b887 936//Performs analysis for both, tracks and particles
937
dc2c3f36 938 AliHBTParticle * part1, * part2;
939 AliHBTParticle * track1, * track2;
940
aeba88d7 941 AliHBTEvent * trackEvent1=0x0,*partEvent1=0x0;
942 AliHBTEvent * trackEvent2=0x0,*partEvent2=0x0;
943 AliHBTEvent * trackEvent3=0x0,*partEvent3=0x0;
dc2c3f36 944
bed069a4 945 AliHBTEvent * rawtrackEvent, *rawpartEvent;//this we get from Reader
dc2c3f36 946
dc2c3f36 947 AliHBTPair * trackpair = new AliHBTPair();
948 AliHBTPair * partpair = new AliHBTPair();
949
bed069a4 950 AliHBTEventBuffer partbuffer(fBufferSize);
951 AliHBTEventBuffer trackbuffer(fBufferSize);
952
953 register UInt_t ii;
dc2c3f36 954
955 trackEvent1 = new AliHBTEvent();
956 partEvent1 = new AliHBTEvent();
957 trackEvent1->SetOwner(kFALSE);
958 partEvent1->SetOwner(kFALSE);;
959
bed069a4 960 Bool_t nocorrfctns = (fNParticleFunctions == 0) && (fNTrackFunctions == 0) && (fNParticleAndTrackFunctions == 0);
961
962 fReader->Rewind();
963
81b7b887 964 Info("ProcessTracksAndParticlesNonIdentAnal","**************************************");
965 Info("ProcessTracksAndParticlesNonIdentAnal","***** NON IDENT MODE ****************");
966 Info("ProcessTracksAndParticlesNonIdentAnal","**************************************");
967
bed069a4 968 for (Int_t i = 0;;i++)//infinite loop
dc2c3f36 969 {
bed069a4 970 if (fReader->Next()) break; //end when no more events available
971
972 rawpartEvent = fReader->GetParticleEvent();
973 rawtrackEvent = fReader->GetTrackEvent();
dc2c3f36 974 if ( (rawpartEvent == 0x0) || (rawtrackEvent == 0x0) ) continue;//in case of any error
bed069a4 975
976 if ( rawpartEvent->GetNumberOfParticles() != rawtrackEvent->GetNumberOfParticles() )
977 {
978 Fatal("ProcessTracksAndParticlesNonIdentAnal",
979 "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
980 i,rawpartEvent->GetNumberOfParticles() , rawtrackEvent->GetNumberOfParticles());
981 return;
982 }
983
dc2c3f36 984 /********************************/
985 /* Filtering out */
986 /********************************/
bed069a4 987 if ( ( (partEvent2==0x0) || (trackEvent2==0x0)) )//in case fBufferSize == 0 and pointers are created do not eneter
988 {
81b7b887 989 partEvent2 = new AliHBTEvent();
990 trackEvent2 = new AliHBTEvent();
bed069a4 991 partEvent2->SetOwner(kTRUE);
992 trackEvent2->SetOwner(kTRUE);
993 }
994
dc2c3f36 995 FilterOut(partEvent1, partEvent2, rawpartEvent, trackEvent1, trackEvent2, rawtrackEvent);
996
997 for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
998 {
81b7b887 999 if ( (j%fDisplayMixingInfo) == 0)
1000 Info("ProcessTracksAndParticlesNonIdentAnal",
1001 "Mixing particle %d from event %d with particles from event %d",j,i,i);
dc2c3f36 1002
1003 part1= partEvent1->GetParticle(j);
1004 track1= trackEvent1->GetParticle(j);
4d91c73a 1005
1006//PID reconstruction imperfections
1007// if( part1->GetPdgCode() != track1->GetPdgCode() )
1008// {
1009// Fatal("ProcessTracksAndParticlesNonIdentAnal",
1010// "Event %d: Particle %d: PID of simulated particle (%d) not the same of reconstructed track (%d)",
1011// i,j, part1->GetPdgCode(),track1->GetPdgCode() );
1012// return;
1013// }
bed069a4 1014
81b7b887 1015 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
e4f2b1da 1016 fParticleMonitorFunctions[ii]->Process(part1);
81b7b887 1017 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
e4f2b1da 1018 fTrackMonitorFunctions[ii]->Process(track1);
81b7b887 1019 for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
e4f2b1da 1020 fParticleAndTrackMonitorFunctions[ii]->Process(track1,part1);
dc2c3f36 1021
bed069a4 1022 if (nocorrfctns) continue;
1023
dc2c3f36 1024 /***************************************/
1025 /****** filling numerators ********/
1026 /****** (particles from event2) ********/
1027 /***************************************/
bed069a4 1028
1029 for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++) //partEvent1 and partEvent2 are particles from the same event but separated to two groups
dc2c3f36 1030 {
1031 part2= partEvent2->GetParticle(k);
d74e6a27 1032 if (part1->GetUID() == part2->GetUID()) continue;//this is the same particle but with different PID
dc2c3f36 1033 partpair->SetParticles(part1,part2);
1034
1035 track2= trackEvent2->GetParticle(k);
1036 trackpair->SetParticles(track1,track2);
1037
777ecb15 1038 if( (fPairCut->PassPairProp(trackpair)) ) //check pair cut
81b7b887 1039 { //do not meets crietria of the pair cut
1040 continue;
1041 }
dc2c3f36 1042 else
1043 {//meets criteria of the pair cut
dc2c3f36 1044 for(ii = 0;ii<fNParticleFunctions;ii++)
1045 fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
1046
1047 for(ii = 0;ii<fNTrackFunctions;ii++)
1048 fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
1049
1050 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
1051 fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(trackpair,partpair);
1052 }
1053 }
5c58441a 1054
81b7b887 1055 if ( fBufferSize == 0) continue;//do not mix diff histograms
dc2c3f36 1056 /***************************************/
1057 /***** Filling denominators *********/
1058 /***************************************/
bed069a4 1059 partbuffer.ResetIter();
1060 trackbuffer.ResetIter();
1061
dc2c3f36 1062 Int_t nmonitor = 0;
1063
bed069a4 1064 while ( (partEvent3 = partbuffer.Next() ) != 0x0)
dc2c3f36 1065 {
bed069a4 1066 trackEvent3 = trackbuffer.Next();
1067
81b7b887 1068 if ( (j%fDisplayMixingInfo) == 0)
1069 Info("ProcessTracksAndParticlesNonIdentAnal",
1070 "Mixing particle %d from event %d with particles from event %d",j,i,i-(++nmonitor));
dc2c3f36 1071
1072 for (Int_t k = 0; k < partEvent3->GetNumberOfParticles() ; k++)
1073 {
1074 part2= partEvent3->GetParticle(k);
1075 partpair->SetParticles(part1,part2);
1b446896 1076
dc2c3f36 1077 track2= trackEvent3->GetParticle(k);
1078 trackpair->SetParticles(track1,track2);
1079
777ecb15 1080 if( (fPairCut->PassPairProp(trackpair)) ) //check pair cut
81b7b887 1081 { //do not meets crietria of the pair cut
1082 continue;
1083 }
dc2c3f36 1084 else
1085 {//meets criteria of the pair cut
1086 UInt_t ii;
1087 for(ii = 0;ii<fNParticleFunctions;ii++)
1088 fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
1089
1090 for(ii = 0;ii<fNTrackFunctions;ii++)
1091 fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
1092
1093 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
1094 fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(trackpair,partpair);
1095 }
1096 }// for particles event2
1097 }//while event2
1098 }//for over particles in event1
1099
bed069a4 1100 partEvent2 = partbuffer.Push(partEvent2);
1101 trackEvent2 = trackbuffer.Push(trackEvent2);
dc2c3f36 1102
1103 }//end of loop over events (1)
bed069a4 1104
1105 partbuffer.SetOwner(kTRUE);
1106 trackbuffer.SetOwner(kTRUE);
1107
b70eb506 1108 delete partEvent1;
1109 delete trackEvent1;
bed069a4 1110 delete partEvent2;
1111 delete trackEvent2;
b70eb506 1112 delete partpair;
1113 delete trackpair;
81b7b887 1114
dc2c3f36 1115}
1b446896 1116/*************************************************************************************/
1117
dc2c3f36 1118void AliHBTAnalysis::ProcessTracksNonIdentAnal()
1119{
2dc7203b 1120//Process Tracks only with non identical mode
dc2c3f36 1121 AliHBTParticle * track1, * track2;
1122
aeba88d7 1123 AliHBTEvent * trackEvent1=0x0;
1124 AliHBTEvent * trackEvent2=0x0;
1125 AliHBTEvent * trackEvent3=0x0;
dc2c3f36 1126
1127 AliHBTEvent * rawtrackEvent;
dc2c3f36 1128
dc2c3f36 1129 AliHBTPair * trackpair = new AliHBTPair();
bed069a4 1130 AliHBTEventBuffer trackbuffer(fBufferSize);
dc2c3f36 1131
81b7b887 1132 register UInt_t ii;
dc2c3f36 1133
1134 trackEvent1 = new AliHBTEvent();
1135 trackEvent1->SetOwner(kFALSE);
1136
bed069a4 1137 fReader->Rewind();
1138
81b7b887 1139 Info("ProcessTracksNonIdentAnal","**************************************");
1140 Info("ProcessTracksNonIdentAnal","***** NON IDENT MODE ****************");
1141 Info("ProcessTracksNonIdentAnal","**************************************");
1142
bed069a4 1143
1144 for (Int_t i = 0;;i++)//infinite loop
dc2c3f36 1145 {
bed069a4 1146 if (fReader->Next()) break; //end when no more events available
1147 rawtrackEvent = fReader->GetTrackEvent();
1148
dc2c3f36 1149 if (rawtrackEvent == 0x0) continue;//in case of any error
1150
1151 /********************************/
1152 /* Filtering out */
1153 /********************************/
bed069a4 1154 if ( trackEvent2==0x0 )//in case fBufferSize == 0 and pointers are created do not eneter
1155 {
81b7b887 1156 trackEvent2 = new AliHBTEvent();
bed069a4 1157 trackEvent2->SetOwner(kTRUE);
1158 }
1159
dc2c3f36 1160 FilterOut(trackEvent1, trackEvent2, rawtrackEvent);
1161
1162 for (Int_t j = 0; j<trackEvent1->GetNumberOfParticles() ; j++)
1163 {
81b7b887 1164 if ( (j%fDisplayMixingInfo) == 0)
1165 Info("ProcessTracksNonIdentAnal",
1166 "Mixing particle %d from event %d with particles from event %d",j,i,i);
dc2c3f36 1167
1168 track1= trackEvent1->GetParticle(j);
1169
81b7b887 1170 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
e4f2b1da 1171 fTrackMonitorFunctions[ii]->Process(track1);
bed069a4 1172
1173 if (fNTrackFunctions == 0x0) continue;
1174
dc2c3f36 1175 /***************************************/
1176 /****** filling numerators ********/
1177 /****** (particles from event2) ********/
1178 /***************************************/
1179 for (Int_t k = 0; k < trackEvent2->GetNumberOfParticles() ; k++)
1180 {
1181 track2= trackEvent2->GetParticle(k);
d74e6a27 1182 if (track1->GetUID() == track2->GetUID()) continue;//this is the same particle but with different PID
dc2c3f36 1183 trackpair->SetParticles(track1,track2);
1184
1185
1186 if( fPairCut->PassPairProp(trackpair)) //check pair cut
1187 { //do not meets crietria of the pair cut
1188 continue;
1189 }
1190 else
1191 {//meets criteria of the pair cut
1192 UInt_t ii;
1193 for(ii = 0;ii<fNTrackFunctions;ii++)
1194 fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
1195 }
1196 }
81b7b887 1197 if ( fBufferSize == 0) continue;//do not mix diff histograms
dc2c3f36 1198 /***************************************/
1199 /***** Filling denominators *********/
1200 /***************************************/
bed069a4 1201 trackbuffer.ResetIter();
dc2c3f36 1202 Int_t nmonitor = 0;
1203
bed069a4 1204 while ( (trackEvent3 = trackbuffer.Next() ) != 0x0)
dc2c3f36 1205 {
1206
81b7b887 1207 if ( (j%fDisplayMixingInfo) == 0)
1208 Info("ProcessTracksNonIdentAnal",
1209 "Mixing particle %d from event %d with particles from event %d",j,i,i-(++nmonitor));
dc2c3f36 1210
1211 for (Int_t k = 0; k < trackEvent3->GetNumberOfParticles() ; k++)
1212 {
1213
1214 track2= trackEvent3->GetParticle(k);
1215 trackpair->SetParticles(track1,track2);
1216
1217
1218 if( fPairCut->PassPairProp(trackpair)) //check pair cut
81b7b887 1219 { //do not meets crietria of the pair cut
1220 continue;
1221 }
dc2c3f36 1222 else
1223 {//meets criteria of the pair cut
81b7b887 1224 for(ii = 0;ii<fNTrackFunctions;ii++)
1225 fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
dc2c3f36 1226 }
1227 }// for particles event2
1228 }//while event2
1229 }//for over particles in event1
1230
bed069a4 1231 trackEvent2 = trackbuffer.Push(trackEvent2);
dc2c3f36 1232
1233 }//end of loop over events (1)
b70eb506 1234
bed069a4 1235 trackbuffer.SetOwner(kTRUE);
b70eb506 1236 delete trackpair;
1237 delete trackEvent1;
bed069a4 1238 delete trackEvent2;
dc2c3f36 1239}
1240/*************************************************************************************/
1241
1242void AliHBTAnalysis::ProcessParticlesNonIdentAnal()
1243{
2dc7203b 1244//process paricles only with non identical mode
dc2c3f36 1245 AliHBTParticle * part1 = 0x0, * part2 = 0x0;
1246
1247 AliHBTEvent * partEvent1 = 0x0;
1248 AliHBTEvent * partEvent2 = 0x0;
1249 AliHBTEvent * partEvent3 = 0x0;
1250
1251 AliHBTEvent * rawpartEvent = 0x0;
1252
dc2c3f36 1253 AliHBTPair * partpair = new AliHBTPair();
bed069a4 1254 AliHBTEventBuffer partbuffer(fBufferSize);
dc2c3f36 1255
bed069a4 1256 register UInt_t ii;
dc2c3f36 1257
1258 partEvent1 = new AliHBTEvent();
bed069a4 1259 partEvent1->SetOwner(kFALSE);
1260
1261 fReader->Rewind();
dc2c3f36 1262
81b7b887 1263 Info("ProcessParticlesNonIdentAnal","**************************************");
1264 Info("ProcessParticlesNonIdentAnal","***** NON IDENT MODE ****************");
1265 Info("ProcessParticlesNonIdentAnal","**************************************");
dc2c3f36 1266
bed069a4 1267 for (Int_t i = 0;;i++)//infinite loop
dc2c3f36 1268 {
bed069a4 1269 if (fReader->Next()) break; //end when no more events available
1270
1271 rawpartEvent = fReader->GetParticleEvent();
dc2c3f36 1272 if ( rawpartEvent == 0x0 ) continue;//in case of any error
1273
1274 /********************************/
1275 /* Filtering out */
1276 /********************************/
bed069a4 1277 if (partEvent2==0x0)//in case fBufferSize == 0 and pointers are created do not eneter
1278 {
b70eb506 1279 partEvent2 = new AliHBTEvent();
bed069a4 1280 partEvent2->SetOwner(kTRUE);
1281 }
1282
dc2c3f36 1283 FilterOut(partEvent1, partEvent2, rawpartEvent);
1284
1285 for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
1286 {
81b7b887 1287 if ( (j%fDisplayMixingInfo) == 0)
1288 Info("ProcessParticlesNonIdentAnal",
1289 "Mixing particle %d from event %d with particles from event %d",j,i,i);
dc2c3f36 1290
1291 part1= partEvent1->GetParticle(j);
1292
bed069a4 1293 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
1294 fParticleMonitorFunctions[ii]->Process(part1);
1295
1296 if (fNParticleFunctions == 0) continue;
1297
dc2c3f36 1298 /***************************************/
1299 /****** filling numerators ********/
1300 /****** (particles from event2) ********/
1301 /***************************************/
1302 for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++)
1303 {
1304 part2= partEvent2->GetParticle(k);
d74e6a27 1305 if (part1->GetUID() == part2->GetUID()) continue;//this is the same particle but with different PID
dc2c3f36 1306 partpair->SetParticles(part1,part2);
1307
1308 if(fPairCut->PassPairProp(partpair) ) //check pair cut
1309 { //do not meets crietria of the pair cut
1310 continue;
1311 }
1312 else
1313 {//meets criteria of the pair cut
dc2c3f36 1314 for(ii = 0;ii<fNParticleFunctions;ii++)
5c58441a 1315 fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
dc2c3f36 1316 }
1317 }
b70eb506 1318 if ( fBufferSize == 0) continue;//do not mix diff histograms
dc2c3f36 1319 /***************************************/
1320 /***** Filling denominators *********/
1321 /***************************************/
bed069a4 1322 partbuffer.ResetIter();
dc2c3f36 1323 Int_t nmonitor = 0;
1324
bed069a4 1325 while ( (partEvent3 = partbuffer.Next() ) != 0x0)
dc2c3f36 1326 {
81b7b887 1327 if ( (j%fDisplayMixingInfo) == 0)
1328 Info("ProcessParticlesNonIdentAnal",
1329 "Mixing particle %d from event %d with particles from event %d",j,i,i-(++nmonitor));
1330
dc2c3f36 1331 for (Int_t k = 0; k < partEvent3->GetNumberOfParticles() ; k++)
1332 {
1333 part2= partEvent3->GetParticle(k);
1334 partpair->SetParticles(part1,part2);
1335
dc2c3f36 1336 if(fPairCut->PassPairProp(partpair) ) //check pair cut
1337 { //do not meets crietria of the pair cut
1338 continue;
1339 }
1340 else
1341 {//meets criteria of the pair cut
dc2c3f36 1342 for(ii = 0;ii<fNParticleFunctions;ii++)
5c58441a 1343 {
5c58441a 1344 fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
1345 }
dc2c3f36 1346 }
1347 }// for particles event2
1348 }//while event2
1349 }//for over particles in event1
bed069a4 1350 partEvent2 = partbuffer.Push(partEvent2);
dc2c3f36 1351 }//end of loop over events (1)
bed069a4 1352
1353 partbuffer.SetOwner(kTRUE);
b70eb506 1354 delete partpair;
1355 delete partEvent1;
bed069a4 1356 delete partEvent2;
dc2c3f36 1357}
1358
1359/*************************************************************************************/
1360void AliHBTAnalysis::FilterOut(AliHBTEvent* outpart1, AliHBTEvent* outpart2, AliHBTEvent* inpart,
1361 AliHBTEvent* outtrack1, AliHBTEvent* outtrack2, AliHBTEvent* intrack)
1362{
1363 //Puts particles accepted as a first particle by global cut in out1
1364 //and as a second particle in out2
1365
1366 AliHBTParticle* part, *track;
1367
1368 outpart1->Reset();
1369 outpart2->Reset();
1370 outtrack1->Reset();
1371 outtrack2->Reset();
1372
1373 AliHBTParticleCut *cut1 = fPairCut->GetFirstPartCut();
1374 AliHBTParticleCut *cut2 = fPairCut->GetSecondPartCut();
1375
1376 Bool_t in1, in2;
1377
1378 for (Int_t i = 0; i < inpart->GetNumberOfParticles(); i++)
1379 {
1380 in1 = in2 = kTRUE;
1381 part = inpart->GetParticle(i);
1382 track = intrack->GetParticle(i);
1383
777ecb15 1384 if ( (cut1->Pass(track)) ) in1 = kFALSE; //if part is rejected by cut1, in1 is false
1385 if ( (cut2->Pass(track)) ) in2 = kFALSE; //if part is rejected by cut2, in2 is false
dc2c3f36 1386
1387 if (gDebug)//to be removed in real analysis
1388 if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
1389 {
1390 //Particle accpted by both cuts
1391 Error("FilterOut","Particle accepted by both cuts");
1392 continue;
1393 }
1394
1395 if (in1)
1396 {
1397 outpart1->AddParticle(part);
1398 outtrack1->AddParticle(track);
1399 continue;
1400 }
1401
1402 if (in2)
1403 {
bed069a4 1404 outpart2->AddParticle(new AliHBTParticle(*part));
1405 outtrack2->AddParticle(new AliHBTParticle(*track));
dc2c3f36 1406 continue;
1407 }
1408 }
dc2c3f36 1409}
1b446896 1410/*************************************************************************************/
81b7b887 1411
dc2c3f36 1412void AliHBTAnalysis::FilterOut(AliHBTEvent* out1, AliHBTEvent* out2, AliHBTEvent* in)
1413{
1414 //Puts particles accepted as a first particle by global cut in out1
1415 //and as a second particle in out2
1416 AliHBTParticle* part;
1417
1418 out1->Reset();
1419 out2->Reset();
1420
1421 AliHBTParticleCut *cut1 = fPairCut->GetFirstPartCut();
1422 AliHBTParticleCut *cut2 = fPairCut->GetSecondPartCut();
1423
1424 Bool_t in1, in2;
1425
1426 for (Int_t i = 0; i < in->GetNumberOfParticles(); i++)
1427 {
1428 in1 = in2 = kTRUE;
1429 part = in->GetParticle(i);
1430
1431 if ( cut1->Pass(part) ) in1 = kFALSE; //if part is rejected by cut1, in1 is false
1432 if ( cut2->Pass(part) ) in2 = kFALSE; //if part is rejected by cut2, in2 is false
1433
1434 if (gDebug)//to be removed in real analysis
1435 if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
1436 {
1437 //Particle accpted by both cuts
1438 Error("FilterOut","Particle accepted by both cuts");
1439 continue;
1440 }
1b446896 1441
dc2c3f36 1442 if (in1)
1443 {
1444 out1->AddParticle(part);
1445 continue;
1446 }
1447
1448 if (in2)
1449 {
1450 out2->AddParticle(part);
1451 continue;
1452 }
1453 }
1454}
1455/*************************************************************************************/
1b446896 1456
dc2c3f36 1457Bool_t AliHBTAnalysis::IsNonIdentAnalysis()
1458{
1459 //checks if it is possible to use special analysis for non identical particles
1460 //it means - in global pair cut first particle id is different than second one
1461 //and both are different from 0
1462 //in the future is possible to perform more sophisticated check
1463 //if cuts have excluding requirements
1464
5c58441a 1465 if (fPairCut->IsEmpty())
1466 return kFALSE;
1467
1468 if (fPairCut->GetFirstPartCut()->IsEmpty())
1469 return kFALSE;
1470
1471 if (fPairCut->GetSecondPartCut()->IsEmpty())
1472 return kFALSE;
dc2c3f36 1473
1474 Int_t id1 = fPairCut->GetFirstPartCut()->GetPID();
1475 Int_t id2 = fPairCut->GetSecondPartCut()->GetPID();
dc2c3f36 1476
5c58441a 1477 if ( (id1==0) || (id2==0) || (id1==id2) )
1478 return kFALSE;
1479
dc2c3f36 1480 return kTRUE;
1481}
9616170a 1482
1483/*************************************************************************************/
1484
1485void AliHBTAnalysis::PressAnyKey()
1486{ //small utility function that helps to make comfortable macros
1487 char c;
1488 int nread = -1;
1489 fcntl(0, F_SETFL, O_NONBLOCK);
1490 ::Info("","Press Any Key to continue ...");
1491 while (nread<1)
1492 {
1493 nread = read(0, &c, 1);
1494 gSystem->ProcessEvents();
1495 }
1496}
1497
1498