Particle, Particle Cut moved to ANALYSIS
[u/mrichter/AliRoot.git] / HBTAN / AliHBTPairCut.cxx
1 /* $Id$ */
2 //____________________________________
3 /////////////////////////////////////////////////////////////////////////
4 //
5 // Class AliHBTPairCut:
6 //
7 // implements cut on the pair of particles
8 // more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
9 // Author: Piotr.Skowronski@cern.ch
10 //-------------------------------------------------------------------
11
12 #include "AliHBTPairCut.h"
13 #include "AliHBTPair.h"
14 #include "AliHBTParticleCut.h"
15 #include "AliHBTTrackPoints.h"
16 #include "AliHBTClusterMap.h"
17
18 ClassImp(AliHBTPairCut)
19 const Int_t AliHBTPairCut::fgkMaxCuts = 50;
20 /**********************************************************/
21
22 AliHBTPairCut::AliHBTPairCut():
23   fNCuts(0)
24 {
25   //constructor
26   fFirstPartCut = new AliHBTEmptyParticleCut(); //empty cuts
27   fSecondPartCut= new AliHBTEmptyParticleCut(); //empty cuts
28     
29   fCuts = new AliHbtBasePairCut*[fgkMaxCuts];
30   for (Int_t i = 0;i<fNCuts;i++)
31    {
32      fCuts[i] = 0x0;
33    }
34 }
35 /**********************************************************/
36
37 AliHBTPairCut::AliHBTPairCut(const AliHBTPairCut& in):
38  TNamed(in)
39 {
40   //copy constructor
41   fCuts = new AliHbtBasePairCut*[fgkMaxCuts];
42   fNCuts = in.fNCuts;
43
44   fFirstPartCut = (AliHBTParticleCut*)in.fFirstPartCut->Clone();
45   fSecondPartCut = (AliHBTParticleCut*)in.fSecondPartCut->Clone();
46  
47   for (Int_t i = 0;i<fNCuts;i++)
48     {
49       fCuts[i] = (AliHbtBasePairCut*)in.fCuts[i]->Clone();//create new object (clone) and rember pointer to it
50     }
51 }
52 /**********************************************************/
53
54 AliHBTPairCut&  AliHBTPairCut::operator=(const AliHBTPairCut& in)
55 {
56   //assignment operator
57   fCuts = new AliHbtBasePairCut*[fgkMaxCuts];
58   fNCuts = in.fNCuts;
59
60   fFirstPartCut = (AliHBTParticleCut*)in.fFirstPartCut->Clone();
61   fSecondPartCut = (AliHBTParticleCut*)in.fSecondPartCut->Clone();
62  
63   for (Int_t i = 0;i<fNCuts;i++)
64     {
65       fCuts[i] = (AliHbtBasePairCut*)in.fCuts[i]->Clone();//create new object (clone) and rember pointer to it
66     }
67   return * this;
68 }
69 /**********************************************************/
70
71 AliHBTPairCut::~AliHBTPairCut()
72 {
73   //destructor
74   if (fFirstPartCut != fSecondPartCut)
75     {
76       delete fSecondPartCut;
77     }
78   delete fFirstPartCut;
79   for (Int_t i = 0;i<fNCuts;i++)
80     {
81       delete fCuts[i];
82     }
83   delete []fCuts;
84
85 /**********************************************************/
86
87 /**********************************************************/
88
89 void AliHBTPairCut::AddBasePairCut(AliHbtBasePairCut* basecut)
90 {
91   //adds the base pair cut (cut on one value)
92   
93   if (!basecut) return;
94   if( fNCuts == (fgkMaxCuts-1) )
95     {
96       Warning("AddBasePairCut","Not enough place for another cut");
97       return;
98     }
99   fCuts[fNCuts++]=basecut;
100 }
101 /**********************************************************/
102
103 Bool_t AliHBTPairCut::Pass(AliHBTPair* pair) const
104 {
105   //methods which checks if given pair meets all criteria of the cut
106   //if it meets returns FALSE
107   //if NOT   returns    TRUE
108   if(!pair) 
109     {
110       Warning("Pass","No Pasaran! We never accept NULL pointers");
111       return kTRUE;
112     }
113   
114   //check particle's cuts
115   if( ( fFirstPartCut->Pass( pair->Particle1()) ) || 
116       ( fSecondPartCut->Pass(pair->Particle2()) )   )
117     {  
118       return kTRUE;
119     }
120   return PassPairProp(pair);
121 }
122 /**********************************************************/
123
124 Bool_t AliHBTPairCut::PassPairProp(AliHBTPair* pair) const
125 {
126   //methods which checks if given pair meets all criteria of the cut
127   //if it meets returns FALSE
128   //if NOT   returns    TRUE
129   //examine all base pair cuts
130   for (Int_t i = 0;i<fNCuts;i++)
131     {
132       if ( (fCuts[i]->Pass(pair)) ) return kTRUE; //if one of the cuts reject, then reject
133     }
134   return kFALSE;
135 }
136 /**********************************************************/
137
138 void AliHBTPairCut::Print()
139 {
140  //Prints the cut
141   for (Int_t i = 0;i<fNCuts;i++)
142     {
143       fCuts[i]->Dump();
144     }
145 }
146 /**********************************************************/
147
148 void AliHBTPairCut::SetFirstPartCut(AliHBTParticleCut* cut)
149 {
150   // set cut for the first particle
151   if(!cut) 
152     {
153       Error("SetFirstPartCut","argument is NULL");
154       return;
155     }
156   delete fFirstPartCut;
157   fFirstPartCut = (AliHBTParticleCut*)cut->Clone();
158   
159 }
160 /**********************************************************/
161
162 void AliHBTPairCut::SetSecondPartCut(AliHBTParticleCut* cut)
163 {
164   // set cut for the second particle
165   if(!cut) 
166     {
167       Error("SetSecondPartCut","argument is NULL");
168       return;
169     }
170   delete fSecondPartCut;
171   fSecondPartCut = (AliHBTParticleCut*)cut->Clone();
172 }
173 /**********************************************************/
174
175 void AliHBTPairCut::SetPartCut(AliHBTParticleCut* cut)
176 {
177   //sets the the same cut on both particles
178   if(!cut) 
179     {
180       Error("SetFirstPartCut","argument is NULL");
181       return;
182     }
183   if (fFirstPartCut == fSecondPartCut) fSecondPartCut = 0x0;
184   
185   delete fFirstPartCut;
186   fFirstPartCut = (AliHBTParticleCut*)cut->Clone();
187   
188   delete fSecondPartCut; //even if null should not be harmful
189   fSecondPartCut = fFirstPartCut;
190 }
191 /**********************************************************/
192
193 void AliHBTPairCut::SetQInvRange(Double_t min, Double_t max)
194 {
195   // set range of accepted invariant masses
196   AliHBTQInvCut* cut= (AliHBTQInvCut*)FindCut(kHbtPairCutPropQInv);
197   if(cut) cut->SetRange(min,max);
198   else fCuts[fNCuts++] = new AliHBTQInvCut(min,max);
199 }
200 /**********************************************************/
201 void AliHBTPairCut::SetQOutCMSLRange(Double_t min, Double_t max)
202 {
203   // set range of accepted QOut in CMS
204   AliHBTQOutLCMSCut* cut= (AliHBTQOutLCMSCut*)FindCut(kHbtPairCutPropQOutLCMS);
205   if(cut) cut->SetRange(min,max);
206   else fCuts[fNCuts++] = new AliHBTQOutLCMSCut(min,max);
207 }
208
209 /**********************************************************/
210 void AliHBTPairCut::SetQSideCMSLRange(Double_t min, Double_t max)
211 {
212   // set range of accepted QSide in CMS
213   AliHBTQSideLCMSCut* cut= (AliHBTQSideLCMSCut*)FindCut(kHbtPairCutPropQSideLCMS);
214   if(cut) cut->SetRange(min,max);
215   else fCuts[fNCuts++] = new AliHBTQSideLCMSCut(min,max);
216 }
217
218 /**********************************************************/
219 void AliHBTPairCut::SetQLongCMSLRange(Double_t min, Double_t max)
220 {
221   // set range of accepted QLong in CMS
222   AliHBTQLongLCMSCut* cut= (AliHBTQLongLCMSCut*)FindCut(kHbtPairCutPropQLongLCMS);
223   if(cut) cut->SetRange(min,max);
224   else fCuts[fNCuts++] = new AliHBTQLongLCMSCut(min,max);
225 }
226
227 /**********************************************************/
228
229 void AliHBTPairCut::SetKtRange(Double_t min, Double_t max)
230 {
231   // set range of accepted Kt (?)
232   AliHBTKtCut* cut= (AliHBTKtCut*)FindCut(kHbtPairCutPropKt);
233   if(cut) cut->SetRange(min,max);
234   else fCuts[fNCuts++] = new AliHBTKtCut(min,max);
235 }
236 /**********************************************************/
237
238 void AliHBTPairCut::SetKStarRange(Double_t min, Double_t max)
239 {
240   // set range of accepted KStar (?)
241   AliHBTKStarCut* cut= (AliHBTKStarCut*)FindCut(kHbtPairCutPropKStar);
242   if(cut) cut->SetRange(min,max);
243   else fCuts[fNCuts++] = new AliHBTKStarCut(min,max);
244 }
245 /**********************************************************/
246
247 void AliHBTPairCut::SetAvSeparationRange(Double_t min, Double_t max)
248 {
249   //sets avarage separation cut ->Anti-Merging cut
250   AliHbtBasePairCut* cut= FindCut(kHbtPairCutPropAvSepar);
251   if(cut) cut->SetRange(min,max);
252   else fCuts[fNCuts++] = new AliHBTAvSeparationCut(min,max);
253 }
254 /**********************************************************/
255
256 void AliHBTPairCut::SetITSSeparation(Int_t layer, Double_t drphi, Double_t dz)
257 {
258   //Anti-Merging Cut for first pixel layer
259   AliHBTITSSeparationCut* cut= dynamic_cast<AliHBTITSSeparationCut*>(FindCut(kHbtPairCutPropPixelSepar));
260   if(cut) 
261    {
262      if (layer == cut->GetLayer())
263       {
264         cut->SetRange(drphi,dz);//In this cut fMin is drphi, and fMax dz
265         return;
266       }
267    }
268   fCuts[fNCuts++] = new AliHBTITSSeparationCut(layer,drphi,dz);
269 //  Info("SetITSSeparation","Added %d at address %#x",fNCuts-1,fCuts[fNCuts-1]);
270 }
271 /**********************************************************/
272
273 void AliHBTPairCut::SetClusterOverlapRange(Double_t min,Double_t max)
274 {
275   //sets cluster overlap factor cut ->Anti-Splitting cut
276   //cluster overlap factor ranges between 
277   // -0.5 (in all padrows both tracks have cluters) 
278   // and 1 (in all padrows one track has cluter and second has not)
279   // When Overlap Factor is 1 this pair of tracks in highly probable to be
280   // splitted track: one particle that is recontructed twise
281   // STAR uses range from -0.5 to 0.6 
282   
283   AliHbtBasePairCut* cut= FindCut(kHbtPairCutPropClOverlap);
284   if(cut) cut->SetRange(min,max);
285   else fCuts[fNCuts++] = new AliHBTCluterOverlapCut(min,max);
286 }
287 /**********************************************************/
288
289 AliHbtBasePairCut* AliHBTPairCut::FindCut(AliHBTPairCutProperty property)
290 {
291   // Find the cut corresponding to "property"
292   for (Int_t i = 0;i<fNCuts;i++)
293     {
294       if (fCuts[i]->GetProperty() == property) 
295         return fCuts[i]; //we found the cut we were searching for
296     }
297   
298   return 0x0; //we did not found this cut
299   
300 }
301 /**********************************************************/
302
303 void AliHBTPairCut::Streamer(TBuffer &b)
304 {
305   // Stream all objects in the array to or from the I/O buffer.
306   
307   UInt_t R__s, R__c;
308   if (b.IsReading()) 
309     {
310       Version_t v = b.ReadVersion(&R__s, &R__c);
311       if (v > -1)
312        {
313           delete fFirstPartCut;
314           delete fSecondPartCut;
315           fFirstPartCut = 0x0;
316           fSecondPartCut = 0x0;
317           TObject::Streamer(b);
318           b >> fFirstPartCut;
319           b >> fSecondPartCut;
320           b >> fNCuts;
321           for (Int_t i = 0;i<fNCuts;i++)
322            {
323              b >> fCuts[i];
324            }
325         }
326       b.CheckByteCount(R__s, R__c,AliHBTPairCut::IsA());
327     } 
328   else 
329     {
330       R__c = b.WriteVersion(AliHBTPairCut::IsA(), kTRUE);
331       TObject::Streamer(b);
332       
333 //      printf("Streamer Cut 1 %#x Cut 2 %#x\n",fFirstPartCut,fSecondPartCut);
334 //      this->Dump();
335 //      fFirstPartCut->Dump();
336       
337       b << fFirstPartCut;
338       b << fSecondPartCut;
339       b << fNCuts;
340       for (Int_t i = 0;i<fNCuts;i++)
341         {
342           b << fCuts[i];
343         }
344       b.SetByteCount(R__c, kTRUE);
345     }
346 }
347 /******************************************************************/
348
349 ClassImp(AliHBTEmptyPairCut)
350   
351 void AliHBTEmptyPairCut::Streamer(TBuffer &b)
352 {
353 //streamer for empty pair cut
354   AliHBTPairCut::Streamer(b);
355 }
356 /******************************************************************/
357
358 ClassImp(AliHbtBasePairCut)
359 ClassImp(AliHBTQInvCut)
360 ClassImp(AliHBTKtCut)
361 ClassImp(AliHBTQSideLCMSCut)
362 ClassImp(AliHBTQOutLCMSCut)
363 ClassImp(AliHBTQLongLCMSCut)
364
365 /******************************************************************/
366 ClassImp(AliHBTAvSeparationCut)
367     
368 Double_t AliHBTAvSeparationCut::GetValue(AliHBTPair* pair) const 
369 {
370   //chacks if avarage distance of two tracks is in given range
371   AliHBTTrackPoints* tpts1 = pair->Particle1()->GetTrackPoints();
372   if ( tpts1 == 0x0)
373    {//it could be simulated pair
374 //     Warning("GetValue","Track 1 does not have Track Points. Pair NOT Passed.");
375      return -1.0;
376    }
377
378   AliHBTTrackPoints* tpts2 = pair->Particle2()->GetTrackPoints();
379   if ( tpts2 == 0x0)
380    {
381 //     Warning("GetValue","Track 2 does not have Track Points. Pair NOT Passed.");
382      return -1.0;
383    }
384    
385   return tpts1->AvarageDistance(*tpts2);
386 }
387 /******************************************************************/
388 ClassImp(AliHBTSeparationCut)
389     
390 Double_t AliHBTSeparationCut::GetValue(AliHBTPair* pair) const 
391 {
392   //chacks if avarage distance of two tracks is in given range
393   AliHBTTrackPoints* tpts1 = pair->Particle1()->GetTrackPoints();
394   if ( tpts1 == 0x0)
395    {//it could be simulated pair
396 //     Warning("GetValue","Track 1 does not have Track Points. Pair NOT Passed.");
397      return -1.0;
398    }
399
400   AliHBTTrackPoints* tpts2 = pair->Particle2()->GetTrackPoints();
401   if ( tpts2 == 0x0)
402    {
403 //     Warning("GetValue","Track 2 does not have Track Points. Pair NOT Passed.");
404      return -1.0;
405    }
406   Float_t x1=0,y1=0,z1=0; 
407   Float_t x2=0,y2=0,z2=0;
408   
409   tpts1->PositionAt(fPoint,x1,y1,z1);
410   tpts2->PositionAt(fPoint,x2,y2,z2);
411   Double_t dx1 = x1 - x2;
412   Double_t dy1 = y1 - y2;
413   Double_t dz1 = z1 - z2;
414   Double_t d = TMath::Sqrt(dx1*dx1 + dy1*dy1 + dz1*dz1);
415   return d;
416 }
417 /******************************************************************/
418
419 ClassImp(AliHBTITSSeparationCut)
420
421 Bool_t AliHBTITSSeparationCut::Pass(AliHBTPair* pair) const
422 {
423  //Checks if two tracks do not cross first pixels too close to each other
424  //If two tracks use the same cluster in pixels they are given
425  //the same position what skews theta angles (both are the same)
426  //These guys create artificial correlation in non-id analyses
427  //which is positive for identical polar angles (Qlong=0) 
428  //and negative for a little bit different theta angle (Qlong=epsilon)
429  //Such tracks "attracks" each other.
430  
431   AliHBTTrackPoints* tpts1 = pair->Particle1()->GetITSTrackPoints();
432   if ( tpts1 == 0x0)
433    {//it could be simulated pair
434      Warning("Pass","Track 1 does not have ITS Track Points. Pair NOT Passed.");
435      return kTRUE;//reject 
436    }
437
438   AliHBTTrackPoints* tpts2 = pair->Particle2()->GetITSTrackPoints();
439   if ( tpts2 == 0x0)
440    {
441      Warning("Pass","Track 2 does not have ITS Track Points. Pair NOT Passed.");
442      return kTRUE;//reject 
443    }
444   Float_t  x1=0.0,y1=0.0,z1=0.0,x2=0.0,y2=0.0,z2=0.0;
445   tpts1->PositionAt(fLayer,x1,y1,z1);
446   tpts2->PositionAt(fLayer,x2,y2,z2);
447   
448 //  Info("Pass","rphi %f z %f",fMin,fMax);
449 //  Info("Pass","P1: %f %f %f", x1,y1,z1);
450 //  Info("Pass","P2: %f %f %f", x2,y2,z2);
451   
452   Double_t dz = TMath::Abs(z1-z2);
453   
454   //fMax encodes treshold valaue of distance in Z
455   if (dz > fMax) return kFALSE;//pair accepted
456   
457   Double_t drphi = TMath::Hypot(x1-x2,y1-y2);
458   
459   //fMin encodes treshold valaue of distance in r-phi
460   if (drphi > fMin) return kFALSE;
461   
462   return kTRUE;//they are too close, rejected
463 }
464 /******************************************************************/
465
466 ClassImp(AliHBTCluterOverlapCut)
467
468 Double_t  AliHBTCluterOverlapCut::GetValue(AliHBTPair* pair) const
469 {
470   //Returns Cluter Overlap Factor
471   //It ranges between -0.5 (in all padrows both tracks have cluters) 
472   // and 1 (in all padrows one track has cluter and second has not)
473   // When Overlap Factor is 1 this pair of tracks in highly probable to be
474   // splitted track: one particle that is recontructed twise
475
476   AliHBTClusterMap* cm1 = pair->Particle1()->GetClusterMap();
477   if ( cm1 == 0x0)
478    {
479      Warning("GetValue","Track 1 does not have Cluster Map. Returning -0.5.");
480      return -.5;
481    }
482
483   AliHBTClusterMap* cm2 = pair->Particle2()->GetClusterMap();
484   if ( cm2 == 0x0)
485    {
486      Warning("GetValue","Track 2 does not have Cluster Map. Returning -0.5.");
487      return -.5;
488    }
489   return cm1->GetOverlapFactor(*cm2);
490 }
491 /******************************************************************/
492 ClassImp(AliHBTOutSideSameSignCut)
493
494 Bool_t AliHBTOutSideSameSignCut::Pass(AliHBTPair *p) const
495 {
496   //returns kTRUE if pair DO NOT meet cut criteria
497   
498   if ( p->GetQOutLCMS()*p->GetQSideLCMS() > 0 ) 
499    {
500      return kFALSE;//accpeted
501    }
502
503   return kTRUE ;//rejected
504 }
505 /******************************************************************/
506 ClassImp(AliHBTOutSideDiffSignCut)
507
508 Bool_t AliHBTOutSideDiffSignCut::Pass(AliHBTPair *p) const
509 {
510   //returns kTRUE if pair DO NOT meet cut criteria
511   
512   if ( p->GetQOutLCMS()*p->GetQSideLCMS() > 0 ) 
513    {
514      return kTRUE;//rejected
515    }
516   
517   return kFALSE;//accepted
518 }
519 /******************************************************************/
520 ClassImp( AliHBTLogicalOperPairCut )
521
522 AliHBTLogicalOperPairCut::AliHBTLogicalOperPairCut():
523  AliHbtBasePairCut(-10e10,10e10,kHbtPairCutPropNone),
524  fFirst(new AliHBTDummyBasePairCut),
525  fSecond(new AliHBTDummyBasePairCut)
526 {
527  //ctor
528 }
529 /******************************************************************/
530
531 AliHBTLogicalOperPairCut::AliHBTLogicalOperPairCut(AliHbtBasePairCut* first, AliHbtBasePairCut* second):
532  AliHbtBasePairCut(-10e10,10e10,kHbtPairCutPropNone),
533  fFirst((first)?(AliHbtBasePairCut*)first->Clone():0x0),
534  fSecond((second)?(AliHbtBasePairCut*)second->Clone():0x0)
535 {
536   //ctor
537   //note that base cuts are copied, not just pointers assigned
538   if ( (fFirst && fSecond) == kFALSE) 
539    {
540      Fatal("AliHBTLogicalOperPairCut","One of parameters is NULL!");
541    }
542 }
543 /******************************************************************/
544
545 AliHBTLogicalOperPairCut::~AliHBTLogicalOperPairCut()
546 {
547   //destructor
548   delete fFirst;
549   delete fSecond;
550 }
551 /******************************************************************/
552
553 Bool_t AliHBTLogicalOperPairCut::AliHBTDummyBasePairCut::Pass(AliHBTPair* /*pair*/)  const
554 {
555   //checks if particles passes properties defined by this cut
556   Warning("Pass","You are using dummy base cut! Probobly some logical cut is not set up properly");
557   return kFALSE;//accept
558 }
559 /******************************************************************/
560
561 void AliHBTLogicalOperPairCut::Streamer(TBuffer &b)
562 {
563   // Stream all objects in the array to or from the I/O buffer.
564   UInt_t R__s, R__c;
565   if (b.IsReading()) 
566    {
567      delete fFirst;
568      delete fSecond;
569      fFirst  = 0x0;
570      fSecond = 0x0;
571
572      b.ReadVersion(&R__s, &R__c);
573      TObject::Streamer(b);
574      b >> fFirst;
575      b >> fSecond;
576      b.CheckByteCount(R__s, R__c,AliHBTLogicalOperPairCut::IsA());
577    } 
578   else 
579    {
580      R__c = b.WriteVersion(AliHBTLogicalOperPairCut::IsA(), kTRUE);
581      TObject::Streamer(b);
582      b << fFirst;
583      b << fSecond;
584      b.SetByteCount(R__c, kTRUE);
585   }
586 }
587
588 /******************************************************************/
589 ClassImp(AliHBTOrPairCut)
590
591 Bool_t AliHBTOrPairCut::Pass(AliHBTPair * p) const
592 {
593   //returns true when rejected 
594   //AND operation is a little bit misleading but is correct
595   //User wants to build logical cuts with natural (positive) logic
596   //while HBTAN use inernally reverse (returns true when rejected)
597   if (fFirst->Pass(p) && fSecond->Pass(p) ) return kTRUE;//rejected (both rejected, returned kTRUE)
598   return kFALSE;//accepted, at least one accepted (returned kFALSE)
599 }
600 /******************************************************************/
601
602 ClassImp(AliHBTAndPairCut)
603
604 Bool_t AliHBTAndPairCut::Pass(AliHBTPair * p)  const
605 {
606   //returns true when rejected 
607   //OR operation is a little bit misleading but is correct
608   //User wants to build logical cuts with natural (positive) logic
609   //while HBTAN use inernally reverse (returns true when rejected)
610   if (fFirst->Pass(p) || fSecond->Pass(p)) return kTRUE;//rejected (any of two rejected(returned kTRUE) )
611   return kFALSE;//accepted (both accepted (returned kFALSE))
612 }
613 /******************************************************************/