Updated to the new way AOD stores particles
[u/mrichter/AliRoot.git] / HBTAN / AliHBTAnalysis.cxx
CommitLineData
0aca58be 1#include "AliHBTAnalysis.h"
bfb09ece 2//_________________________________________________________
3////////////////////////////////////////////////////////////////////////////
4//
5// class AliHBTAnalysis
6//
7// Central Object Of HBTAnalyser:
8// This class performs main looping within HBT Analysis
78d7c6d3 9// User must plug a reader of Type AliReader
bfb09ece 10// User plugs in coorelation and monitor functions
11// as well as monitor functions
12//
13// HBT Analysis Tool, which is integral part of AliRoot,
14// ALICE Off-Line framework:
15//
16// Piotr.Skowronski@cern.ch
17// more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
18//
19////////////////////////////////////////////////////////////////////////////
20//_________________________________________________________
21
78d7c6d3 22
8fba7c63 23#include <TSystem.h>
24#include <TFile.h>
25
78d7c6d3 26#include "AliAOD.h"
27#include "AliAODParticle.h"
28#include "AliAODPairCut.h"
7b6503d6 29#include "AliEventCut.h"
78d7c6d3 30
31#include "AliEventBuffer.h"
32
33#include "AliReader.h"
1b446896 34#include "AliHBTPair.h"
1b446896 35#include "AliHBTFunction.h"
5c58441a 36#include "AliHBTMonitorFunction.h"
81b7b887 37
e4f2b1da 38
1b446896 39ClassImp(AliHBTAnalysis)
40
41const UInt_t AliHBTAnalysis::fgkFctnArraySize = 100;
81b7b887 42const UInt_t AliHBTAnalysis::fgkDefaultMixingInfo = 1000;
43const Int_t AliHBTAnalysis::fgkDefaultBufferSize = 5;
44
45AliHBTAnalysis::AliHBTAnalysis():
090e46d6 46 fProcEvent(0x0),
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]),
7b6503d6 60 fBkgEventCut(0x0),
5994509d 61 fPartBuffer(0x0),
62 fTrackBuffer(0x0),
2dc7203b 63 fBufferSize(2),
64 fDisplayMixingInfo(fgkDefaultMixingInfo),
66d1d1a4 65 fIsOwner(kFALSE),
090e46d6 66 fProcessOption(kSimulatedAndReconstructed),
e92ecbdf 67 fNoCorrfctns(kFALSE),
5994509d 68 fOutputFileName(0x0),
69 fVertexX(0.0),
70 fVertexY(0.0),
71 fVertexZ(0.0)
1b446896 72 {
9616170a 73 //default constructor
66d1d1a4 74
1b446896 75 }
491d1b5d 76/*************************************************************************************/
1b446896 77
81b7b887 78AliHBTAnalysis::AliHBTAnalysis(const AliHBTAnalysis& in):
78d7c6d3 79 AliAnalysis(in),
090e46d6 80 fProcEvent(0x0),
2dc7203b 81 fReader(0x0),
82 fNTrackFunctions(0),
83 fNParticleFunctions(0),
84 fNParticleAndTrackFunctions(0),
85 fNTrackMonitorFunctions(0),
86 fNParticleMonitorFunctions(0),
87 fNParticleAndTrackMonitorFunctions(0),
88 fTrackFunctions(0x0),
89 fParticleFunctions(0x0),
90 fParticleAndTrackFunctions(0x0),
91 fParticleMonitorFunctions(0x0),
92 fTrackMonitorFunctions(0x0),
93 fParticleAndTrackMonitorFunctions(0x0),
7b6503d6 94 fBkgEventCut(0x0),
5994509d 95 fPartBuffer(0x0),
96 fTrackBuffer(0x0),
2dc7203b 97 fBufferSize(fgkDefaultBufferSize),
98 fDisplayMixingInfo(fgkDefaultMixingInfo),
66d1d1a4 99 fIsOwner(kFALSE),
090e46d6 100 fProcessOption(kSimulatedAndReconstructed),
e92ecbdf 101 fNoCorrfctns(kFALSE),
5994509d 102 fOutputFileName(0x0),
103 fVertexX(0.0),
104 fVertexY(0.0),
105 fVertexZ(0.0)
81b7b887 106 {
2dc7203b 107//copy constructor
81b7b887 108 Fatal("AliHBTAnalysis(const AliHBTAnalysis&)","Sensless");
109 }
110/*************************************************************************************/
34914285 111AliHBTAnalysis& AliHBTAnalysis::operator=(const AliHBTAnalysis& /*right*/)
81b7b887 112 {
2dc7203b 113//operator =
81b7b887 114 Fatal("AliHBTAnalysis(const AliHBTAnalysis&)","Sensless");
115 return *this;
116 }
117/*************************************************************************************/
1b446896 118AliHBTAnalysis::~AliHBTAnalysis()
119 {
120 //destructor
121 //note that we do not delete functions itself
122 // they should be deleted by whom where created
123 //we only store pointers, and we use "new" only for pointers array
e7a04795 124
125 if (fIsOwner)
126 {
78d7c6d3 127 if (AliVAODParticle::GetDebug()>5)Info("~AliHBTAnalysis","Is Owner: Attempting to delete functions");
e7a04795 128 DeleteFunctions();
78d7c6d3 129 if (AliVAODParticle::GetDebug()>5)Info("~AliHBTAnalysis","Delete functions done");
e7a04795 130 }
1b446896 131 delete [] fTrackFunctions;
132 delete [] fParticleFunctions;
133 delete [] fParticleAndTrackFunctions;
134
5c58441a 135 delete [] fParticleMonitorFunctions;
136 delete [] fTrackMonitorFunctions;
090e46d6 137 delete [] fParticleAndTrackMonitorFunctions;
5c58441a 138
7b6503d6 139 delete fBkgEventCut;
8fba7c63 140 delete fOutputFileName;
1b446896 141 }
142
143/*************************************************************************************/
7b6503d6 144
78d7c6d3 145Int_t AliHBTAnalysis::ProcessEvent(AliAOD* aodrec, AliAOD* aodsim)
146{
147 //Processes one event
090e46d6 148 if (fProcEvent == 0x0)
78d7c6d3 149 {
a2b9dcc7 150 Error("ProcessEvent","Analysis <<%s>> option not specified.",GetName());
78d7c6d3 151 return 1;
152 }
a74ad849 153 if ( Rejected(aodrec,aodsim) )
154 {
155// Double_t x = 0, y = 0, z = 0;
156// aodrec->GetPrimaryVertex(x,y,z);
157// Info("ProcessEvent","Event has vertex at %f %f %f",x,y,z);
158 Info("ProcessEvent","Nch is %d",aodsim->GetNumberOfCharged());
159 Info("ProcessEvent","%s: Event cut rejected this event",GetName());
160 return 0;
161 }
163cad1f 162
5994509d 163 //Move event to the apparent vertex -> must be after the event cut
164 //It is important for any cut that use any spacial coordiantes,
165 //f.g. anti-merging cut in order to preserve the same bahavior of variables (f.g. distance between tracks)
166 Double_t dvx = 0, dvy = 0, dvz = 0;
167 if (aodrec)
168 {
169 Double_t pvx,pvy,pvz;
170 aodrec->GetPrimaryVertex(pvx,pvy,pvz);
171
172 dvx = fVertexX - pvx;
173 dvy = fVertexY - pvy;
174 dvz = fVertexZ - pvz;
175 aodrec->Move(dvx,dvy,dvz);
176 }
177
178 Int_t result = (this->*fProcEvent)(aodrec,aodsim);
179
180 if (aodrec) aodrec->Move(-dvx,-dvy,-dvz);//move event back to the same spacial coordinates
181
182 return result;
78d7c6d3 183}
184/*************************************************************************************/
185
186Int_t AliHBTAnalysis::Finish()
187{
090e46d6 188//Finishes analysis
78d7c6d3 189 WriteFunctions();
e92ecbdf 190 return 0;
78d7c6d3 191}
192/*************************************************************************************/
e4f2b1da 193
81b7b887 194void AliHBTAnalysis::DeleteFunctions()
195{
e4f2b1da 196 //Deletes all functions added to analysis
5994509d 197
81b7b887 198 UInt_t ii;
199 for(ii = 0;ii<fNParticleFunctions;ii++)
e7a04795 200 {
78d7c6d3 201 if (AliVAODParticle::GetDebug()>5)
90520373 202 {
203 Info("DeleteFunctions","Deleting ParticleFunction %#x",fParticleFunctions[ii]);
204 Info("DeleteFunctions","Deleting ParticleFunction %s",fParticleFunctions[ii]->Name());
205 }
e7a04795 206 delete fParticleFunctions[ii];
207 }
e4f2b1da 208 fNParticleFunctions = 0;
81b7b887 209
210 for(ii = 0;ii<fNTrackFunctions;ii++)
e7a04795 211 {
78d7c6d3 212 if (AliVAODParticle::GetDebug()>5)
e7a04795 213 {
90520373 214 Info("DeleteFunctions","Deleting TrackFunction %#x",fTrackFunctions[ii]);
215 Info("DeleteFunctions","Deleting TrackFunction %s",fTrackFunctions[ii]->Name());
e7a04795 216 }
217 delete fTrackFunctions[ii];
218 }
e4f2b1da 219 fNTrackFunctions = 0;
220
81b7b887 221 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
e7a04795 222 {
78d7c6d3 223 if (AliVAODParticle::GetDebug()>5)
90520373 224 {
225 Info("DeleteFunctions","Deleting ParticleAndTrackFunction %#x",fParticleAndTrackFunctions[ii]);
226 Info("DeleteFunctions","Deleting ParticleAndTrackFunction %s",fParticleAndTrackFunctions[ii]->Name());
227 }
e7a04795 228 delete fParticleAndTrackFunctions[ii];
229 }
e4f2b1da 230 fNParticleAndTrackFunctions = 0;
81b7b887 231
232 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
e7a04795 233 {
78d7c6d3 234 if (AliVAODParticle::GetDebug()>5)
90520373 235 {
236 Info("DeleteFunctions","Deleting ParticleMonitorFunction %#x",fParticleMonitorFunctions[ii]);
237 Info("DeleteFunctions","Deleting ParticleMonitorFunction %s",fParticleMonitorFunctions[ii]->Name());
238 }
e7a04795 239 delete fParticleMonitorFunctions[ii];
240 }
e4f2b1da 241 fNParticleMonitorFunctions = 0;
81b7b887 242
243 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
e7a04795 244 {
78d7c6d3 245 if (AliVAODParticle::GetDebug()>5)
90520373 246 {
247 Info("DeleteFunctions","Deleting TrackMonitorFunction %#x",fTrackMonitorFunctions[ii]);
248 Info("DeleteFunctions","Deleting TrackMonitorFunction %s",fTrackMonitorFunctions[ii]->Name());
249 }
e7a04795 250 delete fTrackMonitorFunctions[ii];
251 }
e4f2b1da 252 fNTrackMonitorFunctions = 0;
81b7b887 253
254 for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
e7a04795 255 {
78d7c6d3 256 if (AliVAODParticle::GetDebug()>5)
90520373 257 {
258 Info("DeleteFunctions","Deleting ParticleAndTrackMonitorFunction %#x",fParticleAndTrackMonitorFunctions[ii]);
259 Info("DeleteFunctions","Deleting ParticleAndTrackMonitorFunction %s",fParticleAndTrackMonitorFunctions[ii]->Name());
260 }
e7a04795 261 delete fParticleAndTrackMonitorFunctions[ii];
262 }
e4f2b1da 263 fNParticleAndTrackMonitorFunctions = 0;
bfb09ece 264
81b7b887 265}
e4f2b1da 266/*************************************************************************************/
267
78d7c6d3 268Int_t AliHBTAnalysis::Init()
e4f2b1da 269{
2dc7203b 270//Initializeation method
271//calls Init for all functions
090e46d6 272
273 //Depending on option and pair cut assigns proper analysis method
274 Bool_t nonid = IsNonIdentAnalysis();
275 switch(fProcessOption)
276 {
277 case kReconstructed:
278 if (nonid) fProcEvent = &AliHBTAnalysis::ProcessRecNonId;
279 else fProcEvent = &AliHBTAnalysis::ProcessRec;
a2b9dcc7 280 SetCutsOnRec();
090e46d6 281 break;
282
283 case kSimulated:
284 if (nonid) fProcEvent = &AliHBTAnalysis::ProcessSimNonId;
285 else fProcEvent = &AliHBTAnalysis::ProcessSim;
a2b9dcc7 286 SetCutsOnSim();
090e46d6 287 break;
288
289 case kSimulatedAndReconstructed:
290 if (nonid) fProcEvent = &AliHBTAnalysis::ProcessRecAndSimNonId;
291 else fProcEvent = &AliHBTAnalysis::ProcessRecAndSim;
292 break;
293 }
e92ecbdf 294
295 if (fPartBuffer == 0x0) fPartBuffer = new AliEventBuffer (fBufferSize);
296 else fPartBuffer->Reset();
297
298 if (fTrackBuffer == 0x0) fTrackBuffer = new AliEventBuffer (fBufferSize);
299 else fTrackBuffer->Reset();
300
301
302 fNoCorrfctns = (fNParticleFunctions == 0) && (fNTrackFunctions == 0) && (fNParticleAndTrackFunctions == 0);
303
e4f2b1da 304 UInt_t ii;
305 for(ii = 0;ii<fNParticleFunctions;ii++)
306 fParticleFunctions[ii]->Init();
307
308 for(ii = 0;ii<fNTrackFunctions;ii++)
309 fTrackFunctions[ii]->Init();
310
311 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
312 fParticleAndTrackFunctions[ii]->Init();
313
314 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
315 fParticleMonitorFunctions[ii]->Init();
316
317 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
318 fTrackMonitorFunctions[ii]->Init();
319
320 for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
321 fParticleAndTrackMonitorFunctions[ii]->Init();
bfb09ece 322
78d7c6d3 323 return 0;
e4f2b1da 324}
325/*************************************************************************************/
326
327void AliHBTAnalysis::ResetFunctions()
328{
329//In case fOwner is true, deletes all functions
330//in other case, just set number of analysis to 0
331 if (fIsOwner) DeleteFunctions();
332 else
333 {
334 fNParticleFunctions = 0;
335 fNTrackFunctions = 0;
336 fNParticleAndTrackFunctions = 0;
337 fNParticleMonitorFunctions = 0;
338 fNTrackMonitorFunctions = 0;
339 fNParticleAndTrackMonitorFunctions = 0;
340 }
341}
342/*************************************************************************************/
e92ecbdf 343Int_t AliHBTAnalysis::ProcessRecAndSim(AliAOD* aodrec, AliAOD* aodsim)
344{
090e46d6 345//Does analysis for both tracks and particles
e92ecbdf 346//mainly for resolution study and analysies with weighting algirithms
e92ecbdf 347
348// cut on particles only -- why?
349// - PID: when we make resolution analysis we want to take only tracks with correct PID
090e46d6 350// We need cut on tracks because there are data characteristic
e92ecbdf 351
352 AliVAODParticle * part1, * part2;
353 AliVAODParticle * track1, * track2;
354
355 AliAOD * trackEvent = aodrec, *partEvent = aodsim;
356 AliAOD* trackEvent1 = new AliAOD();
357 AliAOD* partEvent1 = new AliAOD();
358 partEvent1->SetOwner(kTRUE);
359 trackEvent1->SetOwner(kTRUE);
360
361 AliAOD * trackEvent2,*partEvent2;
362
363// Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
364
365// Int_t nev = fReader->GetNumberOfTrackEvents();
090e46d6 366 static AliHBTPair tpair;
367 static AliHBTPair ppair;
e92ecbdf 368
369 AliHBTPair* trackpair = &tpair;
370 AliHBTPair* partpair = &ppair;
371
372 AliHBTPair * tmptrackpair;//temprary pointers to pairs
373 AliHBTPair * tmppartpair;
374
375 register UInt_t ii;
376
377
378
379 if ( !partEvent || !trackEvent )
380 {
a2b9dcc7 381 Error("ProcessRecAndSim","<<%s>> Can not get event",GetName());
e92ecbdf 382 return 1;
383 }
384
385 if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
386 {
387 Error("ProcessRecAndSim",
388 "Number of simulated particles (%d) not equal to number of reconstructed tracks (%d). Skipping Event.",
389 partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
390 return 2;
391 }
392
393
394 for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
395 {
396 /***************************************/
397 /****** Looping same events ********/
398 /****** filling numerators ********/
399 /***************************************/
400 if ( (j%fDisplayMixingInfo) == 0)
401 Info("ProcessTracksAndParticles",
402 "Mixing particle %d with particles from the same event",j);
403
404 part1= partEvent->GetParticle(j);
405 track1= trackEvent->GetParticle(j);
d50bb49a 406
e92ecbdf 407 Bool_t firstcut = (this->*fkPass1)(part1,track1);
408 if (fBufferSize != 0)
409 if ( (firstcut == kFALSE) || ( (this->*fkPass2)(part1,track1) == kFALSE ) )
410 {
411 //accepted by any cut
412 // we have to copy because reader keeps only one event
413
d50bb49a 414 partEvent1->AddParticle(part1);
415 trackEvent1->AddParticle(track1);
e92ecbdf 416 }
417
418 if (firstcut) continue;
419
420 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
421 fParticleMonitorFunctions[ii]->Process(part1);
422 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
423 fTrackMonitorFunctions[ii]->Process(track1);
424 for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
425 fParticleAndTrackMonitorFunctions[ii]->Process(track1,part1);
426
427 if (fNoCorrfctns) continue;
428
429 for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
430 {
431 part2= partEvent->GetParticle(k);
432 if (part1->GetUID() == part2->GetUID()) continue;
433 partpair->SetParticles(part1,part2);
434
435 track2= trackEvent->GetParticle(k);
436 trackpair->SetParticles(track1,track2);
437
438 if( (this->*fkPass)(partpair,trackpair) ) //check pair cut
439 { //do not meets crietria of the pair cut, try with swapped pairs
440 if( (this->*fkPass)((AliHBTPair*)partpair->GetSwappedPair(),(AliHBTPair*)trackpair->GetSwappedPair()) )
441 continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
442 else
443 { //swaped pair meets all the criteria
444 tmppartpair = (AliHBTPair*)partpair->GetSwappedPair();
445 tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
446 }
447 }
448 else
449 {//meets criteria of the pair cut
450 tmptrackpair = trackpair;
451 tmppartpair = partpair;
452 }
453
454 for(ii = 0;ii<fNParticleFunctions;ii++)
455 fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
456
457 for(ii = 0;ii<fNTrackFunctions;ii++)
458 fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
459
460 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
461 fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair,tmppartpair);
462 //end of 2nd loop over particles from the same event
463 }//for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
464
465 /***************************************/
466 /***** Filling denominators *********/
467 /***************************************/
468 if (fBufferSize == 0) continue;
469
470 fPartBuffer->ResetIter();
471 fTrackBuffer->ResetIter();
472 Int_t m = 0;
473 while (( partEvent2 = fPartBuffer->Next() ))
474 {
475 trackEvent2 = fTrackBuffer->Next();
476
477 m++;
478 if ( (j%fDisplayMixingInfo) == 0)
479 Info("ProcessTracksAndParticles",
480 "Mixing particle %d from current event with particles from event %d",j,-m);
481
482 for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++) // ... on all particles
483 {
484 part2= partEvent2->GetParticle(l);
485 partpair->SetParticles(part1,part2);
486
487 track2= trackEvent2->GetParticle(l);
488 trackpair->SetParticles(track1,track2);
489
490 if( (this->*fkPass)(partpair,trackpair) ) //check pair cut
491 { //do not meets crietria of the
492 if( (this->*fkPass)((AliHBTPair*)partpair->GetSwappedPair(),(AliHBTPair*)trackpair->GetSwappedPair()) )
493 continue;
494 else
495 {
496 tmppartpair = (AliHBTPair*)partpair->GetSwappedPair();
497 tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
498 }
499 }
500 else
501 {//meets criteria of the pair cut
502 tmptrackpair = trackpair;
503 tmppartpair = partpair;
504 }
505
506 for(ii = 0;ii<fNParticleFunctions;ii++)
507 fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
508
509 for(ii = 0;ii<fNTrackFunctions;ii++)
510 fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
511
512 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
513 fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair,tmppartpair);
514 }//for(Int_t l = 0; l<N2;l++) // ... on all particles
515
516 }
517 //end of loop over particles from first event
518 }//for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
519 delete fPartBuffer->Push(partEvent1);
520 delete fTrackBuffer->Push(trackEvent1);
521 //end of loop over events
522 return 0;
523}
524/*************************************************************************************/
e4f2b1da 525
e92ecbdf 526Int_t AliHBTAnalysis::ProcessSim(AliAOD* /*aodrec*/, AliAOD* aodsim)
527{
090e46d6 528 //Does analysis of simulated data
529 AliVAODParticle * part1, * part2;
530
531 AliAOD* partEvent = aodsim;
532 AliAOD* partEvent1 = new AliAOD();
533 partEvent1->SetOwner(kTRUE);
534
535 AliAOD* partEvent2;
536
537 AliHBTPair ppair;
538
539 AliHBTPair* partpair = &ppair;
540
541 AliHBTPair * tmppartpair;
542
543 register UInt_t ii;
544
545
546
547 if ( !partEvent )
548 {
549 Error("ProcessRecAndSim","Can not get event");
550 return 1;
551 }
552
553
554 for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
555 {
556 /***************************************/
557 /****** Looping same events ********/
558 /****** filling numerators ********/
559 /***************************************/
560 if ( (j%fDisplayMixingInfo) == 0)
561 Info("ProcessTracksAndParticles",
562 "Mixing particle %d with particles from the same event",j);
563
564 part1= partEvent->GetParticle(j);
565
cea0a066 566 Bool_t firstcut = fPairCut->GetFirstPartCut()->Rejected(part1);
090e46d6 567
568 if (fBufferSize != 0)
cea0a066 569 if ( (firstcut == kFALSE) || ( fPairCut->GetSecondPartCut()->Rejected(part1) == kFALSE ) )
090e46d6 570 {
571 //accepted by any cut
572 // we have to copy because reader keeps only one event
573
d50bb49a 574 partEvent1->AddParticle(part1);
090e46d6 575 }
576
577 if (firstcut) continue;
578
579 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
580 fParticleMonitorFunctions[ii]->Process(part1);
581
582 if ( fNParticleFunctions == 0 ) continue;
583
584 for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
585 {
586 part2= partEvent->GetParticle(k);
587 if (part1->GetUID() == part2->GetUID()) continue;
588 partpair->SetParticles(part1,part2);
589
cea0a066 590 if(fPairCut->Rejected(partpair)) //check pair cut
090e46d6 591 { //do not meets crietria of the
cea0a066 592 if( fPairCut->Rejected((AliHBTPair*)partpair->GetSwappedPair()) ) continue;
090e46d6 593 else tmppartpair = (AliHBTPair*)partpair->GetSwappedPair();
594 }
595 else
596 {
597 tmppartpair = partpair;
598 }
599
600 for(ii = 0;ii<fNParticleFunctions;ii++)
601 fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
602
603 //end of 2nd loop over particles from the same event
604 }//for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
605
606 /***************************************/
607 /***** Filling denominators *********/
608 /***************************************/
609 if (fBufferSize == 0) continue;
610
611 fPartBuffer->ResetIter();
612 Int_t m = 0;
613 while (( partEvent2 = fPartBuffer->Next() ))
614 {
615 m++;
616 if ( (j%fDisplayMixingInfo) == 0)
617 Info("ProcessParticles",
618 "Mixing particle %d from current event with particles from event %d",j,-m);
619 for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++) // ... on all particles
620 {
621
622 part2= partEvent2->GetParticle(l);
623 partpair->SetParticles(part1,part2);
624
cea0a066 625 if( fPairCut->Rejected(partpair) ) //check pair cut
090e46d6 626 { //do not meets crietria of the
cea0a066 627 if( fPairCut->Rejected((AliHBTPair*)partpair->GetSwappedPair()) )
090e46d6 628 continue;
629 else
630 {
631 tmppartpair = (AliHBTPair*)partpair->GetSwappedPair();
632 }
633 }
634 else
635 {//meets criteria of the pair cut
636 tmppartpair = partpair;
637 }
638
639 for(ii = 0;ii<fNParticleFunctions;ii++)
640 fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
641
642 }//for(Int_t l = 0; l<N2;l++) // ... on all particles
643 }
644 }
645 delete fPartBuffer->Push(partEvent1);
646 //end of loop over events
647 return 0;
e92ecbdf 648}
649/*************************************************************************************/
650Int_t AliHBTAnalysis::ProcessRec(AliAOD* aodrec, AliAOD* /*aodsim*/)
651{
090e46d6 652 //Does analysis of reconstructed data
653 AliVAODParticle * track1, * track2;
654
655 AliAOD* trackEvent = aodrec;
656 AliAOD* trackEvent1 = new AliAOD();
657 trackEvent1->SetOwner(kTRUE);
658
659 AliAOD* trackEvent2;
660
661 AliHBTPair tpair;
662
663 AliHBTPair* trackpair = &tpair;
664
665 AliHBTPair * tmptrackpair;
666
667 register UInt_t ii;
668
669
670 if ( !trackEvent )
671 {
672 Error("ProcessRecAndSim","Can not get event");
673 return 1;
674 }
675
676
677 for (Int_t j = 0; j<trackEvent->GetNumberOfParticles() ; j++)
678 {
679 /***************************************/
680 /****** Looping same events ********/
681 /****** filling numerators ********/
682 /***************************************/
683 if ( (j%fDisplayMixingInfo) == 0)
684 Info("ProcessTracksAndParticles",
685 "Mixing Particle %d with Particles from the same event",j);
686
687 track1= trackEvent->GetParticle(j);
688
cea0a066 689 Bool_t firstcut = fPairCut->GetFirstPartCut()->Rejected(track1);
090e46d6 690
691 if (fBufferSize != 0)
cea0a066 692 if ( (firstcut == kFALSE) || ( fPairCut->GetSecondPartCut()->Rejected(track1) == kFALSE ) )
090e46d6 693 {
694 //accepted by any cut
695 // we have to copy because reader keeps only one event
696
d50bb49a 697 trackEvent1->AddParticle(track1);
090e46d6 698 }
699
700 if (firstcut) continue;
701
702 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
703 fParticleMonitorFunctions[ii]->Process(track1);
704
705 if ( fNParticleFunctions == 0 ) continue;
706
707 for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
708 {
709 track2= trackEvent->GetParticle(k);
710 if (track1->GetUID() == track2->GetUID()) continue;
711 trackpair->SetParticles(track1,track2);
712
cea0a066 713 if(fPairCut->Rejected(trackpair)) //check pair cut
090e46d6 714 { //do not meets crietria of the
cea0a066 715 if( fPairCut->Rejected((AliHBTPair*)trackpair->GetSwappedPair()) ) continue;
090e46d6 716 else tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
717 }
718 else
719 {
720 tmptrackpair = trackpair;
721 }
722
723 for(ii = 0;ii<fNTrackFunctions;ii++)
724 fParticleFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
725
726 //end of 2nd loop over Particles from the same event
727 }//for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
728
729 /***************************************/
730 /***** Filling denominators *********/
731 /***************************************/
732 if (fBufferSize == 0) continue;
733
734 fTrackBuffer->ResetIter();
735 Int_t m = 0;
736 while (( trackEvent2 = fTrackBuffer->Next() ))
737 {
738 m++;
739 if ( (j%fDisplayMixingInfo) == 0)
740 Info("ProcessParticles",
741 "Mixing Particle %d from current event with Particles from event %d",j,-m);
742 for(Int_t l = 0; l<trackEvent2->GetNumberOfParticles();l++) // ... on all Particles
743 {
744
745 track2= trackEvent2->GetParticle(l);
746 trackpair->SetParticles(track1,track2);
747
cea0a066 748 if( fPairCut->Rejected(trackpair) ) //check pair cut
090e46d6 749 { //do not meets crietria of the
cea0a066 750 if( fPairCut->Rejected((AliHBTPair*)trackpair->GetSwappedPair()) )
090e46d6 751 continue;
752 else
753 {
754 tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
755 }
756 }
757 else
758 {//meets criteria of the pair cut
759 tmptrackpair = trackpair;
760 }
761
762 for(ii = 0;ii<fNTrackFunctions;ii++)
763 fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
764
765 }//for(Int_t l = 0; l<N2;l++) // ... on all Particles
766 }
767 }
768 delete fTrackBuffer->Push(trackEvent1);
769 //end of loop over events
770 return 0;
e92ecbdf 771}
772/*************************************************************************************/
773
774Int_t AliHBTAnalysis::ProcessRecAndSimNonId(AliAOD* aodrec, AliAOD* aodsim)
775{
090e46d6 776//Analyzes both reconstructed and simulated data
777 if (aodrec == 0x0)
778 {
779 Error("ProcessTracksAndParticlesNonIdentAnal","Reconstructed event is NULL");
780 return 1;
781 }
782
783 if (aodsim == 0x0)
784 {
785 Error("ProcessTracksAndParticlesNonIdentAnal","Simulated event is NULL");
786 return 1;
787 }
788
789 if ( aodrec->GetNumberOfParticles() != aodsim->GetNumberOfParticles() )
790 {
791 Error("ProcessTracksAndParticlesNonIdentAnal",
792 "Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
793 aodsim->GetNumberOfParticles() , aodrec->GetNumberOfParticles());
794 return 2;
795 }
796
797
798 AliVAODParticle * part1, * part2;
799 AliVAODParticle * track1, * track2;
800
801 static AliAOD aodrec1;
802 static AliAOD aodsim1;
803
804 AliAOD * trackEvent1=&aodrec1,*partEvent1=&aodsim1;//Particle that passes first particle cut, this event
805 AliAOD * trackEvent2=0x0,*partEvent2=0x0;//Particle that passes second particle cut, this event
806 AliAOD * trackEvent3=0x0,*partEvent3=0x0;//Particle that passes second particle cut, events from buffer
807
808 AliAOD* rawtrackEvent = aodrec;//this we get from Reader
809 AliAOD* rawpartEvent = aodsim;//this we get from Reader
810
811 static AliHBTPair tpair;
812 static AliHBTPair ppair;
813
814 AliHBTPair* trackpair = &tpair;
815 AliHBTPair* partpair = &ppair;
816
817 register UInt_t ii;
818
819
820 trackEvent1 = new AliAOD();
821 partEvent1 = new AliAOD();
822 trackEvent1->SetOwner(kFALSE);
823 partEvent1->SetOwner(kFALSE);;
824
825 /********************************/
826 /* Filtering out */
827 /********************************/
828 if ( ( (partEvent2==0x0) || (trackEvent2==0x0)) )//in case fBufferSize == 0 and pointers are created do not eneter
829 {
830 partEvent2 = new AliAOD();
831 trackEvent2 = new AliAOD();
832 partEvent2->SetOwner(kTRUE);
833 trackEvent2->SetOwner(kTRUE);
834 }
835
836 FilterOut(partEvent1, partEvent2, rawpartEvent, trackEvent1, trackEvent2, rawtrackEvent);
837
838 for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
839 {
840 if ( (j%fDisplayMixingInfo) == 0)
841 Info("ProcessTracksAndParticlesNonIdentAnal",
842 "Mixing particle %d from current event with particles from current event",j);
843
844 part1= partEvent1->GetParticle(j);
845 track1= trackEvent1->GetParticle(j);
846
847
848 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
849 fParticleMonitorFunctions[ii]->Process(part1);
850 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
851 fTrackMonitorFunctions[ii]->Process(track1);
852 for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
853 fParticleAndTrackMonitorFunctions[ii]->Process(track1,part1);
854
855 if (fNoCorrfctns) continue;
856
857 /***************************************/
858 /****** filling numerators ********/
859 /****** (particles from event2) ********/
860 /***************************************/
861
862 for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++) //partEvent1 and partEvent2 are particles from the same event but separated to two groups
863 {
864 part2= partEvent2->GetParticle(k);
865 if (part1->GetUID() == part2->GetUID()) continue;//this is the same particle but with different PID
866 partpair->SetParticles(part1,part2);
867
868 track2= trackEvent2->GetParticle(k);
869 trackpair->SetParticles(track1,track2);
870
871 if( (this->*fkPassPairProp)(partpair,trackpair) ) //check pair cut
872 { //do not meets crietria of the pair cut
873 continue;
874 }
875 else
876 {//meets criteria of the pair cut
877 for(ii = 0;ii<fNParticleFunctions;ii++)
878 fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
879
880 for(ii = 0;ii<fNTrackFunctions;ii++)
881 fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
882
883 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
884 fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(trackpair,partpair);
885 }
886 }
887
888 if ( fBufferSize == 0) continue;//do not mix diff histograms
889 /***************************************/
890 /***** Filling denominators *********/
891 /***************************************/
892 fPartBuffer->ResetIter();
893 fTrackBuffer->ResetIter();
894
895 Int_t nmonitor = 0;
896
897 while ( (partEvent3 = fPartBuffer->Next() ) != 0x0)
898 {
899 trackEvent3 = fTrackBuffer->Next();
900
901 if ( (j%fDisplayMixingInfo) == 0)
902 Info("ProcessTracksAndParticlesNonIdentAnal",
903 "Mixing particle %d from current event with particles from event%d",j,-(++nmonitor));
904
905 for (Int_t k = 0; k < partEvent3->GetNumberOfParticles() ; k++)
906 {
907 part2= partEvent3->GetParticle(k);
908 partpair->SetParticles(part1,part2);
909
910 track2= trackEvent3->GetParticle(k);
911 trackpair->SetParticles(track1,track2);
912
913 if( (this->*fkPassPairProp)(partpair,trackpair) ) //check pair cut
914 { //do not meets crietria of the pair cut
915 continue;
916 }
917 else
918 {//meets criteria of the pair cut
919 UInt_t ii;
920 for(ii = 0;ii<fNParticleFunctions;ii++)
921 fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
922
923 for(ii = 0;ii<fNTrackFunctions;ii++)
924 fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
925
926 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
927 fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(trackpair,partpair);
928 }
929 }// for particles event2
930 }//while event2
931 }//for over particles in event1
932
933 delete fPartBuffer->Push(partEvent2);
934 delete fTrackBuffer->Push(trackEvent2);
935
936 return 0;
e92ecbdf 937}
938/*************************************************************************************/
090e46d6 939Int_t AliHBTAnalysis::ProcessSimNonId(AliAOD* /*aodrec*/, AliAOD* aodsim)
e92ecbdf 940{
090e46d6 941//does analysis of simulated (MC) data in non-identical mode
942//i.e. when particles selected by first part. cut are a disjunctive set than particles
943//passed by the second part. cut
944 if (aodsim == 0x0)
945 {
946 return 1;
947 }
948
949
950 AliVAODParticle * part1, * part2;
951
952 static AliAOD aodsim1;
953
954 AliAOD* partEvent1=&aodsim1;//Particle that passes first particle cut, this event
955 AliAOD* partEvent2=0x0;//Particle that passes second particle cut, this event
956 AliAOD* partEvent3=0x0;//Particle that passes second particle cut, events from buffer
957
958 AliAOD* rawpartEvent = aodsim;//this we get from Reader
959
960 static AliHBTPair ppair;
961
962 AliHBTPair* partpair = &ppair;
963
964 register UInt_t ii;
965
966
967 partEvent1 = new AliAOD();
968 partEvent1->SetOwner(kFALSE);;
969
970
971 /********************************/
972 /* Filtering out */
973 /********************************/
974 if (partEvent2==0x0)//in case fBufferSize == 0 and pointers are created do not eneter
975 {
976 partEvent2 = new AliAOD();
977 partEvent2->SetOwner(kTRUE);
978 }
979
980 FilterOut(partEvent1, partEvent2, rawpartEvent);
981
982 for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
983 {
984 if ( (j%fDisplayMixingInfo) == 0)
985 Info("ProcessParticlesNonIdentAnal",
986 "Mixing particle %d from current event with particles from current event",j);
987
988 part1= partEvent1->GetParticle(j);
989
990
991 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
992 fParticleMonitorFunctions[ii]->Process(part1);
993
994 if (fNParticleFunctions == 0) continue;
995
996 /***************************************/
997 /****** filling numerators ********/
998 /****** (particles from event2) ********/
999 /***************************************/
1000
1001 for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++) //partEvent1 and partEvent2 are particles from the same event but separated to two groups
1002 {
1003 part2= partEvent2->GetParticle(k);
1004 if (part1->GetUID() == part2->GetUID()) continue;//this is the same particle but with different PID
1005 partpair->SetParticles(part1,part2);
1006
1007
1008 if(fPairCut->PassPairProp(partpair) ) //check pair cut
1009 { //do not meets crietria of the pair cut
1010 continue;
1011 }
1012 else
1013 {//meets criteria of the pair cut
1014 for(ii = 0;ii<fNParticleFunctions;ii++)
1015 fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
1016 }
1017 }
1018
1019 if ( fBufferSize == 0) continue;//do not mix diff histograms
1020 /***************************************/
1021 /***** Filling denominators *********/
1022 /***************************************/
1023 fPartBuffer->ResetIter();
1024
1025 Int_t nmonitor = 0;
1026
1027 while ( (partEvent3 = fPartBuffer->Next() ) != 0x0)
1028 {
1029
1030 if ( (j%fDisplayMixingInfo) == 0)
1031 Info("ProcessParticlesNonIdentAnal",
1032 "Mixing particle %d from current event with particles from event%d",j,-(++nmonitor));
1033
1034 for (Int_t k = 0; k < partEvent3->GetNumberOfParticles() ; k++)
1035 {
1036 part2= partEvent3->GetParticle(k);
1037 partpair->SetParticles(part1,part2);
1038
1039
1040 if(fPairCut->PassPairProp(partpair) ) //check pair cut
1041 { //do not meets crietria of the pair cut
1042 continue;
1043 }
1044 else
1045 {//meets criteria of the pair cut
1046 for(ii = 0;ii<fNParticleFunctions;ii++)
1047 {
1048 fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
1049 }
1050 }
1051 }// for particles event2
1052 }//while event2
1053 }//for over particles in event1
1054
1055 delete fPartBuffer->Push(partEvent2);
1056
1057 return 0;
e92ecbdf 1058}
1059/*************************************************************************************/
090e46d6 1060Int_t AliHBTAnalysis::ProcessRecNonId(AliAOD* aodrec, AliAOD* /*aodsim*/)
e92ecbdf 1061{
090e46d6 1062//Analyzes both reconstructed and simulated data
1063 if (aodrec == 0x0)
1064 {
1065 return 1;
1066 }
1067
1068 AliVAODParticle * track1, * track2;
1069
1070 static AliAOD aodrec1;
1071
1072 AliAOD * trackEvent1=&aodrec1;//Particle that passes first particle cut, this event
1073 AliAOD * trackEvent2=0x0;//Particle that passes second particle cut, this event
1074 AliAOD * trackEvent3=0x0;//Particle that passes second particle cut, events from buffer
1075
1076 AliAOD* rawtrackEvent = aodrec;//this we get from Reader
1077
1078 static AliHBTPair tpair;
1079
1080 AliHBTPair* trackpair = &tpair;
1081
1082 register UInt_t ii;
1083
1084
1085 trackEvent1 = new AliAOD();
1086 trackEvent1->SetOwner(kFALSE);
1087
1088 /********************************/
1089 /* Filtering out */
1090 /********************************/
1091 if ( trackEvent2==0x0 )//in case fBufferSize == 0 and pointers are created do not eneter
1092 {
1093 trackEvent2 = new AliAOD();
1094 trackEvent2->SetOwner(kTRUE);
1095 }
1096
1097 FilterOut(trackEvent1, trackEvent2, rawtrackEvent);
1098
1099 for (Int_t j = 0; j<trackEvent1->GetNumberOfParticles() ; j++)
1100 {
1101 if ( (j%fDisplayMixingInfo) == 0)
1102 Info("ProcessTracksNonIdentAnal",
1103 "Mixing particle %d from current event with particles from current event",j);
1104
1105 track1= trackEvent1->GetParticle(j);
1106
1107
1108 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
1109 fTrackMonitorFunctions[ii]->Process(track1);
1110
1111 if (fNTrackFunctions == 0x0) continue;
1112
1113 /***************************************/
1114 /****** filling numerators ********/
1115 /****** (particles from event2) ********/
1116 /***************************************/
1117
1118 for (Int_t k = 0; k < trackEvent2->GetNumberOfParticles() ; k++) //partEvent1 and partEvent2 are particles from the same event but separated to two groups
1119 {
1120 track2= trackEvent2->GetParticle(k);
1121 if (track1->GetUID() == track2->GetUID()) continue;//this is the same particle but with different PID
1122 trackpair->SetParticles(track1,track2);
1123
1124
1125 if( fPairCut->PassPairProp(trackpair)) //check pair cut
1126 { //do not meets crietria of the pair cut
1127 continue;
1128 }
1129 else
1130 {//meets criteria of the pair cut
1131 UInt_t ii;
1132 for(ii = 0;ii<fNTrackFunctions;ii++)
1133 fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
1134 }
1135 }
1136
1137 if ( fBufferSize == 0) continue;//do not mix diff histograms
1138 /***************************************/
1139 /***** Filling denominators *********/
1140 /***************************************/
1141 fTrackBuffer->ResetIter();
1142
1143 Int_t nmonitor = 0;
1144
1145 while ( (trackEvent3 = fTrackBuffer->Next() ) != 0x0)
1146 {
1147 if ( (j%fDisplayMixingInfo) == 0)
1148 Info("ProcessTracksNonIdentAnal",
1149 "Mixing particle %d from current event with particles from event%d",j,-(++nmonitor));
1150
1151 for (Int_t k = 0; k < trackEvent3->GetNumberOfParticles() ; k++)
1152 {
1153 track2= trackEvent3->GetParticle(k);
1154 trackpair->SetParticles(track1,track2);
1155
1156 if( fPairCut->PassPairProp(trackpair)) //check pair cut
1157 { //do not meets crietria of the pair cut
1158 continue;
1159 }
1160 else
1161 {//meets criteria of the pair cut
1162 for(ii = 0;ii<fNTrackFunctions;ii++)
1163 fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
1164 }
1165 }// for particles event2
1166 }//while event2
1167 }//for over particles in event1
1168
1169 delete fTrackBuffer->Push(trackEvent2);
1170
1171 return 0;
e92ecbdf 1172}
1173/*************************************************************************************/
090e46d6 1174
1b446896 1175void AliHBTAnalysis::Process(Option_t* option)
1176{
1177 //default option = "TracksAndParticles"
1178 //Main method of the HBT Analysis Package
1179 //It triggers reading with the global cut (default is an empty cut)
1180 //Than it checks options and data which are read
1181 //if everything is OK, then it calls one of the looping methods
1182 //depending on tfReaderhe option
1183 //These methods differs on what thay are looping on
1184 //
1185 // METHOD OPTION
1186 //--------------------------------------------------------------------
1187 //ProcessTracksAndParticles - "TracksAndParticles"
1188 // DEFAULT
1189 // it loops over both, tracks(reconstructed) and particles(simulated)
1190 // all function gethered in all 3 lists are called for each (double)pair
1191 //
1192 //ProcessTracks - "Tracks"
1193 // it loops only on tracks(reconstructed),
1194 // functions ONLY from fTrackFunctions list are called
1195 //
1196 //ProcessParticles - "Particles"
1197 // it loops only on particles(simulated),
1198 // functions ONLY from fParticleAndTrackFunctions list are called
1199 //
1200 //
1201 if (!fReader)
1202 {
01725374 1203 Error("Process","The reader is not set");
1b446896 1204 return;
1205 }
1206
1b446896 1207 const char *oT = strstr(option,"Tracks");
1208 const char *oP = strstr(option,"Particles");
1209
dc2c3f36 1210 Bool_t nonid = IsNonIdentAnalysis();
1211
bed069a4 1212 Init();
1213
1b446896 1214 if(oT && oP)
1215 {
dc2c3f36 1216 if (nonid) ProcessTracksAndParticlesNonIdentAnal();
1217 else ProcessTracksAndParticles();
1b446896 1218 return;
1219 }
1220
1221 if(oT)
1222 {
dc2c3f36 1223 if (nonid) ProcessTracksNonIdentAnal();
1224 else ProcessTracks();
1b446896 1225 return;
1226 }
1227
1228 if(oP)
1229 {
dc2c3f36 1230 if (nonid) ProcessParticlesNonIdentAnal();
1231 else ProcessParticles();
1b446896 1232 return;
1233 }
1234
1235}
1b446896 1236/*************************************************************************************/
491d1b5d 1237
1b446896 1238void AliHBTAnalysis::ProcessTracksAndParticles()
1239{
9616170a 1240//Makes analysis for both tracks and particles
1241//mainly for resolution study and analysies with weighting algirithms
1b446896 1242//In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
1243//the loops are splited
1244
66d1d1a4 1245// cut on particles only -- why?
1246// - PID: when we make resolution analysis we want to take only tracks with correct PID
1247// We need cut on tracks because there are data characteristic to
1b446896 1248
78d7c6d3 1249 AliAOD * trackEvent, *partEvent;
81b7b887 1250
bed069a4 1251 fReader->Rewind();
bed069a4 1252 while (fReader->Next() == kFALSE)
1b446896 1253 {
e92ecbdf 1254 partEvent = fReader->GetEventSim();
78d7c6d3 1255 trackEvent = fReader->GetEventRec();
e92ecbdf 1256 ProcessRecAndSim(trackEvent,partEvent);
1b446896 1257
bed069a4 1258 }//while (fReader->Next() == kFALSE)
4cb6e8f7 1259
1b446896 1260}
1261/*************************************************************************************/
1262
1263void AliHBTAnalysis::ProcessTracks()
1264{
2dc7203b 1265//In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
1b446896 1266//the loops are splited
78d7c6d3 1267 AliAOD * trackEvent;
bed069a4 1268 fReader->Rewind();
bed069a4 1269 while (fReader->Next() == kFALSE)
1b446896 1270 {
78d7c6d3 1271 trackEvent = fReader->GetEventRec();
090e46d6 1272 ProcessRec(trackEvent,0x0);
bed069a4 1273 }//while (fReader->Next() == kFALSE)
1b446896 1274}
1275
1276/*************************************************************************************/
491d1b5d 1277
1b446896 1278void AliHBTAnalysis::ProcessParticles()
1279{
81b7b887 1280//In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
1b446896 1281//the loops are splited
78d7c6d3 1282 AliAOD * partEvent;
bed069a4 1283 fReader->Rewind();
bed069a4 1284 while (fReader->Next() == kFALSE)
1b446896 1285 {
78d7c6d3 1286 partEvent = fReader->GetEventSim();
090e46d6 1287 ProcessSim(0x0,partEvent);
bed069a4 1288 }//while (fReader->Next() == kFALSE)
1b446896 1289}
1b446896 1290/*************************************************************************************/
1291
491d1b5d 1292void AliHBTAnalysis::WriteFunctions()
1b446896 1293{
81b7b887 1294//Calls Write for all defined functions in analysis
1295//== writes all results
8fba7c63 1296 TFile* oututfile = 0x0;
1297 if (fOutputFileName)
1298 {
1299 oututfile = TFile::Open(*fOutputFileName,"update");
1300 }
1b446896 1301 UInt_t ii;
1302 for(ii = 0;ii<fNParticleFunctions;ii++)
90520373 1303 {
78d7c6d3 1304 if (AliVAODParticle::GetDebug()>5)
90520373 1305 {
1306 Info("WriteFunctions","Writing ParticleFunction %#x",fParticleFunctions[ii]);
1307 Info("WriteFunctions","Writing ParticleFunction %s",fParticleFunctions[ii]->Name());
1308 }
1309 fParticleFunctions[ii]->Write();
1310 }
1b446896 1311
1312 for(ii = 0;ii<fNTrackFunctions;ii++)
90520373 1313 {
78d7c6d3 1314 if (AliVAODParticle::GetDebug()>5)
90520373 1315 {
1316 Info("WriteFunctions","Writing TrackFunction %#x",fTrackFunctions[ii]);
1317 Info("WriteFunctions","Writing TrackFunction %s",fTrackFunctions[ii]->Name());
1318 }
1319 fTrackFunctions[ii]->Write();
1320 }
1b446896 1321
1322 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
90520373 1323 {
78d7c6d3 1324 if (AliVAODParticle::GetDebug()>5)
90520373 1325 {
1326 Info("WriteFunctions","Writing ParticleAndTrackFunction %#x",fParticleAndTrackFunctions[ii]);
1327 Info("WriteFunctions","Writing ParticleAndTrackFunction %s",fParticleAndTrackFunctions[ii]->Name());
1328 }
1329 fParticleAndTrackFunctions[ii]->Write();
1330 }
5c58441a 1331
1332 for(ii = 0;ii<fNParticleMonitorFunctions;ii++)
90520373 1333 {
78d7c6d3 1334 if (AliVAODParticle::GetDebug()>5)
90520373 1335 {
1336 Info("WriteFunctions","Writing ParticleMonitorFunction %#x",fParticleMonitorFunctions[ii]);
1337 Info("WriteFunctions","Writing ParticleMonitorFunction %s",fParticleMonitorFunctions[ii]->Name());
1338 }
1339 fParticleMonitorFunctions[ii]->Write();
1340 }
5c58441a 1341
1342 for(ii = 0;ii<fNTrackMonitorFunctions;ii++)
90520373 1343 {
78d7c6d3 1344 if (AliVAODParticle::GetDebug()>5)
90520373 1345 {
1346 Info("WriteFunctions","Writing TrackMonitorFunction %#x",fTrackMonitorFunctions[ii]);
1347 Info("WriteFunctions","Writing TrackMonitorFunction %s",fTrackMonitorFunctions[ii]->Name());
1348 }
1349 fTrackMonitorFunctions[ii]->Write();
1350 }
5c58441a 1351
1352 for(ii = 0;ii<fNParticleAndTrackMonitorFunctions;ii++)
90520373 1353 {
78d7c6d3 1354 if (AliVAODParticle::GetDebug()>5)
90520373 1355 {
1356 Info("WriteFunctions","Writing ParticleAndTrackMonitorFunction %#x",fParticleAndTrackMonitorFunctions[ii]);
1357 Info("WriteFunctions","Writing ParticleAndTrackMonitorFunction %s",fParticleAndTrackMonitorFunctions[ii]->Name());
1358 }
1359 fParticleAndTrackMonitorFunctions[ii]->Write();
1360 }
8fba7c63 1361 delete oututfile;
1362}
1363/*************************************************************************************/
1364
1365void AliHBTAnalysis::SetOutputFileName(const char* fname)
1366{
1367 //Sets fiele name where to dump results,
1368 //if not specified reults are written to gDirectory
1369 if (fname == 0x0)
1370 {
1371 delete fOutputFileName;
1372 fOutputFileName = 0x0;
1373 return;
1374 }
1375 if ( strcmp(fname,"") == 0 )
1376 {
1377 delete fOutputFileName;
1378 fOutputFileName = 0x0;
1379 return;
1380 }
1381 if (fOutputFileName == 0x0) fOutputFileName = new TString(fname);
1382 else *fOutputFileName = fname;
1b446896 1383}
1384/*************************************************************************************/
491d1b5d 1385
78d7c6d3 1386void AliHBTAnalysis::SetGlobalPairCut(AliAODPairCut* cut)
1b446896 1387{
81b7b887 1388//Sets the global cut
1b446896 1389 if (cut == 0x0)
1390 {
1391 Error("AliHBTAnalysis::SetGlobalPairCut","Pointer is NULL. Ignoring");
1392 }
1393 delete fPairCut;
78d7c6d3 1394 fPairCut = (AliAODPairCut*)cut->Clone();
1b446896 1395}
1396
1397/*************************************************************************************/
491d1b5d 1398
27b3fe5d 1399void AliHBTAnalysis::AddTrackFunction(AliHBTOnePairFctn* f)
1b446896 1400{
81b7b887 1401//Adds track function
1b446896 1402 if (f == 0x0) return;
1403 if (fNTrackFunctions == fgkFctnArraySize)
1404 {
1405 Error("AliHBTAnalysis::AddTrackFunction","Can not add this function, not enough place in the array.");
1406 }
1407 fTrackFunctions[fNTrackFunctions] = f;
1408 fNTrackFunctions++;
1409}
491d1b5d 1410/*************************************************************************************/
1411
27b3fe5d 1412void AliHBTAnalysis::AddParticleFunction(AliHBTOnePairFctn* f)
1b446896 1413{
81b7b887 1414//adds particle function
1b446896 1415 if (f == 0x0) return;
1416
1417 if (fNParticleFunctions == fgkFctnArraySize)
1418 {
1419 Error("AliHBTAnalysis::AddParticleFunction","Can not add this function, not enough place in the array.");
1420 }
1421 fParticleFunctions[fNParticleFunctions] = f;
1422 fNParticleFunctions++;
1b446896 1423}
5c58441a 1424/*************************************************************************************/
1425
27b3fe5d 1426void AliHBTAnalysis::AddParticleAndTrackFunction(AliHBTTwoPairFctn* f)
1b446896 1427{
81b7b887 1428//add resolution function
1b446896 1429 if (f == 0x0) return;
1430 if (fNParticleAndTrackFunctions == fgkFctnArraySize)
1431 {
1432 Error("AliHBTAnalysis::AddParticleAndTrackFunction","Can not add this function, not enough place in the array.");
1433 }
1434 fParticleAndTrackFunctions[fNParticleAndTrackFunctions] = f;
1435 fNParticleAndTrackFunctions++;
1436}
5c58441a 1437/*************************************************************************************/
1438
1439void AliHBTAnalysis::AddParticleMonitorFunction(AliHBTMonOneParticleFctn* f)
1440{
81b7b887 1441//add particle monitoring function
5c58441a 1442 if (f == 0x0) return;
1443
1444 if (fNParticleMonitorFunctions == fgkFctnArraySize)
1445 {
1446 Error("AliHBTAnalysis::AddParticleMonitorFunction","Can not add this function, not enough place in the array.");
1447 }
1448 fParticleMonitorFunctions[fNParticleMonitorFunctions] = f;
1449 fNParticleMonitorFunctions++;
1450}
1451/*************************************************************************************/
1b446896 1452
5c58441a 1453void AliHBTAnalysis::AddTrackMonitorFunction(AliHBTMonOneParticleFctn* f)
1454{
81b7b887 1455//add track monitoring function
5c58441a 1456 if (f == 0x0) return;
1b446896 1457
5c58441a 1458 if (fNTrackMonitorFunctions == fgkFctnArraySize)
1459 {
1460 Error("AliHBTAnalysis::AddTrackMonitorFunction","Can not add this function, not enough place in the array.");
1461 }
1462 fTrackMonitorFunctions[fNTrackMonitorFunctions] = f;
1463 fNTrackMonitorFunctions++;
1464}
1b446896 1465/*************************************************************************************/
1466
5c58441a 1467void AliHBTAnalysis::AddParticleAndTrackMonitorFunction(AliHBTMonTwoParticleFctn* f)
1468{
81b7b887 1469//add resolution monitoring function
5c58441a 1470 if (f == 0x0) return;
1471 if (fNParticleAndTrackMonitorFunctions == fgkFctnArraySize)
1472 {
1473 Error("AliHBTAnalysis::AddParticleAndTrackMonitorFunction","Can not add this function, not enough place in the array.");
1474 }
1475 fParticleAndTrackMonitorFunctions[fNParticleAndTrackMonitorFunctions] = f;
1476 fNParticleAndTrackMonitorFunctions++;
1477}
1478
1b446896 1479
5c58441a 1480/*************************************************************************************/
1b446896 1481/*************************************************************************************/
491d1b5d 1482
1b446896 1483Bool_t AliHBTAnalysis::RunCoherencyCheck()
1484{
1485 //Checks if both HBTRuns are similar
1486 //return true if error found
1487 //if they seem to be OK return false
78d7c6d3 1488/*
1489
1b446896 1490 Int_t i;
81b7b887 1491 Info("RunCoherencyCheck","Checking HBT Runs Coherency");
1492
78d7c6d3 1493//When we use non-buffering reader this is a big waste of time -> We need to read all data to check it
1494//and reader is implemented safe in this case anyway
1495// Info("RunCoherencyCheck","Number of events ...");
1496// if (fReader->GetNumberOfPartEvents() == fReader->GetNumberOfTrackEvents() ) //check whether there is the same number of events
1497// {
1498// Info("RunCoherencyCheck","OK. %d found\n",fReader->GetNumberOfTrackEvents());
1499// }
1500// else
1501// { //if not the same - ERROR
1502// Error("RunCoherencyCheck",
1503// "Number of simulated events (%d) is not equal to number of reconstructed events(%d)",
1504// fReader->GetNumberOfPartEvents(),fReader->GetNumberOfTrackEvents());
1505// return kTRUE;
1506// }
1b446896 1507
81b7b887 1508 Info("RunCoherencyCheck","Checking number of Particles AND Particles Types in each event ...");
1b446896 1509
78d7c6d3 1510 AliAOD *partEvent;
1511 AliAOD *trackEvent;
1b446896 1512 for( i = 0; i<fReader->GetNumberOfTrackEvents();i++)
1513 {
78d7c6d3 1514 partEvent= fReader->GetEventSim(i); //gets the "ith" event
1515 trackEvent = fReader->GetEventRec(i);
1b446896 1516
1517 if ( (partEvent == 0x0) && (partEvent == 0x0) ) continue;
1518 if ( (partEvent == 0x0) || (partEvent == 0x0) )
1519 {
78d7c6d3 1520 Error("RunCoherencyCheck",
1b446896 1521 "One event is NULL and the other one not. Event Number %d",i);
1522 return kTRUE;
1523 }
1524
1525 if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
1526 {
78d7c6d3 1527 Error("RunCoherencyCheck",
1b446896 1528 "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
1529 i,partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
1530 return kTRUE;
1531 }
1532 else
1533 for (Int_t j = 0; j<partEvent->GetNumberOfParticles(); j++)
1534 {
1535 if( partEvent->GetParticle(j)->GetPdgCode() != trackEvent->GetParticle(j)->GetPdgCode() )
1536 {
78d7c6d3 1537 Error("RunCoherencyCheck",
1b446896 1538 "Event %d: Particle %d: PID of simulated particle (%d) not the same of reconstructed track (%d)",
1539 i,j, partEvent->GetParticle(j)->GetPdgCode(),trackEvent->GetParticle(j)->GetPdgCode() );
1540 return kTRUE;
1541
1542 }
1543 }
1544 }
81b7b887 1545 Info("RunCoherencyCheck"," Done");
1546 Info("RunCoherencyCheck"," Everything looks OK");
78d7c6d3 1547*/
1b446896 1548 return kFALSE;
1549}
1550
dc2c3f36 1551/*************************************************************************************/
1552
1553void AliHBTAnalysis::ProcessTracksAndParticlesNonIdentAnal()
1554{
81b7b887 1555//Performs analysis for both, tracks and particles
090e46d6 1556 AliAOD* rawtrackEvent, * rawpartEvent;
bed069a4 1557 fReader->Rewind();
090e46d6 1558
81b7b887 1559 Info("ProcessTracksAndParticlesNonIdentAnal","**************************************");
1560 Info("ProcessTracksAndParticlesNonIdentAnal","***** NON IDENT MODE ****************");
1561 Info("ProcessTracksAndParticlesNonIdentAnal","**************************************");
090e46d6 1562
bed069a4 1563 for (Int_t i = 0;;i++)//infinite loop
dc2c3f36 1564 {
bed069a4 1565 if (fReader->Next()) break; //end when no more events available
1566
78d7c6d3 1567 rawpartEvent = fReader->GetEventSim();
1568 rawtrackEvent = fReader->GetEventRec();
bed069a4 1569
090e46d6 1570 ProcessRecAndSimNonId(rawtrackEvent,rawpartEvent);
dc2c3f36 1571 }//end of loop over events (1)
dc2c3f36 1572}
1b446896 1573/*************************************************************************************/
1574
dc2c3f36 1575void AliHBTAnalysis::ProcessTracksNonIdentAnal()
1576{
2dc7203b 1577//Process Tracks only with non identical mode
78d7c6d3 1578 AliAOD * rawtrackEvent;
bed069a4 1579 fReader->Rewind();
1580
81b7b887 1581 Info("ProcessTracksNonIdentAnal","**************************************");
1582 Info("ProcessTracksNonIdentAnal","***** NON IDENT MODE ****************");
1583 Info("ProcessTracksNonIdentAnal","**************************************");
1584
bed069a4 1585 for (Int_t i = 0;;i++)//infinite loop
dc2c3f36 1586 {
bed069a4 1587 if (fReader->Next()) break; //end when no more events available
78d7c6d3 1588 rawtrackEvent = fReader->GetEventRec();
090e46d6 1589 ProcessRecNonId(rawtrackEvent,0x0);
dc2c3f36 1590 }//end of loop over events (1)
dc2c3f36 1591}
1592/*************************************************************************************/
1593
1594void AliHBTAnalysis::ProcessParticlesNonIdentAnal()
1595{
2dc7203b 1596//process paricles only with non identical mode
78d7c6d3 1597 AliAOD * rawpartEvent = 0x0;
bed069a4 1598 fReader->Rewind();
dc2c3f36 1599
81b7b887 1600 Info("ProcessParticlesNonIdentAnal","**************************************");
1601 Info("ProcessParticlesNonIdentAnal","***** NON IDENT MODE ****************");
1602 Info("ProcessParticlesNonIdentAnal","**************************************");
dc2c3f36 1603
bed069a4 1604 for (Int_t i = 0;;i++)//infinite loop
dc2c3f36 1605 {
bed069a4 1606 if (fReader->Next()) break; //end when no more events available
1607
78d7c6d3 1608 rawpartEvent = fReader->GetEventSim();
090e46d6 1609 ProcessSimNonId(0x0,rawpartEvent);
dc2c3f36 1610 }//end of loop over events (1)
dc2c3f36 1611}
1612
1613/*************************************************************************************/
78d7c6d3 1614void AliHBTAnalysis::FilterOut(AliAOD* outpart1, AliAOD* outpart2, AliAOD* inpart,
1615 AliAOD* outtrack1, AliAOD* outtrack2, AliAOD* intrack) const
dc2c3f36 1616{
1617 //Puts particles accepted as a first particle by global cut in out1
1618 //and as a second particle in out2
1619
78d7c6d3 1620 AliVAODParticle* part, *track;
dc2c3f36 1621
1622 outpart1->Reset();
1623 outpart2->Reset();
1624 outtrack1->Reset();
1625 outtrack2->Reset();
1626
dc2c3f36 1627 Bool_t in1, in2;
1628
1629 for (Int_t i = 0; i < inpart->GetNumberOfParticles(); i++)
1630 {
1631 in1 = in2 = kTRUE;
1632 part = inpart->GetParticle(i);
1633 track = intrack->GetParticle(i);
1634
17d74b37 1635 if ( ((this->*fkPass1)(part,track)) ) in1 = kFALSE; //if part is rejected by cut1, in1 is false
1636 if ( ((this->*fkPass2)(part,track)) ) in2 = kFALSE; //if part is rejected by cut2, in2 is false
dc2c3f36 1637
1638 if (gDebug)//to be removed in real analysis
1639 if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
1640 {
1641 //Particle accpted by both cuts
1642 Error("FilterOut","Particle accepted by both cuts");
1643 continue;
1644 }
1645
1646 if (in1)
1647 {
1648 outpart1->AddParticle(part);
1649 outtrack1->AddParticle(track);
1650 continue;
1651 }
1652
1653 if (in2)
1654 {
d50bb49a 1655 outpart2->AddParticle(part);
1656 outtrack2->AddParticle(track);
dc2c3f36 1657 continue;
1658 }
1659 }
dc2c3f36 1660}
1b446896 1661/*************************************************************************************/
81b7b887 1662
78d7c6d3 1663void AliHBTAnalysis::FilterOut(AliAOD* out1, AliAOD* out2, AliAOD* in) const
dc2c3f36 1664{
1665 //Puts particles accepted as a first particle by global cut in out1
1666 //and as a second particle in out2
78d7c6d3 1667 AliVAODParticle* part;
dc2c3f36 1668
1669 out1->Reset();
1670 out2->Reset();
1671
78d7c6d3 1672 AliAODParticleCut *cut1 = fPairCut->GetFirstPartCut();
1673 AliAODParticleCut *cut2 = fPairCut->GetSecondPartCut();
dc2c3f36 1674
1675 Bool_t in1, in2;
1676
1677 for (Int_t i = 0; i < in->GetNumberOfParticles(); i++)
1678 {
1679 in1 = in2 = kTRUE;
1680 part = in->GetParticle(i);
1681
cea0a066 1682 if ( cut1->Rejected(part) ) in1 = kFALSE; //if part is rejected by cut1, in1 is false
1683 if ( cut2->Rejected(part) ) in2 = kFALSE; //if part is rejected by cut2, in2 is false
dc2c3f36 1684
1685 if (gDebug)//to be removed in real analysis
1686 if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
1687 {
1688 //Particle accpted by both cuts
1689 Error("FilterOut","Particle accepted by both cuts");
1690 continue;
1691 }
1b446896 1692
dc2c3f36 1693 if (in1)
1694 {
1695 out1->AddParticle(part);
1696 continue;
1697 }
1698
1699 if (in2)
1700 {
d50bb49a 1701 out2->AddParticle(part);
dc2c3f36 1702 continue;
1703 }
1704 }
1705}
1706/*************************************************************************************/
1b446896 1707
dc2c3f36 1708Bool_t AliHBTAnalysis::IsNonIdentAnalysis()
1709{
1710 //checks if it is possible to use special analysis for non identical particles
1711 //it means - in global pair cut first particle id is different than second one
1712 //and both are different from 0
1713 //in the future is possible to perform more sophisticated check
1714 //if cuts have excluding requirements
1715
5c58441a 1716 if (fPairCut->IsEmpty())
1717 return kFALSE;
1718
1719 if (fPairCut->GetFirstPartCut()->IsEmpty())
1720 return kFALSE;
1721
1722 if (fPairCut->GetSecondPartCut()->IsEmpty())
1723 return kFALSE;
dc2c3f36 1724
1725 Int_t id1 = fPairCut->GetFirstPartCut()->GetPID();
1726 Int_t id2 = fPairCut->GetSecondPartCut()->GetPID();
dc2c3f36 1727
5c58441a 1728 if ( (id1==0) || (id2==0) || (id1==id2) )
1729 return kFALSE;
1730
dc2c3f36 1731 return kTRUE;
1732}
66d1d1a4 1733/*************************************************************************************/
9616170a 1734
5994509d 1735void AliHBTAnalysis::SetApparentVertex(Double_t x, Double_t y, Double_t z)
1736{
1737 //Sets apparent vertex
1738 // All events have to be moved to the same vertex position in order to
1739 // to be able to comare any space positions (f.g. anti-merging)
1740 // This method defines this position
1741
1742 fVertexX = x;
1743 fVertexY = y;
1744 fVertexZ = z;
1745}
1746/*************************************************************************************/
1747
9616170a 1748void AliHBTAnalysis::PressAnyKey()
17d74b37 1749{
1750 //small utility function that helps to make comfortable macros
9616170a 1751 char c;
1752 int nread = -1;
1753 fcntl(0, F_SETFL, O_NONBLOCK);
1754 ::Info("","Press Any Key to continue ...");
1755 while (nread<1)
1756 {
1757 nread = read(0, &c, 1);
1758 gSystem->ProcessEvents();
1759 }
1760}
1761
66d1d1a4 1762/*************************************************************************************/
9616170a 1763