Pair cuts and all functonality connected with them moved to AliAnalysis
[u/mrichter/AliRoot.git] / HBTAN / AliHBTAnalysis.cxx
1 #include "AliHBTAnalysis.h"
2 //_________________________________________________________
3 ////////////////////////////////////////////////////////////////////////////
4 //
5 // class AliHBTAnalysis
6 //
7 // Central Object Of HBTAnalyser: 
8 // This class performs main looping within HBT Analysis
9 // User must plug a reader of Type AliReader
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
22
23 #include <TSystem.h>
24 #include <TFile.h>
25
26 #include "AliAOD.h"
27 #include "AliAODParticle.h"
28 #include "AliAODPairCut.h"
29 #include "AliEventCut.h"
30
31 #include "AliEventBuffer.h"
32
33 #include "AliReader.h"
34 #include "AliHBTPair.h"
35 #include "AliHBTFunction.h"
36 #include "AliHBTMonitorFunction.h"
37  
38
39 ClassImp(AliHBTAnalysis)
40
41 const UInt_t AliHBTAnalysis::fgkFctnArraySize = 100;
42 const UInt_t AliHBTAnalysis::fgkDefaultMixingInfo = 1000;
43 const Int_t  AliHBTAnalysis::fgkDefaultBufferSize = 5;
44
45 AliHBTAnalysis::AliHBTAnalysis():
46   fProcEvent(0x0),
47   fReader(0x0),
48   fNTrackFunctions(0),
49   fNParticleFunctions(0),
50   fNParticleAndTrackFunctions(0),
51   fNTrackMonitorFunctions(0),
52   fNParticleMonitorFunctions(0), 
53   fNParticleAndTrackMonitorFunctions(0),
54   fTrackFunctions ( new AliHBTOnePairFctn* [fgkFctnArraySize]),
55   fParticleFunctions ( new AliHBTOnePairFctn* [fgkFctnArraySize]),
56   fParticleAndTrackFunctions ( new AliHBTTwoPairFctn* [fgkFctnArraySize]),
57   fParticleMonitorFunctions ( new AliHBTMonOneParticleFctn* [fgkFctnArraySize]),    
58   fTrackMonitorFunctions ( new AliHBTMonOneParticleFctn* [fgkFctnArraySize]),    
59   fParticleAndTrackMonitorFunctions ( new AliHBTMonTwoParticleFctn* [fgkFctnArraySize]),    
60   fBkgEventCut(0x0),
61   fPartBuffer(0x0),
62   fTrackBuffer(0x0),
63   fBufferSize(2),
64   fDisplayMixingInfo(fgkDefaultMixingInfo),
65   fIsOwner(kFALSE),
66   fProcessOption(kSimulatedAndReconstructed),
67   fNoCorrfctns(kFALSE),
68   fOutputFileName(0x0),
69   fVertexX(0.0),
70   fVertexY(0.0),
71   fVertexZ(0.0)
72  {
73    //default constructor
74    
75  }
76 /*************************************************************************************/ 
77
78 AliHBTAnalysis::AliHBTAnalysis(const AliHBTAnalysis& in):
79   AliAnalysis(in),
80   fProcEvent(0x0),
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),
94   fBkgEventCut(0x0),
95   fPartBuffer(0x0),
96   fTrackBuffer(0x0),
97   fBufferSize(fgkDefaultBufferSize),
98   fDisplayMixingInfo(fgkDefaultMixingInfo),
99   fIsOwner(kFALSE),
100   fProcessOption(kSimulatedAndReconstructed),
101   fNoCorrfctns(kFALSE),
102   fOutputFileName(0x0),
103   fVertexX(0.0),
104   fVertexY(0.0),
105   fVertexZ(0.0)
106  {
107 //copy constructor
108    Fatal("AliHBTAnalysis(const AliHBTAnalysis&)","Sensless");
109  }
110 /*************************************************************************************/ 
111 AliHBTAnalysis& AliHBTAnalysis::operator=(const AliHBTAnalysis& /*right*/)
112  {
113 //operator =
114    Fatal("AliHBTAnalysis(const AliHBTAnalysis&)","Sensless");
115    return *this;
116  }
117 /*************************************************************************************/ 
118 AliHBTAnalysis::~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
124     
125    if (fIsOwner)
126     {
127       if (AliVAODParticle::GetDebug()>5)Info("~AliHBTAnalysis","Is Owner: Attempting to delete functions");
128       DeleteFunctions();
129       if (AliVAODParticle::GetDebug()>5)Info("~AliHBTAnalysis","Delete functions done");
130     }
131    delete [] fTrackFunctions;
132    delete [] fParticleFunctions;
133    delete [] fParticleAndTrackFunctions;
134    
135    delete [] fParticleMonitorFunctions; 
136    delete [] fTrackMonitorFunctions; 
137    delete [] fParticleAndTrackMonitorFunctions;
138
139    delete fBkgEventCut;
140    delete fOutputFileName;
141  }
142
143 /*************************************************************************************/ 
144
145 Int_t AliHBTAnalysis::ProcessEvent(AliAOD* aodrec, AliAOD* aodsim)
146 {
147   //Processes one event
148   if (fProcEvent == 0x0)
149    {
150      Error("ProcessEvent","Analysis option not specified");
151      return 1;
152    }
153   if ( Pass(aodrec,aodsim) ) return 0;
154   
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;
175 }
176 /*************************************************************************************/ 
177
178 Int_t AliHBTAnalysis::Finish()
179 {
180 //Finishes analysis
181   WriteFunctions();
182   return 0;
183 }
184 /*************************************************************************************/ 
185
186 void AliHBTAnalysis::DeleteFunctions()
187 {
188  //Deletes all functions added to analysis
189  
190  UInt_t ii;
191  for(ii = 0;ii<fNParticleFunctions;ii++)
192   { 
193     if (AliVAODParticle::GetDebug()>5)
194      {
195        Info("DeleteFunctions","Deleting ParticleFunction %#x",fParticleFunctions[ii]);
196        Info("DeleteFunctions","Deleting ParticleFunction %s",fParticleFunctions[ii]->Name());
197      }
198     delete fParticleFunctions[ii];
199   } 
200  fNParticleFunctions = 0;
201                 
202  for(ii = 0;ii<fNTrackFunctions;ii++)
203   { 
204     if (AliVAODParticle::GetDebug()>5)
205      {
206        Info("DeleteFunctions","Deleting TrackFunction %#x",fTrackFunctions[ii]);
207        Info("DeleteFunctions","Deleting TrackFunction %s",fTrackFunctions[ii]->Name());
208      }
209     delete fTrackFunctions[ii];
210   }  
211  fNTrackFunctions = 0;
212  
213  for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
214   { 
215     if (AliVAODParticle::GetDebug()>5)
216      {
217        Info("DeleteFunctions","Deleting ParticleAndTrackFunction %#x",fParticleAndTrackFunctions[ii]);
218        Info("DeleteFunctions","Deleting ParticleAndTrackFunction %s",fParticleAndTrackFunctions[ii]->Name());
219      } 
220     delete fParticleAndTrackFunctions[ii];
221   }  
222  fNParticleAndTrackFunctions = 0;
223  
224  for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
225   { 
226     if (AliVAODParticle::GetDebug()>5)
227      {
228        Info("DeleteFunctions","Deleting ParticleMonitorFunction %#x",fParticleMonitorFunctions[ii]);
229        Info("DeleteFunctions","Deleting ParticleMonitorFunction %s",fParticleMonitorFunctions[ii]->Name());
230      }
231     delete fParticleMonitorFunctions[ii];
232   } 
233  fNParticleMonitorFunctions = 0;
234    
235  for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
236   { 
237     if (AliVAODParticle::GetDebug()>5)
238      {
239        Info("DeleteFunctions","Deleting TrackMonitorFunction %#x",fTrackMonitorFunctions[ii]);
240        Info("DeleteFunctions","Deleting TrackMonitorFunction %s",fTrackMonitorFunctions[ii]->Name());
241      }
242     delete fTrackMonitorFunctions[ii];
243   } 
244  fNTrackMonitorFunctions = 0;
245    
246  for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
247   { 
248     if (AliVAODParticle::GetDebug()>5)
249      {
250       Info("DeleteFunctions","Deleting ParticleAndTrackMonitorFunction %#x",fParticleAndTrackMonitorFunctions[ii]);
251       Info("DeleteFunctions","Deleting ParticleAndTrackMonitorFunction %s",fParticleAndTrackMonitorFunctions[ii]->Name());
252      }
253     delete fParticleAndTrackMonitorFunctions[ii];
254   } 
255  fNParticleAndTrackMonitorFunctions = 0;
256  
257 }
258 /*************************************************************************************/ 
259
260 Int_t AliHBTAnalysis::Init()
261 {
262 //Initializeation method
263 //calls Init for all functions
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   }
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  
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();
312    
313  return 0;
314 }
315 /*************************************************************************************/ 
316
317 void 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 /*************************************************************************************/ 
333 Int_t AliHBTAnalysis::ProcessRecAndSim(AliAOD* aodrec, AliAOD* aodsim)
334 {
335 //Does analysis for both tracks and particles
336 //mainly for resolution study and analysies with weighting algirithms
337   
338 // cut on particles only -- why?
339 // - PID: when we make resolution analysis we want to take only tracks with correct PID
340 // We need cut on tracks because there are data characteristic
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();
356   static AliHBTPair tpair;
357   static AliHBTPair ppair;
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 /*************************************************************************************/ 
515
516 Int_t AliHBTAnalysis::ProcessSim(AliAOD* /*aodrec*/, AliAOD* aodsim)
517 {
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;
638 }
639 /*************************************************************************************/ 
640 Int_t AliHBTAnalysis::ProcessRec(AliAOD* aodrec, AliAOD* /*aodsim*/)
641 {
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;
761 }
762 /*************************************************************************************/ 
763      
764 Int_t AliHBTAnalysis::ProcessRecAndSimNonId(AliAOD* aodrec, AliAOD* aodsim)
765 {
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;
927 }
928 /*************************************************************************************/ 
929 Int_t AliHBTAnalysis::ProcessSimNonId(AliAOD* /*aodrec*/, AliAOD* aodsim)
930 {
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;
1048 }
1049 /*************************************************************************************/ 
1050 Int_t AliHBTAnalysis::ProcessRecNonId(AliAOD* aodrec, AliAOD* /*aodsim*/)
1051 {
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;
1162 }
1163 /*************************************************************************************/ 
1164
1165 void 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   {
1193    Error("Process","The reader is not set");
1194    return;
1195   }
1196  
1197  const char *oT = strstr(option,"Tracks");
1198  const char *oP = strstr(option,"Particles");
1199  
1200  Bool_t nonid = IsNonIdentAnalysis();
1201  
1202  Init();
1203  
1204  if(oT && oP)
1205   { 
1206     if (nonid) ProcessTracksAndParticlesNonIdentAnal();
1207     else ProcessTracksAndParticles();
1208     return;
1209   }
1210
1211  if(oT)
1212   {
1213     if (nonid) ProcessTracksNonIdentAnal();
1214     else ProcessTracks();
1215     return;
1216   }
1217  
1218  if(oP)
1219   {
1220     if (nonid) ProcessParticlesNonIdentAnal();
1221     else ProcessParticles();
1222     return;
1223   }
1224  
1225 }
1226 /*************************************************************************************/ 
1227
1228 void AliHBTAnalysis::ProcessTracksAndParticles()
1229 {
1230 //Makes analysis for both tracks and particles
1231 //mainly for resolution study and analysies with weighting algirithms
1232 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
1233 //the loops are splited
1234   
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 
1238   
1239   AliAOD * trackEvent, *partEvent;
1240
1241   fReader->Rewind();
1242   while (fReader->Next() == kFALSE)
1243     {
1244       partEvent = fReader->GetEventSim();
1245       trackEvent = fReader->GetEventRec();
1246       ProcessRecAndSim(trackEvent,partEvent);
1247
1248     }//while (fReader->Next() == kFALSE)
1249     
1250
1251 /*************************************************************************************/
1252
1253 void AliHBTAnalysis::ProcessTracks()
1254 {
1255 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
1256 //the loops are splited
1257   AliAOD * trackEvent;
1258   fReader->Rewind();
1259   while (fReader->Next() == kFALSE)
1260     {
1261       trackEvent = fReader->GetEventRec();
1262       ProcessRec(trackEvent,0x0);
1263     }//while (fReader->Next() == kFALSE)
1264 }
1265
1266 /*************************************************************************************/
1267
1268 void AliHBTAnalysis::ProcessParticles()
1269 {
1270 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
1271 //the loops are splited
1272   AliAOD * partEvent;
1273   fReader->Rewind();
1274   while (fReader->Next() == kFALSE)
1275     {
1276       partEvent = fReader->GetEventSim();
1277       ProcessSim(0x0,partEvent);
1278     }//while (fReader->Next() == kFALSE)
1279 }
1280 /*************************************************************************************/
1281
1282 void AliHBTAnalysis::WriteFunctions()
1283 {
1284 //Calls Write for all defined functions in analysis
1285 //== writes all results
1286   TFile* oututfile = 0x0;
1287   if (fOutputFileName)
1288    {
1289      oututfile = TFile::Open(*fOutputFileName,"update");
1290    }
1291   UInt_t ii;
1292   for(ii = 0;ii<fNParticleFunctions;ii++)
1293    {
1294      if (AliVAODParticle::GetDebug()>5)
1295       {
1296         Info("WriteFunctions","Writing ParticleFunction %#x",fParticleFunctions[ii]);
1297         Info("WriteFunctions","Writing ParticleFunction %s",fParticleFunctions[ii]->Name());
1298       }
1299      fParticleFunctions[ii]->Write();
1300    }
1301                  
1302   for(ii = 0;ii<fNTrackFunctions;ii++)
1303    {
1304      if (AliVAODParticle::GetDebug()>5)
1305       {
1306         Info("WriteFunctions","Writing TrackFunction %#x",fTrackFunctions[ii]);
1307         Info("WriteFunctions","Writing TrackFunction %s",fTrackFunctions[ii]->Name());
1308       }
1309      fTrackFunctions[ii]->Write();
1310    }
1311                  
1312   for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
1313    {
1314      if (AliVAODParticle::GetDebug()>5)
1315       {
1316         Info("WriteFunctions","Writing ParticleAndTrackFunction %#x",fParticleAndTrackFunctions[ii]);
1317         Info("WriteFunctions","Writing ParticleAndTrackFunction %s",fParticleAndTrackFunctions[ii]->Name());
1318       } 
1319      fParticleAndTrackFunctions[ii]->Write();
1320    }
1321
1322   for(ii = 0;ii<fNParticleMonitorFunctions;ii++)
1323    {
1324      if (AliVAODParticle::GetDebug()>5)
1325       {
1326         Info("WriteFunctions","Writing ParticleMonitorFunction %#x",fParticleMonitorFunctions[ii]);
1327         Info("WriteFunctions","Writing ParticleMonitorFunction %s",fParticleMonitorFunctions[ii]->Name());
1328       }
1329      fParticleMonitorFunctions[ii]->Write();
1330    }
1331
1332   for(ii = 0;ii<fNTrackMonitorFunctions;ii++)
1333    {
1334      if (AliVAODParticle::GetDebug()>5)
1335       {
1336         Info("WriteFunctions","Writing TrackMonitorFunction %#x",fTrackMonitorFunctions[ii]);
1337         Info("WriteFunctions","Writing TrackMonitorFunction %s",fTrackMonitorFunctions[ii]->Name());
1338       }
1339      fTrackMonitorFunctions[ii]->Write();
1340    }
1341
1342   for(ii = 0;ii<fNParticleAndTrackMonitorFunctions;ii++)
1343    {
1344      if (AliVAODParticle::GetDebug()>5)
1345       {
1346        Info("WriteFunctions","Writing ParticleAndTrackMonitorFunction %#x",fParticleAndTrackMonitorFunctions[ii]);
1347        Info("WriteFunctions","Writing ParticleAndTrackMonitorFunction %s",fParticleAndTrackMonitorFunctions[ii]->Name());
1348       }
1349      fParticleAndTrackMonitorFunctions[ii]->Write();
1350    }
1351   delete oututfile; 
1352 }
1353 /*************************************************************************************/
1354
1355 void 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;
1373 }
1374 /*************************************************************************************/
1375
1376 void AliHBTAnalysis::SetGlobalPairCut(AliAODPairCut* cut)
1377 {
1378 //Sets the global cut
1379   if (cut == 0x0)
1380    {
1381      Error("AliHBTAnalysis::SetGlobalPairCut","Pointer is NULL. Ignoring");
1382    }
1383   delete fPairCut;
1384   fPairCut = (AliAODPairCut*)cut->Clone();
1385 }
1386
1387 /*************************************************************************************/
1388
1389 void AliHBTAnalysis::AddTrackFunction(AliHBTOnePairFctn* f)
1390 {
1391 //Adds track function
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 }
1400 /*************************************************************************************/ 
1401
1402 void AliHBTAnalysis::AddParticleFunction(AliHBTOnePairFctn* f)
1403 {
1404 //adds particle function
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++;
1413 }
1414 /*************************************************************************************/ 
1415
1416 void AliHBTAnalysis::AddParticleAndTrackFunction(AliHBTTwoPairFctn* f)
1417 {
1418 //add resolution function
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 }
1427 /*************************************************************************************/ 
1428
1429 void AliHBTAnalysis::AddParticleMonitorFunction(AliHBTMonOneParticleFctn* f)
1430 {
1431 //add particle monitoring function
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 /*************************************************************************************/ 
1442
1443 void AliHBTAnalysis::AddTrackMonitorFunction(AliHBTMonOneParticleFctn* f)
1444 {
1445 //add track monitoring function
1446   if (f == 0x0) return;
1447
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 }
1455 /*************************************************************************************/ 
1456
1457 void AliHBTAnalysis::AddParticleAndTrackMonitorFunction(AliHBTMonTwoParticleFctn* f)
1458 {
1459 //add resolution monitoring function
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
1469
1470 /*************************************************************************************/ 
1471 /*************************************************************************************/  
1472
1473 Bool_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
1478 /* 
1479
1480  Int_t i;  
1481  Info("RunCoherencyCheck","Checking HBT Runs Coherency");
1482
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 //  }
1497  
1498  Info("RunCoherencyCheck","Checking number of Particles AND Particles Types in each event ...");
1499  
1500  AliAOD *partEvent;
1501  AliAOD *trackEvent;
1502  for( i = 0; i<fReader->GetNumberOfTrackEvents();i++)
1503   {
1504     partEvent= fReader->GetEventSim(i); //gets the "ith" event 
1505     trackEvent = fReader->GetEventRec(i);
1506     
1507     if ( (partEvent == 0x0) && (partEvent == 0x0) ) continue;
1508     if ( (partEvent == 0x0) || (partEvent == 0x0) )
1509      {
1510        Error("RunCoherencyCheck",
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      {
1517        Error("RunCoherencyCheck",
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          {
1527            Error("RunCoherencyCheck",
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   }
1535  Info("RunCoherencyCheck","  Done");
1536  Info("RunCoherencyCheck","  Everything looks OK");
1537 */ 
1538  return kFALSE;
1539 }
1540
1541 /*************************************************************************************/  
1542  
1543 void AliHBTAnalysis::ProcessTracksAndParticlesNonIdentAnal()
1544 {
1545 //Performs analysis for both, tracks and particles
1546   AliAOD* rawtrackEvent, * rawpartEvent;
1547   fReader->Rewind();
1548  
1549   Info("ProcessTracksAndParticlesNonIdentAnal","**************************************");
1550   Info("ProcessTracksAndParticlesNonIdentAnal","*****  NON IDENT MODE ****************");
1551   Info("ProcessTracksAndParticlesNonIdentAnal","**************************************");
1552
1553   for (Int_t i = 0;;i++)//infinite loop
1554     {
1555       if (fReader->Next()) break; //end when no more events available
1556       
1557       rawpartEvent  = fReader->GetEventSim();
1558       rawtrackEvent = fReader->GetEventRec();
1559       
1560       ProcessRecAndSimNonId(rawtrackEvent,rawpartEvent);
1561   }//end of loop over events (1)
1562 }
1563 /*************************************************************************************/  
1564  
1565 void AliHBTAnalysis::ProcessTracksNonIdentAnal()
1566 {
1567 //Process Tracks only with non identical mode
1568   AliAOD * rawtrackEvent;
1569   fReader->Rewind();
1570   
1571   Info("ProcessTracksNonIdentAnal","**************************************");
1572   Info("ProcessTracksNonIdentAnal","*****  NON IDENT MODE ****************");
1573   Info("ProcessTracksNonIdentAnal","**************************************");
1574   
1575   for (Int_t i = 0;;i++)//infinite loop
1576     {
1577       if (fReader->Next()) break; //end when no more events available
1578       rawtrackEvent = fReader->GetEventRec();
1579       ProcessRecNonId(rawtrackEvent,0x0);
1580   }//end of loop over events (1)
1581 }
1582 /*************************************************************************************/  
1583
1584 void AliHBTAnalysis::ProcessParticlesNonIdentAnal()
1585 {
1586 //process paricles only with non identical mode
1587   AliAOD * rawpartEvent = 0x0;
1588   fReader->Rewind();
1589   
1590   Info("ProcessParticlesNonIdentAnal","**************************************");
1591   Info("ProcessParticlesNonIdentAnal","*****  NON IDENT MODE ****************");
1592   Info("ProcessParticlesNonIdentAnal","**************************************");
1593
1594   for (Int_t i = 0;;i++)//infinite loop
1595     {
1596       if (fReader->Next()) break; //end when no more events available
1597       
1598       rawpartEvent  = fReader->GetEventSim();
1599       ProcessSimNonId(0x0,rawpartEvent);
1600   }//end of loop over events (1)
1601 }
1602
1603 /*************************************************************************************/  
1604 void AliHBTAnalysis::FilterOut(AliAOD* outpart1, AliAOD* outpart2, AliAOD* inpart,
1605                      AliAOD* outtrack1, AliAOD* outtrack2, AliAOD* intrack) const
1606 {
1607  //Puts particles accepted as a first particle by global cut in out1
1608  //and as a second particle in out2
1609
1610   AliVAODParticle* part, *track;
1611
1612   outpart1->Reset();
1613   outpart2->Reset();
1614   outtrack1->Reset();
1615   outtrack2->Reset();
1616   
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      
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
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       {
1645         outpart2->AddParticle(new AliAODParticle(*part));
1646         outtrack2->AddParticle(new AliAODParticle(*track));
1647         continue;
1648       }
1649    }
1650 }
1651 /*************************************************************************************/  
1652
1653 void AliHBTAnalysis::FilterOut(AliAOD* out1, AliAOD* out2, AliAOD* in) const
1654 {
1655  //Puts particles accepted as a first particle by global cut in out1
1656  //and as a second particle in out2
1657   AliVAODParticle* part;
1658   
1659   out1->Reset();
1660   out2->Reset();
1661   
1662   AliAODParticleCut *cut1 = fPairCut->GetFirstPartCut();
1663   AliAODParticleCut *cut2 = fPairCut->GetSecondPartCut();
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       }
1682
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 /*************************************************************************************/ 
1697
1698 Bool_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  
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;
1714  
1715  Int_t id1 = fPairCut->GetFirstPartCut()->GetPID();
1716  Int_t id2 = fPairCut->GetSecondPartCut()->GetPID();
1717
1718  if ( (id1==0) || (id2==0) || (id1==id2) ) 
1719    return kFALSE;
1720  
1721  return kTRUE;
1722 }
1723 /*************************************************************************************/ 
1724
1725 void 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
1738 void AliHBTAnalysis::PressAnyKey()
1739
1740  //small utility function that helps to make comfortable macros
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
1752 /*************************************************************************************/ 
1753