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 = 10;
  64   //                                   0       1              2              3               4     5      6       7       8         9
  65   const TString atitle[nvar]={"TrackEta","JetPt", "TrackTPCsig", "Log10TrackP", "Log10TrackPt",  "z",  "xi",  "pdg",  "tof",  "hmpid"};
  66   const Int_t nbins[nvar]   ={         4,     20,          2400,          Nx(),             50,   25,    50,      7,      7,        7};
  67   const Double_t xmins[nvar]={         0,      0,             0,        Xmin(),         Xmin(),    0,     0,   -3.5,   -3.5,     -3.5};
  68   const Double_t xmaxs[nvar]={       0.9,    100,           200,        Xmax(),         Xmax(),    1,     6,    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 Bool_t AliIDFFUtils::HMPIDAcceptance(const AliAODTrack *track)
  88 {
  89   //
  90   //check HMPID acceptance
  91   //From S. Pochybova
  92   //
  93
  94   Double_t tEta = TMath::Abs(track->Eta());
  95   Double_t tPhi = track->Phi();
  96   if(tPhi < 0.) tPhi += TMath::TwoPi();
  97   if(tPhi > TMath::TwoPi()) tPhi -= TMath::TwoPi();
  98
  99   if(tEta > 0.46){
 100     return kFALSE;
 101   }
 102
 103   if(tPhi < 0.08 || tPhi > 1.12){
 104     return kFALSE;
 105   }
 106
 107   return kTRUE;
 108 }
 109
 110 Bool_t AliIDFFUtils::HMPIDQA(const AliAODTrack *track)
 111 {
 112   //
 113   //check HMPID PID quality
 114   //From S. Pochybova
 115   //
 116
 117   if(track->GetHMPIDsignal() <= 0.){
 118     return kFALSE;
 119   }
 120
 121   //check track-quality cuts
 122   //dist_(mip-trk)
 123   //track variables
 124   Float_t tX, tY, tTh, tPh;
 125   //mip variables
 126   Float_t mpX, mpY;
 127   Int_t mpQ, mpNph;
 128   track->GetHMPIDtrk(tX, tY, tTh, tPh);
 129   track->GetHMPIDmip(mpX, mpY, mpQ, mpNph);
 130   const Double_t dist = TMath::Sqrt((tX-mpX)*(tX-mpX)+(tY-mpY)*(tY-mpY));
 131
 132   //taking the pass 2 case for the moment
 133   if(dist > 1){
 134     //Printf("Track did not pass the distance cut");
 135     return kFALSE;
 136   }
 137
 138   //cut on charge
 139   //have to check if this also varies with pass
 140   //taking the pass 2 case for the moment
 141   if(mpQ < 120){
 142     //Printf("Track did not pass the MipQ cut");
 143     return kFALSE;
 144   }
 145
 146   return kTRUE;
 147 }
 148
 149 Int_t AliIDFFUtils::HMPIDType(const AliAODTrack * track)
 150 {
 151   //
 152   //return the (locally defined) particle type judged by HMPID
 153   //From S. Pochybova
 154   //
 155
 156   if(!HMPIDAcceptance(track))
 157     return kNOTACCEPTED;
 158
 159   if(!HMPIDQA(track))
 160     return kNOINFO;
 161
 162   Double_t nsigma[]={-999,-999,-999};
 163   nsigma[kPION]     = fPid->NumberOfSigmasHMPID( track, AliPID::kPion);
 164   nsigma[kKAON]     = fPid->NumberOfSigmasHMPID( track, AliPID::kKaon);
 165   nsigma[kPROTON]   = fPid->NumberOfSigmasHMPID( track, AliPID::kProton);
 166
 167   const Double_t inclusion=2;
 168   const Double_t exclusion=3;
 169
 170   if(TMath::Abs(nsigma[kPION])<inclusion && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])>exclusion) return kPION;
 171   if(TMath::Abs(nsigma[kPION])>exclusion && TMath::Abs(nsigma[kKAON])<inclusion && TMath::Abs(nsigma[kPROTON])>exclusion) return kKAON;
 172   if(TMath::Abs(nsigma[kPION])>exclusion && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])<inclusion) return kPROTON;
 173
 174   return kNOTSELECTED;
 175 }
 176
 177 Int_t AliIDFFUtils::TOFType(const AliAODTrack * trackptr, const Int_t tofmode)
 178 {
 179   //
 180   //return the (locally defined) particle type judged by TOF
 181   //
 182
 183   //check kTOFout, kTIME, mismatch
 184   const AliPIDResponse::EDetPidStatus tofstatus =  fPid->CheckPIDStatus(AliPIDResponse::kTOF, trackptr);
 185   if(tofstatus != AliPIDResponse::kDetPidOk)
 186     return kNOINFO;
 187
 188   Double_t nsigma[]={-999,-999,-999};
 189   nsigma[kPION]     = fPid->NumberOfSigmasTOF( trackptr, AliPID::kPion);
 190   nsigma[kKAON]     = fPid->NumberOfSigmasTOF( trackptr, AliPID::kKaon);
 191   nsigma[kPROTON]   = fPid->NumberOfSigmasTOF( trackptr, AliPID::kProton);
 192
 193   Double_t inclusion=-999;
 194   Double_t exclusion=-999;
 195   if(tofmode==1){
 196     inclusion = 1.5;
 197     exclusion = 3;
 198   }
 199   else if(tofmode==2){
 200     inclusion = 2;
 201     exclusion = 3;
 202   }
 203   else if(tofmode==3){
 204     inclusion = 2.5;
 205     exclusion = 3;
 206   }
 207   else{
 208     printf("AliIDFFUtils::TOFType bad tofmode ! %d\n", tofmode); exit(1);
 209   }
 210
 211   if(TMath::Abs(nsigma[kPION])<inclusion && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])>exclusion) return kPION;
 212   if(TMath::Abs(nsigma[kPION])>exclusion && TMath::Abs(nsigma[kKAON])<inclusion && TMath::Abs(nsigma[kPROTON])>exclusion) return kKAON;
 213   if(TMath::Abs(nsigma[kPION])>exclusion && TMath::Abs(nsigma[kKAON])>exclusion && TMath::Abs(nsigma[kPROTON])<inclusion) return kPROTON;
 214
 215   return kNOTSELECTED;
 216 }
 217
 218
 219 void AliIDFFUtils::FillTHn(THnSparseD * hh, Double_t jetpt, const AliAODTrack * track,  AliAODEvent *aodevent, const Int_t tofmode) //AliMCEvent * mcevent)
 220 {
 221   //
 222   //fill variables
 223   //
 224
 225   Int_t mcpdg = kNOINFO;
 226
 227   TClonesArray *tca = dynamic_cast<TClonesArray*>(aodevent->FindListObject(AliAODMCParticle::StdBranchName()));
 228   if(tca){
 229     const Int_t mclabel = TMath::Abs(track->GetLabel());
 230
 231     AliAODMCParticle* gentrack = dynamic_cast<AliAODMCParticle*> (tca->At(mclabel));
 232     if(gentrack){
 233       //printf("null gentrack in AliIDFFUtils::FillTHn mclabel %d\n", mclabel);
 234       mcpdg = PDG2Type(TMath::Abs(gentrack->GetPdgCode()));
 235     }
 236   }
 237
 238   //===========================================================================================================
 239
 240   const Int_t ktof = TOFType(track, tofmode);
 241   const Int_t khmpid = HMPIDType(track);
 242
 243   //===========================================================================================================
 244
 245   const Double_t eps = 1e-6;
 246   const Double_t tracketa = TMath::Abs(track->Eta());
 247   Double_t tracktpc = track->GetTPCsignal();
 248   if(tracktpc>=200)
 249     tracktpc = 200 - eps;
 250   const Double_t tracklogmom = TMath::Log10(track->P());
 251   const Double_t tracklogpt  = TMath::Log10(track->Pt());
 252
 253   Double_t zz = -999;
 254   Double_t xi = -999;
 255   if(jetpt<-990){
 256     jetpt = zz = xi = eps;
 257   }
 258   //from Oliver
 259   else if(track->Pt()>(1-eps)*jetpt && track->Pt()<(1+eps)*jetpt){ // case z=1 : move entry to last histo bin <1
 260     zz  = 1-eps;
 261     xi = eps;
 262   }
 263   else if(jetpt>0){
 264     zz = track->Pt()/jetpt;
 265     xi = TMath::Log(1/zz);
 266   }
 267
 268   //===========================================================================================================
 269
 270   const Double_t vars[]={tracketa, jetpt, tracktpc, tracklogmom, tracklogpt, zz, xi, Double_t(mcpdg), Double_t(ktof), Double_t(khmpid)};
 271
 272   hh->Fill(vars);
 273 }
 274
 275 //_________________________________________________________________________________
 276 Bool_t AliIDFFUtils::TPCCutPIDN(const AliAODTrack * track)
 277 {
 278   //
 279   //TPC Cut kPIDN
 280   //
 281
 282   //assuming this cut is particle type independent, then it is fine
 283   //need further investigation
 284   if(track->GetTPCsignalN()<60){
 285     return kFALSE;
 286   }
 287
 288   return kTRUE;
 289 }
 290
 291 //_________________________________________________________________________________
 292 Bool_t AliIDFFUtils::TPCCutMIGeo(const AliAODTrack * track, const AliVEvent* evt, TTreeStream * streamer)
 293 {
 294   //
 295   //TPC Cut MIGeo
 296   //
 297
 298   Short_t sign = track->Charge();
 299   Double_t xyz[50];
 300   Double_t pxpypz[50];
 301   Double_t cv[100];
 302
 303   track->GetXYZ(xyz);
 304   track->GetPxPyPz(pxpypz);
 305
 306   AliExternalTrackParam * par = new AliExternalTrackParam(xyz, pxpypz, cv, sign);
 307   AliESDtrack dummy;
 308
 309   Double_t varGeom = dummy.GetLengthInActiveZone(par,3,236, evt->GetMagneticField(), 0,0);
 310   Double_t varNcr  = track->GetTPCClusterInfo(3,1);
 311   Double_t varNcls = track->GetTPCsignalN();
 312
 313   Bool_t cutGeom     = varGeom > 1.*(130-5*TMath::Abs(1./track->Pt()));
 314   Bool_t cutNcr      = varNcr  > 0.85*(130-5*TMath::Abs(1./track->Pt()));
 315   Bool_t cutNcls     = varNcls > 0.7*(130-5*TMath::Abs(1./track->Pt()));
 316
 317   Bool_t kout = cutGeom && cutNcr && cutNcls;
 318
 319   if(streamer){
 320     Double_t dedx = track->GetTPCsignal();
 321
 322     (*streamer)<<"tree"<<
 323       "param.="<< par<<
 324
 325       "varGeom="<<varGeom<<
 326       "varNcr="<<varNcr<<
 327       "varNcls="<<varNcls<<
 328
 329       "cutGeom="<<cutGeom<<
 330       "cutNcr="<<cutNcr<<
 331       "cutNcls="<<cutNcls<<
 332
 333       "kout="<<kout<<
 334       "dedx="<<dedx<<
 335
 336       "\n";
 337   }
 338
 339   delete par;
 340
 341   return kout;
 342 }
 343