PWGJE/UserTasks/AliIDFFUtils.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //Utils for identified fragmentation function (IDFF) analysis
  17 //Author: Xianguo Lu (xianguo.lu@cern.ch)
  18
  19 #include "AliIDFFUtils.h"
  20 #include "AliAODEvent.h"
  21 #include "AliAODMCParticle.h"
  22 #include "AliESDtrack.h"
  23
  24 #include "TClonesArray.h"
  25
  26 ClassImp(AliIDFFUtils);
  27
  28 AliPIDResponse * AliIDFFUtils::fPid=0x0;
  29
  30 Int_t AliIDFFUtils::PDG2Type(const Int_t pdg)
  31 {
  32   //
  33   //conversion from pdg code to local definition of particle type
  34   //
  35
  36   Int_t itype = kNOTSELECTED;
  37   switch(pdg){
  38   case 11:
  39     itype = kELECTRON;
  40     break;
  41   case 211:
  42     itype = kPION;
  43     break;
  44   case 2212:
  45     itype = kPROTON;
  46     break;
  47   case 321:
  48     itype = kKAON;
  49     break;
  50   default:
  51     itype = kNOTSELECTED;
  52     break;
  53   }
  54   return itype;
  55 }
  56
  57 THnSparseD *AliIDFFUtils::GetTHn(const TString name)
  58 {
  59   //
  60   //get THnSparseD
  61   //
  62
  63   const Int_t nvar = 11;
  64   //                                   0       1              2              3               4     5      6       7       8         9     10
  65   const TString atitle[nvar]={"TrackEta","JetPt", "TrackTPCsig", "Log10TrackP", "Log10TrackPt",  "z",  "xi",  "pdg",  "comb",   "tof",  "tpc"};
  66   const Int_t nbins[nvar]   ={         4,     15,          1200,          Nx(),             50,   30,    60,      7,      7,        7,      7};
  67   const Double_t xmins[nvar]={         0,      5,             0,        Xmin(),         Xmin(),    0,     0,   -3.5,   -3.5,     -3.5,   -3.5};
  68   const Double_t xmaxs[nvar]={       0.9,     20,           200,        Xmax(),         Xmax(),  1.2,     6,    3.5,    3.5,      3.5,    3.5};
  69
  70   THnSparseD * hh = new THnSparseD(name,"", nvar, nbins, xmins, xmaxs);
  71   for(Int_t ia=0; ia<nvar; ia++){
  72     hh->GetAxis(ia)->SetTitle(atitle[ia]);
  73   }
  74
  75   TAxis * ax = 0x0;
  76   Int_t nb = 0;
  77
  78   //0: eta bin
  79   ax = hh->GetAxis(0);
  80   nb = ax->GetNbins();
  81   const Double_t etas[]={0, 0.2, 0.4, 0.6, 0.9};
  82   ax->Set(nb, etas);
  83
  84   return hh;
  85 }
  86
  87 /*
  88 Bool_t AliIDFFUtils::HMPIDAcceptance(const AliAODTrack *track)
  89 {
  90   //
  91   //check HMPID acceptance
  92   //From S. Pochybova
  93   //
  94
  95   Double_t tEta = TMath::Abs(track->Eta());
  96   Double_t tPhi = track->Phi();
  97   if(tPhi < 0.) tPhi += TMath::TwoPi();
  98   if(tPhi > TMath::TwoPi()) tPhi -= TMath::TwoPi();
  99
 100   if(tEta > 0.46){
 101     return kFALSE;
 102   }
 103
 104   if(tPhi < 0.08 || tPhi > 1.12){
 105     return kFALSE;
 106   }
 107
 108   return kTRUE;
 109 }
 110
 111 Bool_t AliIDFFUtils::HMPIDQA(const AliAODTrack *track)
 112 {
 113   //
 114   //check HMPID PID quality
 115   //From S. Pochybova
 116   //
 117
 118   if(track->GetHMPIDsignal() <= 0.){
 119     return kFALSE;
 120   }
 121
 122   //check track-quality cuts
 123   //dist_(mip-trk)
 124   //track variables
 125   Float_t tX, tY, tTh, tPh;
 126   //mip variables
 127   Float_t mpX, mpY;
 128   Int_t mpQ, mpNph;
 129   track->GetHMPIDtrk(tX, tY, tTh, tPh);
 130   track->GetHMPIDmip(mpX, mpY, mpQ, mpNph);
 131   const Double_t dist = TMath::Sqrt((tX-mpX)*(tX-mpX)+(tY-mpY)*(tY-mpY));
 132
 133   //taking the pass 2 case for the moment
 134   if(dist > 1){
 135     //Printf("Track did not pass the distance cut");
 136     return kFALSE;
 137   }
 138
 139   //cut on charge
 140   //have to check if this also varies with pass
 141   //taking the pass 2 case for the moment
 142   if(mpQ < 120){
 143     //Printf("Track did not pass the MipQ cut");
 144     return kFALSE;
 145   }
 146
 147   return kTRUE;
 148 }
 149
 150 Int_t AliIDFFUtils::HMPIDType(const AliAODTrack * track)
 151 {
 152   //
 153   //return the (locally defined) particle type judged by HMPID
 154   //From S. Pochybova
 155   //
 156
 157   if(!HMPIDAcceptance(track))
 158     return kNOTACCEPTED;
 159
 160   if(!HMPIDQA(track))
 161     return kNOINFO;
 162
 163   Double_t nsigma[]={-999,-999,-999};
 164   nsigma[kPION]     = fPid->NumberOfSigmasHMPID( track, AliPID::kPion);
 165   nsigma[kKAON]     = fPid->NumberOfSigmasHMPID( track, AliPID::kKaon);
 166   nsigma[kPROTON]   = fPid->NumberOfSigmasHMPID( track, AliPID::kProton);
 167
 168   const Double_t inclusion=2;
 169   const Double_t exclusion=3;
 170
 171   if(TMath::Abs(nsigma[kPION])<inclusion && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])>exclusion) return kPION;
 172   if(TMath::Abs(nsigma[kPION])>exclusion && TMath::Abs(nsigma[kKAON])<inclusion && TMath::Abs(nsigma[kPROTON])>exclusion) return kKAON;
 173   if(TMath::Abs(nsigma[kPION])>exclusion && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])<inclusion) return kPROTON;
 174
 175   return kNOTSELECTED;
 176 }
 177 */
 178
 179
 180 Int_t AliIDFFUtils::TPCType(const AliAODTrack * trackptr)
 181 {
 182   //
 183   //return the (locally defined) particle type judged by TPC
 184   //use tofmode for TPC mode
 185   //
 186
 187   const AliPIDResponse::EDetPidStatus tpcstatus =  fPid->CheckPIDStatus(AliPIDResponse::kTPC, trackptr);
 188   if(tpcstatus != AliPIDResponse::kDetPidOk)
 189     return kNOINFO;
 190
 191   Double_t nsigma[]={-999,-999,-999, -999};
 192   nsigma[kPION]     = fPid->NumberOfSigmasTPC( trackptr, AliPID::kPion);
 193   nsigma[kKAON]     = fPid->NumberOfSigmasTPC( trackptr, AliPID::kKaon);
 194   nsigma[kPROTON]   = fPid->NumberOfSigmasTPC( trackptr, AliPID::kProton);
 195   nsigma[kELECTRON] = fPid->NumberOfSigmasTPC( trackptr, AliPID::kElectron);
 196
 197   //so that the effective region is really low momentum
 198   const Double_t inclusion=5;
 199   const Double_t exclusion=5;
 200
 201   //don't destroy TPC signal shape below 120
 202   const Double_t maxsig = 150;
 203   if(trackptr->GetTPCsignal()> maxsig){
 204     if(TMath::Abs(nsigma[kPION])<inclusion && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])>exclusion && TMath::Abs(nsigma[kELECTRON])>exclusion) return kPION;
 205     if(TMath::Abs(nsigma[kPION])>exclusion && TMath::Abs(nsigma[kKAON])<inclusion && TMath::Abs(nsigma[kPROTON])>exclusion && TMath::Abs(nsigma[kELECTRON])>exclusion) return kKAON;
 206     if(TMath::Abs(nsigma[kPION])>exclusion && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])<inclusion && TMath::Abs(nsigma[kELECTRON])>exclusion) return kPROTON;
 207     if(TMath::Abs(nsigma[kPION])>exclusion && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])>exclusion && TMath::Abs(nsigma[kELECTRON])<inclusion) return kELECTRON;
 208   }
 209
 210   return kNOTSELECTED;
 211 }
 212
 213 Int_t AliIDFFUtils::TOFType(const AliAODTrack * trackptr, const Int_t tofmode)
 214 {
 215   //
 216   //return the (locally defined) particle type judged by TOF
 217   //
 218
 219   //check kTOFout, kTIME, mismatch
 220   const AliPIDResponse::EDetPidStatus tofstatus =  fPid->CheckPIDStatus(AliPIDResponse::kTOF, trackptr);
 221   if(tofstatus != AliPIDResponse::kDetPidOk)
 222     return kNOINFO;
 223
 224   Double_t nsigma[]={-999,-999,-999, -999};
 225   nsigma[kPION]     = fPid->NumberOfSigmasTOF( trackptr, AliPID::kPion);
 226   nsigma[kKAON]     = fPid->NumberOfSigmasTOF( trackptr, AliPID::kKaon);
 227   nsigma[kPROTON]   = fPid->NumberOfSigmasTOF( trackptr, AliPID::kProton);
 228   nsigma[kELECTRON] = fPid->NumberOfSigmasTOF( trackptr, AliPID::kElectron);
 229
 230   Double_t inclusion=-999;
 231   Double_t exclusion=-999;
 232   if(tofmode == 1){
 233     inclusion = 2;
 234     exclusion = 2;
 235   }
 236   else if(tofmode == 2){
 237     inclusion = 2;
 238     exclusion = 3;
 239   }
 240   else if(tofmode == 3){
 241     inclusion = 3;
 242     exclusion = 3;
 243   }
 244   else if(tofmode == 4){
 245     inclusion = 3;
 246     exclusion = 4;
 247   }
 248   else{
 249     printf("AliIDFFUtils::TOFType bad tofmode ! %d\n", tofmode); exit(1);
 250   }
 251
 252   const Bool_t kpassEle = kTRUE;
 253   /*
 254   const Double_t cutEle = 1.5;
 255   //tofmode = 1x then require electron exclusion cut
 256   if( tofmode == 4 ){
 257     if(TMath::Abs(nsigma[kELECTRON])> cutEle ){
 258       kpassEle = kTRUE;
 259     }
 260     else{
 261       kpassEle = kFALSE;
 262     }
 263   }
 264   */
 265
 266   //cut on electron for pion because the precision of pion is good and the contamination of electron can not be ignored
 267   //+1 exclusion sigma in electron/pion to enforce better purity, otherwise not only pion, but also kaon is bias for jet pt 5-10 GeV/c
 268   if(TMath::Abs(nsigma[kPION])<inclusion     && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])>exclusion && kpassEle) return kPION;
 269   if(TMath::Abs(nsigma[kPION])>exclusion     && TMath::Abs(nsigma[kKAON])<inclusion && TMath::Abs(nsigma[kPROTON])>exclusion && kpassEle) return kKAON;
 270   if(TMath::Abs(nsigma[kPION])>exclusion     && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])<inclusion && kpassEle) return kPROTON;
 271   if(TMath::Abs(nsigma[kPION])>exclusion+0.5 && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])>exclusion && TMath::Abs(nsigma[kELECTRON])<inclusion) return kELECTRON;
 272
 273   return kNOTSELECTED;
 274 }
 275
 276 Int_t AliIDFFUtils::CombineTPCTOF(const Int_t ktpc, const Int_t ktof)
 277 {
 278   //tpc and tof, if <0 only noinfo or notselected
 279   if(ktpc == ktof)
 280     return ktpc;
 281   else if(ktpc < 0 && ktof >= 0 )
 282     return ktof;
 283   else if(ktof < 0 && ktpc >= 0)
 284     return ktpc;
 285   else
 286     return kNOTACCEPTED;
 287 }
 288
 289 void AliIDFFUtils::FillTHn(THnSparseD * hh, Double_t jetpt, const AliAODTrack * track,  AliAODEvent *aodevent, const Int_t tofmode) //AliMCEvent * mcevent)
 290 {
 291   //
 292   //fill variables
 293   //
 294
 295   Int_t mcpdg = kNOINFO;
 296
 297   TClonesArray *tca = dynamic_cast<TClonesArray*>(aodevent->FindListObject(AliAODMCParticle::StdBranchName()));
 298   if(tca){
 299     const Int_t mclabel = TMath::Abs(track->GetLabel());
 300
 301     AliAODMCParticle* gentrack = dynamic_cast<AliAODMCParticle*> (tca->At(mclabel));
 302     if(gentrack){
 303       //printf("null gentrack in AliIDFFUtils::FillTHn mclabel %d\n", mclabel);
 304       mcpdg = PDG2Type(TMath::Abs(gentrack->GetPdgCode()));
 305     }
 306   }
 307
 308   //===========================================================================================================
 309
 310   //use tofmode for tpcmode
 311   const Int_t ktpc = TPCType(track);
 312   const Int_t ktof = TOFType(track, tofmode);
 313
 314   //fake kcomb by pretpc+pretof
 315   const Int_t kcomb = CombineTPCTOF(ktpc, ktof);
 316
 317   //const Int_t khmpid = HMPIDType(track);
 318
 319   //===========================================================================================================
 320
 321   const Double_t eps = 1e-6;
 322   const Double_t tracketa = TMath::Abs(track->Eta());
 323   Double_t tracktpc = track->GetTPCsignal();
 324   if(tracktpc>=200)
 325     tracktpc = 200 - eps;
 326   const Double_t tracklogmom = TMath::Log10(track->P());
 327   const Double_t tracklogpt  = TMath::Log10(track->Pt());
 328
 329   Double_t zz = -999;
 330   Double_t xi = -999;
 331   if(jetpt<-990){
 332     jetpt = zz = xi = eps;
 333   }
 334   //from Oliver
 335   else if(track->Pt()>(1-eps)*jetpt && track->Pt()<(1+eps)*jetpt){ // case z=1 : move entry to last histo bin <1
 336     zz  = 1-eps;
 337     xi = eps;
 338   }
 339   else if(jetpt>0){
 340     zz = track->Pt()/jetpt;
 341     xi = TMath::Log(1/zz);
 342   }
 343
 344   //===========================================================================================================
 345
 346   const Double_t vars[]={tracketa, jetpt, tracktpc, tracklogmom, tracklogpt, zz, xi, (Double_t)mcpdg, (Double_t)kcomb, (Double_t)ktof, (Double_t)ktpc};
 347
 348   hh->Fill(vars);
 349 }
 350
 351 //_________________________________________________________________________________
 352 Bool_t AliIDFFUtils::TPCCutPIDN(const AliAODTrack * track)
 353 {
 354   //
 355   //TPC Cut kPIDN
 356   //
 357
 358   //assuming this cut is particle type independent, then it is fine
 359   //need further investigation
 360   if(track->GetTPCsignalN()<60){
 361     return kFALSE;
 362   }
 363
 364   return kTRUE;
 365 }
 366
 367 //_________________________________________________________________________________
 368 Bool_t AliIDFFUtils::TPCCutMIGeo(const AliAODTrack * track, const AliVEvent* evt, TTreeStream * streamer)
 369 {
 370   //
 371   //TPC Cut MIGeo
 372   //
 373
 374   Short_t sign = track->Charge();
 375   Double_t xyz[50];
 376   Double_t pxpypz[50];
 377   Double_t cv[100];
 378
 379   track->GetXYZ(xyz);
 380   track->GetPxPyPz(pxpypz);
 381
 382   AliExternalTrackParam * par = new AliExternalTrackParam(xyz, pxpypz, cv, sign);
 383   AliESDtrack dummy;
 384
 385   Double_t varGeom = dummy.GetLengthInActiveZone(par,3,236, evt->GetMagneticField(), 0,0);
 386   Double_t varNcr  = track->GetTPCClusterInfo(3,1);
 387   Double_t varNcls = track->GetTPCsignalN();
 388
 389   Bool_t cutGeom     = varGeom > 1.*(130-5*TMath::Abs(1./track->Pt()));
 390   Bool_t cutNcr      = varNcr  > 0.85*(130-5*TMath::Abs(1./track->Pt()));
 391   Bool_t cutNcls     = varNcls > 0.7*(130-5*TMath::Abs(1./track->Pt()));
 392
 393   Bool_t kout = cutGeom && cutNcr && cutNcls;
 394
 395   if(streamer){
 396     Double_t dedx = track->GetTPCsignal();
 397
 398     (*streamer)<<"tree"<<
 399       "param.="<< par<<
 400
 401       "varGeom="<<varGeom<<
 402       "varNcr="<<varNcr<<
 403       "varNcls="<<varNcls<<
 404
 405       "cutGeom="<<cutGeom<<
 406       "cutNcr="<<cutNcr<<
 407       "cutNcls="<<cutNcls<<
 408
 409       "kout="<<kout<<
 410       "dedx="<<dedx<<
 411
 412       "\n";
 413   }
 414
 415   delete par;
 416
 417   return kout;
 418 }
 419