PWGHF/correlationHF/AliHFAssociatedTrackCuts.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2010, 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 /* $Id$ */
  17
  18 ////////////////////////////////////////////////////////////////////////
  19 //
  20 // Base class for cuts on Associated tracks for HF Correlation analysis
  21 //
  22 // Author: S.Bjelogrlic (Utrecht) sandro.bjelogrlic@cern.ch
  23 ////////////////////////////////////////////////////////////////////////
  24 #include <Riostream.h>
  25 #include "AliHFAssociatedTrackCuts.h"
  26 #include "AliAODPidHF.h"
  27 #include "AliESDtrackCuts.h"
  28 #include "AliESDtrack.h"
  29 #include "AliAODv0.h"
  30 #include "AliAODVertex.h"
  31 #include "AliAODMCParticle.h"
  32 #include "TString.h"
  33
  34 using std::cout;
  35 using std::endl;
  36
  37 ClassImp(AliHFAssociatedTrackCuts)
  38
  39 //--------------------------------------------------------------------------
  40 AliHFAssociatedTrackCuts::AliHFAssociatedTrackCuts():
  41 AliAnalysisCuts(),
  42 fESDTrackCuts(0),
  43 fPidObj(0),
  44
  45 fPoolMaxNEvents(0),
  46 fPoolMinNTracks(0),
  47 fMinEventsToMix(0),
  48 fNzVtxBins(0),
  49 fNzVtxBinsDim(0),
  50 fZvtxBins(0),
  51 fNCentBins(0),
  52 fNCentBinsDim(0),
  53 fCentBins(0),
  54
  55 fNTrackCuts(0),
  56 fAODTrackCuts(0),
  57 fTrackCutsNames(0),
  58 fNvZeroCuts(0),
  59 fAODvZeroCuts(0),
  60 fvZeroCutsNames(0),
  61 fBit(0),
  62 fCharge(0)
  63
  64 {
  65         //
  66         //default constructor
  67         //
  68         //
  69         //default constructor
  70         //
  71
  72 }
  73
  74 //--------------------------------------------------------------------------
  75 AliHFAssociatedTrackCuts::AliHFAssociatedTrackCuts(const char* name, const char* title):
  76 AliAnalysisCuts(name,title),
  77 fESDTrackCuts(0),
  78 fPidObj(0),
  79
  80 fPoolMaxNEvents(0),
  81 fPoolMinNTracks(0),
  82 fMinEventsToMix(0),
  83 fNzVtxBins(0),
  84 fNzVtxBinsDim(0),
  85 fZvtxBins(0),
  86 fNCentBins(0),
  87 fNCentBinsDim(0),
  88 fCentBins(0),
  89
  90 fNTrackCuts(0),
  91 fAODTrackCuts(0),
  92 fTrackCutsNames(0),
  93 fNvZeroCuts(0),
  94 fAODvZeroCuts(0),
  95 fvZeroCutsNames(0),
  96 fBit(0),
  97 fCharge(0)
  98
  99 {
 100         //
 101         //default constructor
 102         //
 103
 104 }
 105 //--------------------------------------------------------------------------
 106 AliHFAssociatedTrackCuts::AliHFAssociatedTrackCuts(const AliHFAssociatedTrackCuts &source) :
 107 AliAnalysisCuts(source),
 108 fESDTrackCuts(source.fESDTrackCuts),
 109 fPidObj(source.fPidObj),
 110
 111 fPoolMaxNEvents(source.fPoolMaxNEvents),
 112 fPoolMinNTracks(source.fPoolMinNTracks),
 113 fMinEventsToMix(source.fMinEventsToMix),
 114 fNzVtxBins(source.fNzVtxBins),
 115 fNzVtxBinsDim(source.fNzVtxBinsDim),
 116 fZvtxBins(source.fZvtxBins),
 117 fNCentBins(source.fNCentBins),
 118 fNCentBinsDim(source.fNCentBinsDim),
 119 fCentBins(source.fCentBins),
 120
 121 fNTrackCuts(source.fNTrackCuts),
 122 fAODTrackCuts(source.fAODTrackCuts),
 123 fTrackCutsNames(source.fTrackCutsNames),
 124 fNvZeroCuts(source.fNvZeroCuts),
 125 fAODvZeroCuts(source.fAODvZeroCuts),
 126 fvZeroCutsNames(source.fvZeroCutsNames),
 127 fBit(source.fBit),
 128 fCharge(source.fCharge)
 129 {
 130         //
 131         // copy constructor
 132         //
 133
 134
 135         AliInfo("AliHFAssociatedTrackCuts::Copy constructor ");
 136         if(source.fESDTrackCuts) AddTrackCuts(source.fESDTrackCuts);
 137         if(source.fAODTrackCuts) SetAODTrackCuts(source.fAODTrackCuts);
 138         if(source.fAODvZeroCuts) SetAODvZeroCuts(source.fAODvZeroCuts);
 139         if(source.fPidObj) SetPidHF(source.fPidObj);
 140 }
 141 //--------------------------------------------------------------------------
 142 AliHFAssociatedTrackCuts &AliHFAssociatedTrackCuts::operator=(const AliHFAssociatedTrackCuts &source)
 143 {
 144         //
 145         // assignment operator
 146         //
 147         if(&source == this) return *this;
 148
 149         AliAnalysisCuts::operator=(source);
 150         fESDTrackCuts=source.fESDTrackCuts;
 151         fPidObj=source.fPidObj;
 152         fNTrackCuts=source.fNTrackCuts;
 153         fAODTrackCuts=source.fAODTrackCuts;
 154         fTrackCutsNames=source.fTrackCutsNames;
 155         fNvZeroCuts=source.fNvZeroCuts;
 156         fAODvZeroCuts=source.fAODvZeroCuts;
 157         fvZeroCutsNames=source.fvZeroCutsNames;
 158         fBit=source.fBit;
 159         fCharge=source.fCharge;
 160
 161         return *this;
 162
 163 }
 164
 165
 166 //--------------------------------------------------------------------------
 167 AliHFAssociatedTrackCuts::~AliHFAssociatedTrackCuts()
 168 {
 169         if(fESDTrackCuts) {delete fESDTrackCuts; fESDTrackCuts = 0;}
 170         if(fPidObj) {delete fPidObj; fPidObj = 0;}
 171         if(fZvtxBins) {delete[] fZvtxBins; fZvtxBins=0;}
 172         if(fCentBins) {delete[] fCentBins; fCentBins=0;}
 173         if(fAODTrackCuts) {delete[] fAODTrackCuts; fAODTrackCuts=0;}
 174         if(fTrackCutsNames) {delete[] fTrackCutsNames; fTrackCutsNames=0;}
 175         if(fAODvZeroCuts){delete[] fAODvZeroCuts; fAODvZeroCuts=0;}
 176         if(fvZeroCutsNames) {delete[] fvZeroCutsNames; fvZeroCutsNames=0;}
 177
 178 }
 179 //--------------------------------------------------------------------------
 180 Bool_t AliHFAssociatedTrackCuts::IsInAcceptance()
 181 {
 182         printf("Careful: method AliHFAssociatedTrackCuts::IsInAcceptance is not implemented yet \n");
 183         return kFALSE;
 184 }
 185 //--------------------------------------------------------------------------
 186 Bool_t AliHFAssociatedTrackCuts::IsHadronSelected(AliAODTrack * track)
 187 {
 188         AliESDtrack esdtrack(track);
 189         if(!fESDTrackCuts->IsSelected(&esdtrack)) return kFALSE;
 190
 191         if(fBit && !track->TestFilterBit(fBit)) return kFALSE; // check the filter bit
 192
 193         return kTRUE;
 194
 195 }
 196
 197 //--------------------------------------------------------------------------
 198 Bool_t AliHFAssociatedTrackCuts::CheckHadronKinematic(Double_t pt, Double_t d0)
 199 {
 200
 201
 202
 203         if(pt < fAODTrackCuts[0]) return kFALSE;
 204         if(pt > fAODTrackCuts[1]) return kFALSE;
 205         if(d0 < fAODTrackCuts[2]) return kFALSE;
 206         if(d0 > fAODTrackCuts[3]) return kFALSE;
 207
 208         return kTRUE;
 209
 210
 211 }
 212 //--------------------------------------------------------------------------
 213
 214 Bool_t AliHFAssociatedTrackCuts::Charge(Short_t charge, AliAODTrack* track)
 215 {// charge is the charge to compare to (for example, a daughter of a D meson)
 216
 217         if(!fCharge) return kTRUE; // if fCharge is set to 0 (no selection on the charge), returns always true
 218         if(track->Charge()!= fCharge*charge) return kFALSE;
 219         return kTRUE;
 220 }
 221
 222 //--------------------------------------------------------------------------
 223 Bool_t AliHFAssociatedTrackCuts::CheckKaonCompatibility(AliAODTrack * track, Bool_t useMc, TClonesArray* mcArray, Int_t method)
 224 {
 225         Bool_t isKaon = kFALSE;
 226
 227         if(useMc) { // on MC
 228                 Int_t hadLabel = track->GetLabel();
 229                 if(hadLabel < 0) return kFALSE;
 230                 AliAODMCParticle* hadron = dynamic_cast<AliAODMCParticle*>(mcArray->At(hadLabel));
 231                 Int_t pdg = TMath::Abs(hadron->GetPdgCode());
 232                 if (pdg == 321) isKaon = kTRUE;
 233         }
 234
 235         if(!useMc) { // on DATA
 236           switch(method) {
 237           case(1): {
 238                 Bool_t isKTPC=kFALSE;
 239                 Bool_t isPiTPC=kFALSE;
 240                 Bool_t isPTPC=kFALSE;
 241                 Bool_t isKTOF=kFALSE;
 242                 Bool_t isPiTOF=kFALSE;
 243                 Bool_t isPTOF=kFALSE;
 244
 245                 Bool_t KaonHyp = kFALSE;
 246                 Bool_t PionHyp = kFALSE;
 247                 Bool_t ProtonHyp = kFALSE;
 248
 249                 if(fPidObj->CheckStatus(track,"TOF")) {
 250                         isKTOF=fPidObj->IsKaonRaw(track,"TOF");
 251                         isPiTOF=fPidObj->IsPionRaw(track,"TOF");
 252                         isPTOF=fPidObj->IsProtonRaw(track,"TOF");
 253                 }
 254                 if(fPidObj->CheckStatus(track,"TPC")){
 255                         isKTPC=fPidObj->IsKaonRaw(track,"TPC");
 256                         isPiTPC=fPidObj->IsPionRaw(track,"TPC");
 257                         isPTPC=fPidObj->IsProtonRaw(track,"TPC");
 258                 }
 259
 260                 if (isKTOF && isKTPC) KaonHyp = kTRUE;
 261                 if (isPiTOF && isPiTPC) PionHyp = kTRUE;
 262                 if (isPTOF && isPTPC) ProtonHyp = kTRUE;
 263
 264                 if(KaonHyp && !PionHyp && !ProtonHyp) isKaon = kTRUE;
 265                 break;
 266               }
 267       case(2): {
 268                 if(fPidObj->MakeRawPid(track,3)>=1) isKaon = kTRUE;
 269                 break;
 270               }
 271       }
 272         }
 273
 274         return isKaon;
 275
 276 }
 277 //--------------------------------------------------------------------------
 278 Bool_t AliHFAssociatedTrackCuts::IsKZeroSelected(AliAODv0 *vzero, AliAODVertex *vtx1)
 279 {
 280
 281         if(vzero->DcaV0Daughters()>fAODvZeroCuts[0]) return kFALSE;
 282         if(vzero->Chi2V0()>fAODvZeroCuts[1]) return kFALSE;
 283         if(vzero->DecayLength(vtx1) < fAODvZeroCuts[2]) return kFALSE;
 284         if(vzero->DecayLength(vtx1) > fAODvZeroCuts[3]) return kFALSE;
 285         if(vzero->OpenAngleV0() > fAODvZeroCuts[4]) return kFALSE;
 286         if(vzero->Pt() < fAODvZeroCuts[5]) return kFALSE;
 287         if(TMath::Abs(vzero->Eta()) > fAODvZeroCuts[6]) return kFALSE;
 288
 289
 290         return kTRUE;
 291 }
 292 //--------------------------------------------------------------------------
 293 Int_t AliHFAssociatedTrackCuts::IsMCpartFromHF(Int_t label, TClonesArray*mcArray){
 294   // Check origin in MC
 295
 296   AliAODMCParticle* mcParticle;
 297   Int_t pdgCode = -1;
 298   if (label<0) return 0;
 299   Bool_t isCharmy = kFALSE;
 300   Bool_t isBeauty = kFALSE;
 301   Bool_t isD = kFALSE;
 302   Bool_t isB = kFALSE;
 303
 304   while(pdgCode!=2212){ // loops back to the collision to check the particle source
 305
 306     mcParticle = dynamic_cast<AliAODMCParticle*>(mcArray->At(label));
 307     if(!mcParticle) {AliError("NO MC PARTICLE"); break;}
 308     pdgCode =  TMath::Abs(mcParticle->GetPdgCode());
 309
 310     label = mcParticle->GetMother();
 311
 312
 313     if((pdgCode>=400 && pdgCode <500) || (pdgCode>=4000 && pdgCode<5000 )) isD = kTRUE;
 314     if((pdgCode>=500 && pdgCode <600) || (pdgCode>=5000 && pdgCode<6000 )) {isD = kFALSE; isB = kTRUE;}
 315
 316
 317     if(pdgCode == 4) isCharmy = kTRUE;
 318     if(pdgCode == 5) {isBeauty = kTRUE; isCharmy = kFALSE;}
 319     if(label<0) break;
 320
 321   }
 322   if (isCharmy && isD) return 44; // the first 4(5) indicates the charm/beauty quark, the second the charmy meson
 323   if (isBeauty && isD) return 54;
 324   if (isBeauty && isB) return 55;
 325   return 0;
 326 }
 327
 328 //--------------------------------------------------------------------------
 329 Bool_t AliHFAssociatedTrackCuts::InvMassDstarRejection(AliAODRecoDecayHF2Prong* d, AliAODTrack *track, Int_t hypD0) const {
 330         //
 331         // Calculates invmass of track+D0 and rejects if compatible with D*
 332         // (to remove pions from D*)
 333         //
 334         Double_t nsigma = 3.;
 335
 336         Double_t mD0, mD0bar;
 337         d->InvMassD0(mD0,mD0bar);
 338
 339         Double_t invmassDstar1 = 0, invmassDstar2 = 0;
 340         Double_t e1Pi = d->EProng(0,211), e2K = d->EProng(1,321); //hyp 1 (pi,K) - D0
 341         Double_t e1K = d->EProng(0,321), e2Pi = d->EProng(1,211); //hyp 2 (K,pi) - D0bar
 342         Double_t psum2 = (d->Px()+track->Px())*(d->Px()+track->Px())
 343         +(d->Py()+track->Py())*(d->Py()+track->Py())
 344         +(d->Pz()+track->Pz())*(d->Pz()+track->Pz());
 345
 346         switch(hypD0) {
 347                 case 1:
 348                         invmassDstar1 = TMath::Sqrt(pow(e1Pi+e2K+track->E(0.1396),2.)-psum2);
 349                         if ((TMath::Abs(invmassDstar1-mD0)-0.14543) < nsigma*800.*pow(10.,-6.)) return kFALSE;
 350                         break;
 351                 case 2:
 352                         invmassDstar2 = TMath::Sqrt(pow(e2Pi+e1K+track->E(0.1396),2.)-psum2);
 353                         if ((TMath::Abs(invmassDstar2-mD0bar)-0.14543) < nsigma*800.*pow(10.,-6.)) return kFALSE;
 354                         break;
 355                 case 3:
 356                         invmassDstar1 = TMath::Sqrt(pow(e1Pi+e2K+track->E(0.1396),2.)-psum2);
 357                         invmassDstar2 = TMath::Sqrt(pow(e2Pi+e1K+track->E(0.1396),2.)-psum2);
 358                         if ((TMath::Abs(invmassDstar1-mD0)-0.14543) < nsigma*800.*pow(10.,-6.)) return kFALSE;
 359                         if ((TMath::Abs(invmassDstar2-mD0bar)-0.14543) < nsigma*800.*pow(10.,-6.)) return kFALSE;
 360                         break;
 361         }
 362
 363         return kTRUE;
 364 }
 365 //________________________________________________________
 366 void AliHFAssociatedTrackCuts::SetAODTrackCuts(Float_t *cutsarray)
 367 {
 368         cout << "Check 1" << endl;
 369         if(!fAODTrackCuts) fAODTrackCuts = new Float_t[fNTrackCuts];
 370         cout << "Check 2" << endl;
 371         for(Int_t i =0; i<fNTrackCuts; i++){
 372                 cout << "Check 2." << i << endl;
 373                 fAODTrackCuts[i] = cutsarray[i];
 374         }
 375         cout << "Check 3" << endl;
 376         SetTrackCutsNames();
 377         cout << "Check 4" << endl;
 378         return;
 379 }
 380 //________________________________________________________
 381 void AliHFAssociatedTrackCuts::SetTrackCutsNames(/*TString *namearray*/){
 382
 383         fTrackCutsNames = new TString[4];
 384         fTrackCutsNames[0]= "associated track:: pt min [GeV/c]................: ";
 385         fTrackCutsNames[1]= "associated track:: pt max [GeV/c]................: ";
 386         fTrackCutsNames[2]= "associated track:: d0 min [cm]...................: ";
 387         fTrackCutsNames[3]= "associated track:: d0 max [cm]...................: ";
 388
 389
 390
 391         return;
 392 }
 393 //--------------------------------------------------------------------------
 394 void AliHFAssociatedTrackCuts::SetAODvZeroCuts(Float_t *cutsarray)
 395 {
 396
 397
 398         if(!fAODvZeroCuts) fAODvZeroCuts = new Float_t[fNvZeroCuts];
 399         for(Int_t i =0; i<fNvZeroCuts; i++){
 400                 fAODvZeroCuts[i] = cutsarray[i];
 401         }
 402         SetvZeroCutsNames();
 403         return;
 404 }
 405 //--------------------------------------------------------------------------
 406 void AliHFAssociatedTrackCuts::SetvZeroCutsNames(/*TString *namearray*/){
 407
 408         fvZeroCutsNames = new TString[7];
 409         fvZeroCutsNames[0] = "vZero:: max DCA between two daughters [cm].......: ";
 410         fvZeroCutsNames[1] = "vZero:: max fit Chi Square between two daughters.: ";
 411         fvZeroCutsNames[2] = "vZero:: min decay length [cm]....................: ";
 412         fvZeroCutsNames[3] = "vZero:: max decay length [cm]....................: ";
 413         fvZeroCutsNames[4] = "vZero:: max opening angle between daughters [rad]: ";
 414         fvZeroCutsNames[5] = "vZero:: pt min [Gev/c]...........................: ";
 415         fvZeroCutsNames[6] = "vZero:: |Eta| range <............................: ";
 416
 417
 418         return;
 419 }
 420
 421 //--------------------------------------------------------------------------
 422 void AliHFAssociatedTrackCuts::PrintAll()
 423 {
 424         printf("\nCuts for the associated track: \n \n");
 425
 426         printf("ITS Refit........................................: %s\n",fESDTrackCuts->GetRequireITSRefit() ? "Yes" : "No");
 427         printf("TPC Refit........................................: %s\n",fESDTrackCuts->GetRequireTPCRefit() ? "Yes" : "No");
 428         printf("ITS SA...........................................: %s\n",fESDTrackCuts->GetRequireITSStandAlone() ? "Yes" : "No");
 429         printf("TPC SA...........................................: %s\n",fESDTrackCuts->GetRequireTPCStandAlone() ? "Yes" : "No");
 430         printf("Min number of ITS clusters.......................: %d\n",fESDTrackCuts->GetMinNClustersITS());
 431         printf("Min number of TPC clusters.......................: %d\n",fESDTrackCuts->GetMinNClusterTPC());
 432         Int_t spd = fESDTrackCuts->GetClusterRequirementITS(AliESDtrackCuts::kSPD);
 433         if(spd==0) cout <<  "SPD..............................................: kOff"  << endl;
 434         if(spd==1) cout <<  "SPD..............................................: kNone"  << endl;
 435         if(spd==2) cout <<  "SPD..............................................: kAny"  << endl;
 436         if(spd==3) cout <<  "SPD..............................................: kFirst"  << endl;
 437         if(spd==4) cout <<  "SPD..............................................: kOnlyFirst"  << endl;
 438         if(spd==5) cout <<  "SPD..............................................: kSecond"  << endl;
 439         if(spd==6) cout <<  "SPD..............................................: kOnlySecond"  << endl;
 440         if(spd==7) cout <<  "SPD..............................................: kBoth"  << endl;
 441
 442         cout <<  "Filter Bit.......................................: " << fBit  << endl;
 443         cout <<  "Charge...........................................: " << fCharge  << endl;
 444
 445         for(Int_t j=0;j<fNTrackCuts;j++){
 446                 cout << fTrackCutsNames[j] << fAODTrackCuts[j] << endl;
 447         }
 448
 449         printf("\nCuts for the K0 candidates: \n \n");
 450         for(Int_t k=0;k<fNvZeroCuts;k++){
 451                 cout << fvZeroCutsNames[k] <<  fAODvZeroCuts[k] << endl;
 452         }
 453         cout << " " << endl;
 454         PrintPoolParameters();
 455
 456 }
 457
 458 //--------------------------------------------------------------------------
 459 void AliHFAssociatedTrackCuts::PrintPoolParameters()
 460 {
 461         printf("\nEvent Pool settings: \n \n");
 462
 463         printf("Number of zVtx Bins: %d\n", fNzVtxBins);
 464         printf("\nzVtx Bins:\n");
 465         //Double_t zVtxbinLims[fNzVtxBins+1] = fNzVtxBins;
 466         for(Int_t k=0; k<fNzVtxBins; k++){
 467                 printf("Bin %d..............................................: %.1f - %.1f cm\n", k, fZvtxBins[k], fZvtxBins[k+1]);
 468         }
 469
 470         printf("\nNumber of Centrality(multiplicity) Bins: %d\n", fNCentBins);
 471         printf("\nCentrality(multiplicity) Bins:\n");
 472         for(Int_t k=0; k<fNCentBins; k++){
 473                 printf("Bin %d..............................................: %.1f - %.1f\n", k, fCentBins[k], fCentBins[k+1]);
 474         }
 475
 476
 477
 478 }
 479
 480