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