PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoESDTrackCut.cxx

   1 /***************************************************************************
   2  *
   3  * $Id$
   4  *
   5  *
   6  ***************************************************************************
   7  *
   8  *
   9  *
  10  *
  11  ***************************************************************************
  12  *
  13  * $Log$
  14  * Revision 1.3  2007/05/22 09:01:42  akisiel
  15  * Add the possibiloity to save cut settings in the ROOT file
  16  *
  17  * Revision 1.2  2007/05/21 10:38:25  akisiel
  18  * More coding rule conformance
  19  *
  20  * Revision 1.1  2007/05/16 10:25:06  akisiel
  21  * Making the directory structure of AliFemtoUser flat. All files go into one common directory
  22  *
  23  * Revision 1.4  2007/05/03 09:46:10  akisiel
  24  * Fixing Effective C++ warnings
  25  *
  26  * Revision 1.3  2007/04/27 07:25:59  akisiel
  27  * Make revisions needed for compilation from the main AliRoot tree
  28  *
  29  * Revision 1.1.1.1  2007/04/25 15:38:41  panos
  30  * Importing the HBT code dir
  31  *
  32  * Revision 1.4  2007-04-03 16:00:08  mchojnacki
  33  * Changes to iprove memory managing
  34  *
  35  * Revision 1.3  2007/03/13 15:30:03  mchojnacki
  36  * adding reader for simulated data
  37  *
  38  * Revision 1.2  2007/03/08 14:58:03  mchojnacki
  39  * adding some alice stuff
  40  *
  41  * Revision 1.1.1.1  2007/03/07 10:14:49  mchojnacki
  42  * First version on CVS
  43  *
  44  **************************************************************************/
  45
  46 #include "AliFemtoESDTrackCut.h"
  47 #include <cstdio>
  48
  49 #ifdef __ROOT__
  50 ClassImp(AliFemtoESDTrackCut)
  51 #endif
  52
  53
  54 // electron
  55 // 0.13 - 1.8
  56 // 0       7.594129e-02    8.256141e-03
  57 // 1       -5.535827e-01   8.170825e-02
  58 // 2       1.728591e+00    3.104210e-01
  59 // 3       -2.827893e+00   5.827802e-01
  60 // 4       2.503553e+00    5.736207e-01
  61 // 5       -1.125965e+00   2.821170e-01
  62 // 6       2.009036e-01    5.438876e-02
  63
  64 // pion
  65 // 0.13 - 2.0
  66 // 0       1.063457e+00    8.872043e-03
  67 // 1       -4.222208e-01   2.534402e-02
  68 // 2       1.042004e-01    1.503945e-02
  69
  70 // kaon
  71 // 0.18 - 2.0
  72 // 0       -7.289406e-02   1.686074e-03
  73 // 1       4.415666e-01    1.143939e-02
  74 // 2       -2.996790e-01   1.840964e-02
  75 // 3       6.704652e-02    7.783990e-03
  76
  77 // proton
  78 // 0.26 - 2.0
  79 // 0       -3.730200e-02   2.347311e-03
  80 // 1       1.163684e-01    1.319316e-02
  81 // 2       8.354116e-02    1.997948e-02
  82 // 3       -4.608098e-02   8.336400e-03
  83
  84
  85 AliFemtoESDTrackCut::AliFemtoESDTrackCut() :
  86     fCharge(0),
  87     fLabel(0),
  88     fStatus(0),
  89     fminTPCclsF(0),
  90     fminITScls(0),
  91     fMaxITSchiNdof(1000.0),
  92     fMaxTPCchiNdof(1000.0),
  93     fMaxSigmaToVertex(1000.0),
  94     fNTracksPassed(0),
  95     fNTracksFailed(0),
  96     fRemoveKinks(kFALSE),
  97     fMostProbable(0)
  98 {
  99   // Default constructor
 100   fNTracksPassed = fNTracksFailed = 0;
 101   fCharge = 0;  // takes both charges 0
 102   fPt[0]=0.0;              fPt[1] = 100.0;//100
 103   fRapidity[0]=-2;       fRapidity[1]=2;//-2 2
 104   fPidProbElectron[0]=-1;fPidProbElectron[1]=2;
 105   fPidProbPion[0]=-1;    fPidProbPion[1]=2;
 106   fPidProbKaon[0]=-1;fPidProbKaon[1]=2;
 107   fPidProbProton[0]=-1;fPidProbProton[1]=2;
 108   fPidProbMuon[0]=-1;fPidProbMuon[1]=2;
 109   fLabel=false;
 110   fStatus=0;
 111   fminTPCclsF=0;
 112   fminITScls=0;
 113 }
 114 //------------------------------
 115 AliFemtoESDTrackCut::~AliFemtoESDTrackCut(){
 116   /* noop */
 117 }
 118 //------------------------------
 119 bool AliFemtoESDTrackCut::Pass(const AliFemtoTrack* track)
 120 {
 121   // test the particle and return
 122   // true if it meets all the criteria
 123   // false if it doesn't meet at least one of the criteria
 124   float tMost[5];
 125
 126   //cout<<"AliFemtoESD  cut"<<endl;
 127   //cout<<fPidProbPion[0]<<" < pi ="<<track->PidProbPion()<<" <"<<fPidProbPion[1]<<endl;
 128   if (fStatus!=0)
 129     {
 130       //cout<<" status "<<track->Label()<<" "<<track->Flags()<<" "<<track->TPCnclsF()<<" "<<track->ITSncls()<<endl;
 131       if ((track->Flags()&fStatus)!=fStatus)
 132         {
 133           //      cout<<track->Flags()<<" "<<fStatus<<" no go through status"<<endl;
 134           return false;
 135         }
 136
 137     }
 138   if (fRemoveKinks) {
 139     if ((track->KinkIndex(0)) || (track->KinkIndex(1)) || (track->KinkIndex(2)))
 140       return false;
 141   }
 142   if (fminTPCclsF>track->TPCnclsF())
 143     {
 144       //cout<<" No go because TPC Number of ClsF"<<fminTPCclsF<< " "<<track->TPCnclsF()<<endl;
 145       return false;
 146     }
 147   if (fminTPCncls>track->TPCncls())
 148     {
 149       //cout<<" No go because TPC Number of ClsF"<<fminTPCclsF<< " "<<track->TPCnclsF()<<endl;
 150       return false;
 151     }
 152   if (fminITScls>track->ITSncls())
 153     {
 154       //cout<<" No go because ITS Number of Cls"<<fminITScls<< " "<<track->ITSncls()<<endl;
 155       return false;
 156     }
 157
 158   if (fMaxSigmaToVertex < track->SigmaToVertex()) {
 159     return false;
 160   }
 161
 162   if (track->ITSncls() > 0)
 163     if ((track->ITSchi2()/track->ITSncls()) > fMaxITSchiNdof) {
 164       return false;
 165     }
 166
 167   if (track->TPCncls() > 0)
 168     if ((track->TPCchi2()/track->TPCncls()) > fMaxTPCchiNdof) {
 169       return false;
 170     }
 171
 172   if (fLabel)
 173     {
 174       //cout<<"labels"<<endl;
 175       if(track->Label()<0)
 176         {
 177           fNTracksFailed++;
 178           //   cout<<"No Go Through the cut"<<endl;
 179           //  cout<<fLabel<<" Label="<<track->Label()<<endl;
 180           return false;
 181         }
 182     }
 183   if (fCharge!=0)
 184     {
 185       //cout<<"AliFemtoESD  cut ch "<<endl;
 186       //cout<<fCharge<<" Charge="<<track->Charge()<<endl;
 187       if (track->Charge()!= fCharge)
 188         {
 189           fNTracksFailed++;
 190           //  cout<<"No Go Through the cut"<<endl;
 191           // cout<<fCharge<<" Charge="<<track->Charge()<<endl;
 192           return false;
 193         }
 194     }
 195   float tEnergy = ::sqrt(track->P().mag2()+fMass*fMass);
 196   float tRapidity = 0.5*::log((tEnergy+track->P().z())/(tEnergy-track->P().z()));
 197   float tPt = ::sqrt((track->P().x())*(track->P().x())+(track->P().y())*(track->P().y()));
 198   if ((tRapidity<fRapidity[0])||(tRapidity>fRapidity[1]))
 199     {
 200       fNTracksFailed++;
 201       //cout<<"No Go Through the cut"<<endl;
 202       //cout<<fRapidity[0]<<" < Rapidity ="<<tRapidity<<" <"<<fRapidity[1]<<endl;
 203       return false;
 204     }
 205   if ((tPt<fPt[0])||(tPt>fPt[1]))
 206     {
 207       fNTracksFailed++;
 208       //cout<<"No Go Through the cut"<<endl;
 209       //cout<<fPt[0]<<" < Pt ="<<Pt<<" <"<<fPt[1]<<endl;
 210       return false;
 211     }
 212 //   cout << "Track has pids: "
 213 //        << track->PidProbElectron() << " "
 214 //        << track->PidProbMuon() << " "
 215 //        << track->PidProbPion() << " "
 216 //        << track->PidProbKaon() << " "
 217 //        << track->PidProbProton() << " "
 218 //        << track->PidProbElectron()+track->PidProbMuon()+track->PidProbPion()+track->PidProbKaon()+track->PidProbProton() << endl;
 219
 220
 221   if ((track->PidProbElectron()<fPidProbElectron[0])||(track->PidProbElectron()>fPidProbElectron[1]))
 222     {
 223       fNTracksFailed++;
 224       //cout<<"No Go Through the cut"<<endl;
 225       //cout<<fPidProbElectron[0]<<" < e ="<<track->PidProbElectron()<<" <"<<fPidProbElectron[1]<<endl;
 226       return false;
 227     }
 228   if ((track->PidProbPion()<fPidProbPion[0])||(track->PidProbPion()>fPidProbPion[1]))
 229     {
 230       fNTracksFailed++;
 231       //cout<<"No Go Through the cut"<<endl;
 232       //cout<<fPidProbPion[0]<<" < pi ="<<track->PidProbPion()<<" <"<<fPidProbPion[1]<<endl;
 233       return false;
 234     }
 235   if ((track->PidProbKaon()<fPidProbKaon[0])||(track->PidProbKaon()>fPidProbKaon[1]))
 236     {
 237       fNTracksFailed++;
 238       //cout<<"No Go Through the cut"<<endl;
 239       //cout<<fPidProbKaon[0]<<" < k ="<<track->PidProbKaon()<<" <"<<fPidProbKaon[1]<<endl;
 240       return false;
 241     }
 242   if ((track->PidProbProton()<fPidProbProton[0])||(track->PidProbProton()>fPidProbProton[1]))
 243     {
 244       fNTracksFailed++;
 245       //cout<<"No Go Through the cut"<<endl;
 246       //cout<<fPidProbProton[0]<<" < p  ="<<track->PidProbProton()<<" <"<<fPidProbProton[1]<<endl;
 247       return false;
 248     }
 249   if ((track->PidProbMuon()<fPidProbMuon[0])||(track->PidProbMuon()>fPidProbMuon[1]))
 250     {
 251       fNTracksFailed++;
 252       //cout<<"No Go Through the cut"<<endl;
 253       //cout<<fPidProbMuon[0]<<" <  mi="<<track->PidProbMuon()<<" <"<<fPidProbMuon[1]<<endl;
 254       return false;
 255     }
 256
 257   if (fMostProbable) {
 258     tMost[0] = track->PidProbElectron()*PidFractionElectron(track->P().mag());
 259     tMost[1] = 0.0;
 260     tMost[2] = track->PidProbPion()*PidFractionPion(track->P().mag());
 261     tMost[3] = track->PidProbKaon()*PidFractionKaon(track->P().mag());
 262     tMost[4] = track->PidProbProton()*PidFractionProton(track->P().mag());
 263     int imost=0;
 264     float ipidmax = 0.0;
 265     for (int ip=0; ip<5; ip++)
 266       if (tMost[ip] > ipidmax) { ipidmax = tMost[ip]; imost = ip; };
 267     if (imost != fMostProbable) return false;
 268   }
 269
 270   // cout<<"Go Through the cut"<<endl;
 271   // cout<<fLabel<<" Label="<<track->Label()<<endl;
 272   // cout<<fCharge<<" Charge="<<track->Charge()<<endl;
 273   // cout<<fPt[0]<<" < Pt ="<<Pt<<" <"<<fPt[1]<<endl;
 274   //cout<<fRapidity[0]<<" < Rapidity ="<<tRapidity<<" <"<<fRapidity[1]<<endl;
 275   //cout<<fPidProbElectron[0]<<" <  e="<<track->PidProbElectron()<<" <"<<fPidProbElectron[1]<<endl;
 276   //cout<<fPidProbPion[0]<<" <  pi="<<track->PidProbPion()<<" <"<<fPidProbPion[1]<<endl;
 277   //cout<<fPidProbKaon[0]<<" <  k="<<track->PidProbKaon()<<" <"<<fPidProbKaon[1]<<endl;
 278   //cout<<fPidProbProton[0]<<" <  p="<<track->PidProbProton()<<" <"<<fPidProbProton[1]<<endl;
 279   //cout<<fPidProbMuon[0]<<" <  mi="<<track->PidProbMuon()<<" <"<<fPidProbMuon[1]<<endl;
 280   fNTracksPassed++ ;
 281   return true;
 282
 283
 284 }
 285 //------------------------------
 286 AliFemtoString AliFemtoESDTrackCut::Report()
 287 {
 288   // Prepare report from the execution
 289   string tStemp;
 290   char tCtemp[100];
 291   sprintf(tCtemp,"Particle mass:\t%E\n",this->Mass());
 292   tStemp=tCtemp;
 293   sprintf(tCtemp,"Particle charge:\t%d\n",fCharge);
 294   tStemp+=tCtemp;
 295   sprintf(tCtemp,"Particle pT:\t%E - %E\n",fPt[0],fPt[1]);
 296   tStemp+=tCtemp;
 297   sprintf(tCtemp,"Particle rapidity:\t%E - %E\n",fRapidity[0],fRapidity[1]);
 298   tStemp+=tCtemp;
 299   sprintf(tCtemp,"Number of tracks which passed:\t%ld  Number which failed:\t%ld\n",fNTracksPassed,fNTracksFailed);
 300   tStemp += tCtemp;
 301   AliFemtoString returnThis = tStemp;
 302   return returnThis;
 303 }
 304 TList *AliFemtoESDTrackCut::ListSettings()
 305 {
 306   // return a list of settings in a writable form
 307   TList *tListSetttings = new TList();
 308   char buf[200];
 309   snprintf(buf, 200, "AliFemtoESDTrackCut.mass=%lf", this->Mass());
 310   tListSetttings->AddLast(new TObjString(buf));
 311
 312   snprintf(buf, 200, "AliFemtoESDTrackCut.charge=%i", fCharge);
 313   tListSetttings->AddLast(new TObjString(buf));
 314   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobpion.minimum=%lf", fPidProbPion[0]);
 315   tListSetttings->AddLast(new TObjString(buf));
 316   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobpion.maximum=%lf", fPidProbPion[1]);
 317   tListSetttings->AddLast(new TObjString(buf));
 318   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobkaon.minimum=%lf", fPidProbKaon[0]);
 319   tListSetttings->AddLast(new TObjString(buf));
 320   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobkaon.maximum=%lf", fPidProbKaon[1]);
 321   tListSetttings->AddLast(new TObjString(buf));
 322   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobproton.minimum=%lf", fPidProbProton[0]);
 323   tListSetttings->AddLast(new TObjString(buf));
 324   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobproton.maximum=%lf", fPidProbProton[1]);
 325   tListSetttings->AddLast(new TObjString(buf));
 326   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobelectron.minimum=%lf", fPidProbElectron[0]);
 327   tListSetttings->AddLast(new TObjString(buf));
 328   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobelectron.maximum=%lf", fPidProbElectron[1]);
 329   tListSetttings->AddLast(new TObjString(buf));
 330   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobMuon.minimum=%lf", fPidProbMuon[0]);
 331   tListSetttings->AddLast(new TObjString(buf));
 332   snprintf(buf, 200, "AliFemtoESDTrackCut.pidprobMuon.maximum=%lf", fPidProbMuon[1]);
 333   tListSetttings->AddLast(new TObjString(buf));
 334   snprintf(buf, 200, "AliFemtoESDTrackCut.minimumtpcclusters=%i", fminTPCclsF);
 335   tListSetttings->AddLast(new TObjString(buf));
 336   snprintf(buf, 200, "AliFemtoESDTrackCut.minimumitsclusters=%i", fminTPCclsF);
 337   tListSetttings->AddLast(new TObjString(buf));
 338   snprintf(buf, 200, "AliFemtoESDTrackCut.pt.minimum=%lf", fPt[0]);
 339   tListSetttings->AddLast(new TObjString(buf));
 340   snprintf(buf, 200, "AliFemtoESDTrackCut.pt.maximum=%lf", fPt[1]);
 341   tListSetttings->AddLast(new TObjString(buf));
 342   snprintf(buf, 200, "AliFemtoESDTrackCut.rapidity.minimum=%lf", fRapidity[0]);
 343   tListSetttings->AddLast(new TObjString(buf));
 344   snprintf(buf, 200, "AliFemtoESDTrackCut.rapidity.maximum=%lf", fRapidity[1]);
 345   tListSetttings->AddLast(new TObjString(buf));
 346   snprintf(buf, 200, "AliFemtoESDTrackCut.removekinks=%i", fRemoveKinks);
 347   tListSetttings->AddLast(new TObjString(buf));
 348   snprintf(buf, 200, "AliFemtoESDTrackCut.maxitschindof=%lf", fMaxITSchiNdof);
 349   tListSetttings->AddLast(new TObjString(buf));
 350   snprintf(buf, 200, "AliFemtoESDTrackCut.maxtpcchindof=%lf", fMaxTPCchiNdof);
 351   tListSetttings->AddLast(new TObjString(buf));
 352   snprintf(buf, 200, "AliFemtoESDTrackCut.maxsigmatovertex=%lf", fMaxSigmaToVertex);
 353   tListSetttings->AddLast(new TObjString(buf));
 354   if (fMostProbable) {
 355     if (fMostProbable == 2)
 356       snprintf(buf, 200, "AliFemtoESDTrackCut.mostprobable=%s", "Pion");
 357     if (fMostProbable == 3)
 358       snprintf(buf, 200, "AliFemtoESDTrackCut.mostprobable=%s", "Kaon");
 359     if (fMostProbable == 4)
 360       snprintf(buf, 200, "AliFemtoESDTrackCut.mostprobable=%s", "Proton");
 361     tListSetttings->AddLast(new TObjString(buf));
 362   }
 363   return tListSetttings;
 364 }
 365 void AliFemtoESDTrackCut::SetRemoveKinks(const bool& flag)
 366 {
 367   fRemoveKinks = flag;
 368 }
 369
 370                             // electron
 371 // 0.13 - 1.8
 372 // 0       7.594129e-02    8.256141e-03
 373 // 1       -5.535827e-01   8.170825e-02
 374 // 2       1.728591e+00    3.104210e-01
 375 // 3       -2.827893e+00   5.827802e-01
 376 // 4       2.503553e+00    5.736207e-01
 377 // 5       -1.125965e+00   2.821170e-01
 378 // 6       2.009036e-01    5.438876e-02
 379 float AliFemtoESDTrackCut::PidFractionElectron(float mom) const
 380 {
 381   // Provide a parameterized fraction of electrons dependent on momentum
 382   if (mom<0.13) return 0.0;
 383   if (mom>1.8) return 0.0;
 384   return (7.594129e-02
 385           -5.535827e-01*mom
 386           +1.728591e+00*mom*mom
 387           -2.827893e+00*mom*mom*mom
 388           +2.503553e+00*mom*mom*mom*mom
 389           -1.125965e+00*mom*mom*mom*mom*mom
 390           +2.009036e-01*mom*mom*mom*mom*mom*mom);
 391 }
 392
 393 // pion
 394 // 0.13 - 2.0
 395 // 0       1.063457e+00    8.872043e-03
 396 // 1       -4.222208e-01   2.534402e-02
 397 // 2       1.042004e-01    1.503945e-02
 398 float AliFemtoESDTrackCut::PidFractionPion(float mom) const
 399 {
 400   // Provide a parameterized fraction of pions dependent on momentum
 401   if (mom<0.13) return 0.0;
 402   if (mom>2.0) return 0.0;
 403   return ( 1.063457e+00
 404            -4.222208e-01*mom
 405            +1.042004e-01*mom*mom);
 406 }
 407
 408 // kaon
 409 // 0.18 - 2.0
 410 // 0       -7.289406e-02   1.686074e-03
 411 // 1       4.415666e-01    1.143939e-02
 412 // 2       -2.996790e-01   1.840964e-02
 413 // 3       6.704652e-02    7.783990e-03
 414 float AliFemtoESDTrackCut::PidFractionKaon(float mom) const
 415 {
 416   // Provide a parameterized fraction of kaons dependent on momentum
 417   if (mom<0.18) return 0.0;
 418   if (mom>2.0) return 0.0;
 419   return (-7.289406e-02
 420           +4.415666e-01*mom
 421           -2.996790e-01*mom*mom
 422           +6.704652e-02*mom*mom*mom);
 423 }
 424
 425 // proton
 426 // 0.26 - 2.0
 427 // 0       -3.730200e-02   2.347311e-03
 428 // 1       1.163684e-01    1.319316e-02
 429 // 2       8.354116e-02    1.997948e-02
 430 // 3       -4.608098e-02   8.336400e-03
 431 float AliFemtoESDTrackCut::PidFractionProton(float mom) const
 432 {
 433   // Provide a parameterized fraction of protons dependent on momentum
 434   if (mom<0.26) return  0.0;
 435   if (mom>2.0) return 0.0;
 436   return (-3.730200e-02
 437           +1.163684e-01*mom
 438           +8.354116e-02*mom*mom
 439           -4.608098e-02*mom*mom*mom);
 440 }