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