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