PWGCF/Correlations/DPhi/DiHadronPID/AliTrackDiHadronPID.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2008, 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 // -----------------------------------------------------------------------
  17 //  Track class for the DiHadronPID analysis.
  18 // -----------------------------------------------------------------------
  19 //  Author: Misha Veldhoen (misha.veldhoen@cern.ch)
  20
  21 #include "AliTrackDiHadronPID.h"
  22
  23 #include "AliAODVertex.h"
  24 #include "AliPID.h"
  25 #include "AliTPCPIDResponse.h"
  26
  27 ClassImp(AliTrackDiHadronPID);
  28
  29 Double_t AliTrackDiHadronPID::fSigmaTOFStd = 80.;       // Should perhaps be replaced with a
  30 Double_t AliTrackDiHadronPID::fSigmaTPCStd = 3.5;   // function later.
  31
  32 // -----------------------------------------------------------------------
  33 AliTrackDiHadronPID::AliTrackDiHadronPID():
  34         TObject(),
  35         fAODTrack(0x0),
  36         fAODGlobalTrack(0x0),
  37         fAODEvent(0x0),
  38         fAODMCParticle(0x0),
  39         fPIDResponse(0x0),
  40         fBasicInfoAvailable(kFALSE),
  41         fFlagsAvailable(kFALSE),
  42         fDCAInfoAvailable(kFALSE),
  43         fITSInfoAvailable(kFALSE),
  44         fTPCInfoAvailable(kFALSE),
  45         fTOFInfoAvailable(kFALSE),
  46         fMCInfoAvailable(kFALSE),
  47         fPt(-999.),
  48         fEta(-999.),
  49         fPhi(-999.),
  50         fFlags(0),
  51         fFilterMap(0),
  52         fID(0),
  53         fLabel(0),
  54         fCharge(0),
  55         fNclsTPC(-999),
  56         fDCAz(-999.),
  57         fDCAxy(-999.),
  58         fTOFsignal(-999.),
  59         fTOFMatchingStatus(-1),
  60         fTPCsignal(-999.),
  61         fTPCmomentum(-999.),
  62         fITSClusterMap(0),
  63         fMCPt(-999.),
  64         fMCEta(-999.),
  65         fMCPhi(-999.),
  66         fMCY(-999.),
  67         fPdgCode(0),
  68         fIsPhysicalPrimary(kFALSE),
  69         fIsSecondaryFromWeakDecay(kFALSE),
  70         fIsSecondaryFromMaterial(kFALSE),
  71         fDebug(0)
  72
  73 {
  74
  75         //
  76         // Default Constructor.
  77         //
  78
  79         if (fDebug > 2) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
  80
  81         for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
  82                 fTOFsignalMinusExpected[iSpecies] = -999.;
  83                 fTOFNsigma[iSpecies] = -999.;
  84                 fTPCsignalMinusExpected[iSpecies] = -999.;
  85                 fTPCNsigma[iSpecies] = -999.;
  86                 fY[iSpecies] = -999.;
  87         }
  88
  89         for (Int_t iITSlayer = 0; iITSlayer < 6; iITSlayer++) {
  90                 fITSHits[iITSlayer] = kFALSE;
  91         }
  92
  93         for (Int_t iN = 0; iN < 3; ++iN) {
  94                 fTOFLabel[iN] = -1;     // Same convention as in ESDs
  95         }
  96
  97 }
  98
  99 // -----------------------------------------------------------------------
 100 AliTrackDiHadronPID::AliTrackDiHadronPID(AliAODTrack* track, AliAODTrack* globaltrack, AliAODMCParticle* mcparticle, AliPIDResponse* pidresponse):
 101         TObject(),
 102         fAODTrack(0x0),
 103         fAODGlobalTrack(0x0),
 104         fAODEvent(0x0),
 105         fAODMCParticle(0x0),
 106         fPIDResponse(0x0),
 107         fBasicInfoAvailable(kFALSE),
 108         fFlagsAvailable(kFALSE),
 109         fDCAInfoAvailable(kFALSE),
 110         fITSInfoAvailable(kFALSE),
 111         fTPCInfoAvailable(kFALSE),
 112         fTOFInfoAvailable(kFALSE),
 113         fMCInfoAvailable(kFALSE),
 114         fPt(-999.),
 115         fEta(-999.),
 116         fPhi(-999.),
 117         fFlags(0),
 118         fFilterMap(0),
 119         fID(0),
 120         fLabel(0),
 121         fCharge(0),
 122         fNclsTPC(-999),
 123         fDCAz(-999.),
 124         fDCAxy(-999.),
 125         fTOFsignal(-999.),
 126         fTOFMatchingStatus(-1),
 127         fTPCsignal(-999.),
 128         fTPCmomentum(-999.),
 129         fITSClusterMap(0),
 130         fMCPt(-999.),
 131         fMCEta(-999.),
 132         fMCPhi(-999.),
 133         fMCY(-999.),
 134         fPdgCode(0),
 135         fIsPhysicalPrimary(kFALSE),
 136         fIsSecondaryFromWeakDecay(kFALSE),
 137         fIsSecondaryFromMaterial(kFALSE),
 138         fDebug(0)
 139 {
 140
 141         //
 142         // Constructor.
 143         //
 144
 145         if (fDebug > 2) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 146
 147         for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
 148                 fTOFsignalMinusExpected[iSpecies] = -999.;
 149                 fTOFNsigma[iSpecies] = -999.;
 150                 fTPCsignalMinusExpected[iSpecies] = -999.;
 151                 fTPCNsigma[iSpecies] = -999.;
 152                 fY[iSpecies] = -999.;
 153         }
 154
 155         for (Int_t iITSlayer = 0; iITSlayer < 6; iITSlayer++) {
 156                 fITSHits[iITSlayer] = kFALSE;
 157         }
 158
 159         for (Int_t iN = 0; iN < 3; ++iN) {
 160                 fTOFLabel[iN] = -1;     // Same convention as in ESDs
 161         }
 162
 163         if (track) {
 164                 fAODTrack = track;
 165                 fAODEvent = const_cast<AliAODEvent*>(track->GetAODEvent());
 166         }
 167         if (globaltrack) fAODGlobalTrack = globaltrack;
 168         if (mcparticle) fAODMCParticle = mcparticle;
 169         if (pidresponse) fPIDResponse = pidresponse;
 170
 171         // Copy AOD Track info.
 172         if (fAODTrack) {
 173                 CopyBasicTrackInfo();
 174         } else {
 175                 AliError("No Track Supplied.");
 176         }
 177
 178         // Copy the rest of the track parameters if the filtermap is nonzero.
 179         // If fFiltermap == 0, then propagation to the DCA will result in a floating point error.
 180         if (fFilterMap) {
 181
 182                 // Find the Global Track.
 183                 if (fID >= 0) fAODGlobalTrack = fAODTrack;
 184
 185                 // Copy DCA and PID info.
 186                 if (fAODGlobalTrack) {
 187                         CopyFlags();
 188                         if (fAODEvent) CopyDCAInfo();
 189                         else AliError("Couln't find AOD Event.");
 190                         CopyITSInfo();
 191                         if (fPIDResponse) CopyTPCInfo();
 192                         CopyTOFInfo();
 193                 } else {
 194                         AliError("Couldn't find Global Track.");
 195                 }
 196
 197                 // Copy MC info.
 198                 if (fAODMCParticle) {
 199                         CopyMCInfo();
 200                 }
 201
 202         }
 203
 204         // Test
 205         /*      Double_t sigmaTOFProton = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(fAODTrack, AliPID::kProton));
 206                 if ( sigmaTOFProton < 1.0) {cout<<"tofsigmabelowone: "<<sigmaTOFProton<<endl;}
 207
 208                 Double_t sigmaTPCProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(fAODTrack, AliPID::kProton));
 209                 if ( sigmaTPCProton < 1.0) {cout<<"tpcsigmabelowone: "<<sigmaTPCProton<<endl;}*/
 210 }
 211
 212 // -----------------------------------------------------------------------
 213 Bool_t AliTrackDiHadronPID::CopyBasicTrackInfo() {
 214
 215         //
 216         // Copies everything available in every AOD track.
 217         //
 218
 219         if (fDebug > 2) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 220
 221         fPt = fAODTrack->Pt();
 222         fEta = fAODTrack->Eta();
 223         fPhi = fAODTrack->Phi();
 224
 225         fY[0] = fAODTrack->Y(AliAODTrack::kPion);
 226         fY[1] = fAODTrack->Y(AliAODTrack::kKaon);
 227         fY[2] = fAODTrack->Y(AliAODTrack::kProton);
 228
 229         //fFlags = fAODTrack->GetFlags(); // FLAGS MUST BE TAKEN FROM GLOBAL TRACKS.
 230         fFilterMap = fAODTrack->GetFilterMap();
 231
 232         fID = fAODTrack->GetID();
 233         fLabel = fAODTrack->GetLabel();
 234
 235         fCharge = fAODTrack->Charge();
 236         fNclsTPC = fAODTrack->GetTPCNcls();
 237
 238         fBasicInfoAvailable = kTRUE;
 239         return fBasicInfoAvailable;
 240
 241 }
 242
 243 // -----------------------------------------------------------------------
 244 Bool_t AliTrackDiHadronPID::CopyFlags() {
 245
 246         //
 247         // Copies Flags (properly stored in global track)
 248         //
 249
 250         if (fDebug > 2) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 251
 252         // Copy Flags
 253         fFlags = fAODGlobalTrack->GetFlags();
 254 /*
 255         // Is TOF mismatch?
 256         if (AliAODTrack::kTOFmismatch&fFlags) {
 257                 fTOFMatchingStatus = kTRUE;
 258                 //cout<<"Found TOF mismatch!"<<endl;
 259         }
 260         else fTOFMatchingStatus = kFALSE;
 261 */
 262         fFlagsAvailable = kTRUE;
 263         return fFlagsAvailable;
 264
 265 }
 266
 267 // -----------------------------------------------------------------------
 268 Bool_t AliTrackDiHadronPID::CopyDCAInfo() {
 269
 270         //
 271         // Copies DCA info. (only stored in a global track)
 272         //
 273
 274         if (fDebug > 2) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 275
 276         if (fAODGlobalTrack->IsMuonTrack()) return kFALSE;
 277
 278         // Propagate track to DCA.
 279         Double_t PosAtDCA[2] = {-999,-999};
 280     Double_t covar[3] = {-999,-999,-999};
 281     //cout<<fAODTrack<<" "<<fAODGlobalTrack<<endl;
 282     Bool_t propagate = fAODGlobalTrack->PropagateToDCA(fAODEvent->GetPrimaryVertex(),fAODEvent->GetMagneticField(),100.,PosAtDCA,covar);
 283
 284     if (propagate) {
 285         fDCAxy = PosAtDCA[0];
 286         fDCAz = PosAtDCA[1];
 287     } else {
 288         //AliError("Could not propagate track to DCA.");
 289     }
 290
 291     if (propagate) fDCAInfoAvailable = kTRUE;
 292     return fDCAInfoAvailable;
 293
 294 }
 295
 296 // -----------------------------------------------------------------------
 297 Bool_t AliTrackDiHadronPID::CopyITSInfo() {
 298
 299         //
 300         // Copies ITS info.
 301         //
 302
 303         if (fDebug > 2) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 304
 305     // Get the ITS clustermap
 306     fITSClusterMap = fAODGlobalTrack->GetITSClusterMap();
 307
 308     // Copy the ITS hits.
 309     for (Int_t iITSlayer = 0; iITSlayer < 6; iITSlayer++) {
 310                 fITSHits[iITSlayer] = fAODGlobalTrack->HasPointOnITSLayer(iITSlayer);
 311         }
 312
 313     fITSInfoAvailable = kTRUE;
 314     return fITSInfoAvailable;
 315
 316 }
 317
 318 // -----------------------------------------------------------------------
 319 Bool_t AliTrackDiHadronPID::CopyTPCInfo() {
 320
 321         //
 322         // Copies TPC info. (needs global track and pid response)
 323         //
 324
 325         if (fDebug > 2) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 326
 327     // Get TPC signal.
 328     fTPCsignal = fAODGlobalTrack->GetTPCsignal();
 329
 330     // Compute expected TPC signal under pi/K/p mass assumption.
 331     AliTPCPIDResponse& TPCPIDResponse = fPIDResponse->GetTPCResponse();
 332     fTPCmomentum = fAODGlobalTrack->GetTPCmomentum();
 333
 334         fTPCsignalMinusExpected[0] = fTPCsignal - TPCPIDResponse.GetExpectedSignal(fTPCmomentum,AliPID::kPion);
 335         fTPCsignalMinusExpected[1] = fTPCsignal - TPCPIDResponse.GetExpectedSignal(fTPCmomentum,AliPID::kKaon);
 336         fTPCsignalMinusExpected[2] = fTPCsignal - TPCPIDResponse.GetExpectedSignal(fTPCmomentum,AliPID::kProton);
 337
 338         fTPCNsigma[0] = fPIDResponse->NumberOfSigmasTPC(fAODGlobalTrack, AliPID::kPion);
 339         fTPCNsigma[1] = fPIDResponse->NumberOfSigmasTPC(fAODGlobalTrack, AliPID::kKaon);
 340         fTPCNsigma[2] = fPIDResponse->NumberOfSigmasTPC(fAODGlobalTrack, AliPID::kProton);
 341
 342     fTPCInfoAvailable = kTRUE;
 343     return fTPCInfoAvailable;
 344
 345 }
 346
 347 // -----------------------------------------------------------------------
 348 Bool_t AliTrackDiHadronPID::CopyTOFInfo() {
 349
 350         //
 351         // Copies TOF info. (needs global track)
 352         //
 353
 354         if (fDebug > 2) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 355
 356     // Get TOF signal.
 357     fTOFsignal = fAODGlobalTrack->GetTOFsignal();
 358
 359     // Compute expected TOF signal under pi/K/p mass assumption.
 360     Double_t times[AliPID::kSPECIES];
 361     fAODGlobalTrack->GetIntegratedTimes(times);
 362     fTOFsignalMinusExpected[0] = fTOFsignal - times[AliPID::kPion];
 363         fTOFsignalMinusExpected[1] = fTOFsignal - times[AliPID::kKaon];
 364         fTOFsignalMinusExpected[2] = fTOFsignal - times[AliPID::kProton];
 365
 366         fTOFNsigma[0] = fPIDResponse->NumberOfSigmasTOF(fAODGlobalTrack, AliPID::kPion);
 367         fTOFNsigma[1] = fPIDResponse->NumberOfSigmasTOF(fAODGlobalTrack, AliPID::kKaon);
 368         fTOFNsigma[2] = fPIDResponse->NumberOfSigmasTOF(fAODGlobalTrack, AliPID::kProton);
 369
 370         fTOFInfoAvailable = kTRUE;
 371
 372         // Q: what do the different TOF labels mean?
 373         // It seems that in AOD090 the TOF labels aren't copied properly.
 374         //Int_t TOFlabeltmp[3] = {0};
 375         fAODGlobalTrack->GetTOFLabel(fTOFLabel);
 376         //for (Int_t iN = 0; iN < 3; ++iN) {fTOFLabel[iN] = TOFlabeltmp[iN];}
 377         /*
 378         if (fTOFLabel[1] == fLabel || fTOFLabel[2] == fLabel) {
 379                 cout<<"fLabel = " << fLabel << " fTOFLabel =  {" << fTOFLabel[0] << "," << fTOFLabel[1] << "," << fTOFLabel[2] <<"}"<<endl;
 380         }
 381         */
 382         // The following will only work in an AOD production with the fTOFlabels set.
 383         // If it wasn't set, then every track will be labeled as no match.
 384         if (fTOFLabel[0] == -1) {fTOFMatchingStatus = 2;}                       // TPC Track was not matched to any TOF hit.
 385         else if (fLabel == fTOFLabel[0]) {fTOFMatchingStatus = 0;}      // TPC Track was correctly matched to a TOF hit.
 386         else {fTOFMatchingStatus = 1;}                                                          // TPC Track was mismatched.
 387
 388         return fTOFInfoAvailable;
 389
 390 }
 391
 392 // -----------------------------------------------------------------------
 393 Bool_t AliTrackDiHadronPID::CopyMCInfo() {
 394
 395         // Copies MC info (needs an MC track with the same label)
 396
 397         // Check if the label of the current track matches the label of the
 398         // generated particle. Note that the label of the AOD track can be
 399         // negative. This means that the quality of this track is not awesome,
 400         // but that it does correspond to the MC particle.
 401
 402         if (fDebug > 2) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 403
 404         /*
 405         if (fAODMCParticle->Label() != TMath::Abs(fAODTrack->GetLabel())) {
 406                 cout<<"Label of supplied MC particle and reconstructed track do not match."<<endl;
 407                 return kFALSE;
 408         }
 409         */
 410         // Note: It seems like the Label of the AOD track points to the INDEX of the
 411         // MCPArticle, not to the label (See AliAnalysisTaskCompareAODTrackCuts.cxx)
 412
 413         fMCPt = fAODMCParticle->Pt();
 414         fMCEta = fAODMCParticle->Eta();
 415         fMCPhi = fAODMCParticle->Phi();
 416         fMCY = fAODMCParticle->Y();
 417         fPdgCode = fAODMCParticle->PdgCode();
 418
 419         TClonesArray* mcArray = 0x0;
 420         mcArray = dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 421         if (!mcArray) {
 422         AliFatal("No MC array found in the AOD.");
 423     }
 424
 425         // Primary particle
 426         if ( fAODMCParticle->IsPhysicalPrimary() ){
 427                 fIsPhysicalPrimary = kTRUE;
 428         } else {
 429                 // Safety check for mother existence.
 430                 if (fAODMCParticle->GetMother() >= 0){
 431
 432                         Int_t mcMotherPDG = -999;
 433                         Int_t firstInt = -999;
 434
 435                         AliAODMCParticle* mcMother = (AliAODMCParticle*) mcArray->At(TMath::Abs(fAODMCParticle->GetMother()));
 436                         mcMotherPDG = TMath::Abs(mcMother->GetPdgCode());
 437
 438                         // Need a way to get the first intiger, for now Marek's method:
 439                         firstInt = Int_t (mcMotherPDG/ TMath::Power(10, Int_t(TMath::Log10(mcMotherPDG))));
 440                         // cout<<"Mother PDG: "<<mcMotherPDG<<"; Firt integer: "<<firstInt<<endl;
 441
 442                         // Weak decay
 443                         if( firstInt == 3){
 444                                 fIsSecondaryFromWeakDecay = kTRUE;
 445                         // Material decay
 446                         } else {
 447                                 fIsSecondaryFromMaterial = kTRUE;
 448                         }
 449                 }
 450         }
 451
 452         fMCInfoAvailable = kTRUE;
 453         return fMCInfoAvailable;
 454
 455 }