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