PWG2/KINK/AliResonanceKinkLikeSign.cxx

   1 /**************************************************************************
   2  * Author: Paraskevi Ganoti, University of Athens (pganoti@phys.uoa.gr)   *
   3  * Contributors are mentioned in the code where appropriate.              *
   4  *                                                                        *
   5  * Permission to use, copy, modify and distribute this software and its   *
   6  * documentation strictly for non-commercial purposes is hereby granted   *
   7  * without fee, provided that the above copyright notice appears in all   *
   8  * copies and that both the copyright notice and this permission notice   *
   9  * appear in the supporting documentation. The authors make no claims     *
  10  * about the suitability of this software for any purpose. It is          *
  11  * provided "as is" without express or implied warranty.                  *
  12  **************************************************************************/
  13
  14 //----------------------------------------------------------------------------------------------------------------
  15 //                        class AliResonanceKinkLikeSign
  16 //        Example of an analysis task for producing a like-sign background for resonances having at least one
  17 //        kaon-kink in their decay products.
  18 //        Background is calculated from a negative kaon kink and a negative track.
  19 //-----------------------------------------------------------------------------------------------------------------
  20
  21 #include "TChain.h"
  22 #include "TTree.h"
  23 #include "TVector3.h"
  24 #include "TF1.h"
  25 #include "TH1D.h"
  26 #include <TDatabasePDG.h>
  27 #include <TParticlePDG.h>
  28 #include "TLorentzVector.h"
  29 #include "AliAnalysisTaskSE.h"
  30 #include "AliAnalysisManager.h"
  31
  32 #include "AliESDInputHandler.h"
  33
  34 #include "AliResonanceKinkLikeSign.h"
  35 #include "AliESDkink.h"
  36
  37 ClassImp(AliResonanceKinkLikeSign)
  38
  39 //________________________________________________________________________
  40 AliResonanceKinkLikeSign::AliResonanceKinkLikeSign()
  41   : AliAnalysisTaskSE(), fESD(0), fListOfHistos(0), f1(0), f2(0), fNegKaonLikeSign(0)
  42 {
  43   // Constructor
  44 }
  45
  46 //________________________________________________________________________
  47 AliResonanceKinkLikeSign::AliResonanceKinkLikeSign(const char *name)
  48   : AliAnalysisTaskSE(name), fESD(0), fListOfHistos(0), f1(0), f2(0), fNegKaonLikeSign(0)
  49
  50 {
  51   // Constructor
  52
  53   // Define input and output slots here
  54   // Input slot #0 works with a TChain
  55   DefineInput(0, TChain::Class());
  56   DefineOutput(1, TList::Class());
  57 }
  58
  59 //________________________________________________________________________
  60 void AliResonanceKinkLikeSign::ConnectInputData(Option_t *)
  61 {
  62   // Connect ESD or AOD here
  63   // Called once
  64
  65   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
  66   if (!tree) {
  67     Printf("ERROR: Could not read chain from input slot 0");
  68   } else {
  69
  70     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
  71
  72     if (!esdH) {
  73       Printf("ERROR: Could not get ESDInputHandler");
  74     } else
  75       fESD = esdH->GetEvent();
  76   }
  77 }
  78
  79 //________________________________________________________________________
  80 void AliResonanceKinkLikeSign::CreateOutputObjects()
  81 {
  82   // Create histograms
  83   // Called once
  84
  85    f1=new TF1("f1","((atan([0]*[1]*(1.0/(sqrt((x^2)*(1.0-([1]^2))-([0]^2)*([1]^2))))))*180.)/[2]",1.1,10.0);
  86    f1->SetParameter(0,0.493677);
  87    f1->SetParameter(1,0.9127037);
  88    f1->SetParameter(2,TMath::Pi());
  89
  90    f2=new TF1("f2","((atan([0]*[1]*(1.0/(sqrt((x^2)*(1.0-([1]^2))-([0]^2)*([1]^2))))))*180.)/[2]",0.1,10.0);
  91    f2->SetParameter(0,0.13957018);
  92    f2->SetParameter(1,0.2731374);
  93    f2->SetParameter(2,TMath::Pi());
  94
  95    //OpenFile(1);  //uncomment for proof analysis
  96
  97     // for K*0(892)
  98    fNegKaonLikeSign=new TH1D("fNegKaonLikeSign"," ", 60,0.6,1.2);
  99
 100    // for phi(1020)
 101  //  fNegKaonLikeSign=new TH1D("fNegKaonLikeSign"," ", 70,0.99,1.088);
 102
 103
 104    fListOfHistos=new TList();
 105    fListOfHistos->Add(fNegKaonLikeSign);
 106
 107 }
 108
 109 //________________________________________________________________________
 110 void AliResonanceKinkLikeSign::Exec(Option_t *)
 111 {
 112   // Main loop
 113   // Called for each event
 114
 115    Double_t daughter1Mass=TDatabasePDG::Instance()->GetParticle(kKPlus)->Mass();
 116    Double_t daughter2Mass=TDatabasePDG::Instance()->GetParticle(kPiPlus)->Mass();
 117
 118   if (!fESD) {
 119     Printf("ERROR: fESD not available");
 120     return;
 121   }
 122
 123   const AliESDVertex* vertex = GetEventVertex(fESD);
 124   if(!vertex) return;
 125
 126   Double_t ptrackpos[3], ptrackneg[3];
 127
 128   TLorentzVector p4pos;
 129   TLorentzVector p4neg;
 130   TLorentzVector p4comb;
 131
 132   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 133     AliESDtrack* trackpos = fESD->GetTrack(iTracks);
 134     if (!trackpos) {
 135       Printf("ERROR: Could not receive track %d", iTracks);
 136       continue;
 137     }
 138
 139     if (trackpos->GetSign() > 0) continue;
 140
 141     trackpos->GetPxPyPz(ptrackpos);
 142
 143     Float_t nSigmaToVertex = GetSigmaToVertex(trackpos);
 144
 145     Float_t bpos[2];
 146     Float_t bCovpos[3];
 147     trackpos->GetImpactParameters(bpos,bCovpos);
 148
 149     if (bCovpos[0]<=0 || bCovpos[2]<=0) {
 150      Printf("Estimated b resolution lower or equal zero!");
 151      bCovpos[0]=0; bCovpos[2]=0;
 152     }
 153
 154     Float_t dcaToVertexXYpos = bpos[0];
 155     Float_t dcaToVertexZpos = bpos[1];
 156
 157     if(nSigmaToVertex>=4) continue;
 158     if((dcaToVertexXYpos>3.0)||(dcaToVertexZpos>3.0)) continue;
 159
 160     TVector3 posTrackMom(ptrackpos[0],ptrackpos[1],ptrackpos[2]);
 161
 162     if(posTrackMom.Perp()<=0.25) continue;
 163
 164     //Uncomment the following block if the Like Sign is made of K- kink + negative track
 165
 166     //Int_t indexKink=trackpos->GetKinkIndex(0);
 167     //Int_t kaonKinkFlag=0;
 168 //     if(indexKink<0){
 169 //
 170 //         AliESDkink *kink=fESD->GetKink(TMath::Abs(IndexKink)-1);
 171 //      const TVector3 motherMfromKink(kink->GetMotherP());
 172 //      const TVector3 daughterMKink(kink->GetDaughterP());
 173 //      Float_t qt=kink->GetQt();
 174 //
 175 //         Double_t maxDecAngKmu=f1->Eval(motherMfromKink.Mag(),0.,0.,0.);
 176 //         Double_t maxDecAngpimu=f2->Eval(motherMfromKink.Mag(),0.,0.,0.);
 177 //
 178 //         Float_t kinkAngle=TMath::RadToDeg()*kink->GetAngle(2);
 179 //
 180 //      Float_t energyDaughterMu=TMath::Sqrt(daughterMKink.Mag()*daughterMKink.Mag()+0.105658*0.105658);
 181 //         Float_t p1XM= motherMfromKink.Px();
 182 //         Float_t p1YM= motherMfromKink.Py();
 183 //         Float_t p1ZM= motherMfromKink.Pz();
 184 //         Float_t p2XM= daughterMKink.Px();
 185 //         Float_t p2YM= daughterMKink.Py();
 186 //         Float_t p2ZM= daughterMKink.Pz();
 187 //         Float_t p3Daughter=TMath::Sqrt(((p1XM-p2XM)*(p1XM-p2XM))+((p1YM-p2YM)*(p1YM-p2YM))+((p1ZM-p2ZM)*(p1ZM-p2ZM)));
 188 //         Double_t invariantMassKmu= TMath::Sqrt((energyDaughterMu+p3Daughter)*(energyDaughterMu+p3Daughter)-motherMfromKink.Mag()*motherMfromKink.Mag());
 189 //
 190 //         if((kinkAngle>maxDecAngpimu)&&(qt>0.05)&&(qt<0.25)&&((kink->GetR()>110.)&&(kink->GetR()<230.))&&(TMath::Abs(posTrackMom.Eta())<1.1)&&(invariantMassKmu<0.6)) {
 191 //
 192 //           if(posTrackMom.Mag()<=1.1) {
 193 //         kaonKinkFlag=1;
 194 //        }
 195 //        else
 196 //        if (kinkAngle<maxDecAngKmu) {
 197 //         kaonKinkFlag=1;
 198 //        }
 199 //      }
 200 //
 201 //     }  //End Kink Information
 202 //
 203 //     if(kaonKinkFlag==0) continue;
 204 //     if(kaonKinkFlag==1) p4pos.SetVectM(posTrackMom,daughter1Mass);
 205
 206       // Comment the following statements till the "for" if the Like Sign of K- kink + negative track is needed
 207
 208       UInt_t status=trackpos->GetStatus();
 209       if((status&AliESDtrack::kTPCrefit)==0) continue;
 210       if(trackpos->GetTPCclusters(0)<50) continue;
 211       if((trackpos->GetTPCchi2()/trackpos->GetTPCclusters(0))>3.5) continue;
 212       Double_t extCovPos[15];
 213       trackpos->GetExternalCovariance(extCovPos);
 214       if(extCovPos[0]>2) continue;
 215       if(extCovPos[2]>2) continue;
 216       if(extCovPos[5]>0.5) continue;
 217       if(extCovPos[9]>0.5) continue;
 218       if(extCovPos[14]>2) continue;
 219
 220        p4pos.SetVectM(posTrackMom,daughter1Mass);
 221
 222       for (Int_t j=0; j<fESD->GetNumberOfTracks(); j++) {
 223         if(iTracks==j) continue;
 224         AliESDtrack* trackneg=fESD->GetTrack(j);
 225         if (trackneg->GetSign() > 0) continue;
 226
 227         trackneg->GetPxPyPz(ptrackneg);
 228         Float_t negSigmaToVertex = GetSigmaToVertex(trackneg);
 229
 230         Float_t bneg[2];
 231         Float_t bCovneg[3];
 232         trackneg->GetImpactParameters(bneg,bCovneg);
 233         if (bCovneg[0]<=0 || bCovneg[2]<=0) {
 234           Printf("Estimated b resolution lower or equal zero!");
 235           bCovneg[0]=0; bCovneg[2]=0;
 236         }
 237
 238         Float_t dcaToVertexXYneg = bneg[0];
 239         Float_t dcaToVertexZneg = bneg[1];
 240
 241         if(negSigmaToVertex>=4) continue;
 242         if((dcaToVertexXYneg>3.0)||(dcaToVertexZneg>3.0)) continue;
 243
 244         TVector3 negTrackMom(ptrackneg[0],ptrackneg[1],ptrackneg[2]);
 245
 246         if(negTrackMom.Perp()<=0.25) continue;
 247
 248         UInt_t statusneg=trackneg->GetStatus();
 249
 250         if((statusneg&AliESDtrack::kTPCrefit)==0) continue;
 251
 252         if(trackneg->GetTPCclusters(0)<50) continue;
 253         if((trackneg->GetTPCchi2()/trackneg->GetTPCclusters(0))>3.5) continue;
 254         Double_t extCovneg[15];
 255         trackneg->GetExternalCovariance(extCovneg);
 256         if(extCovneg[0]>2) continue;
 257         if(extCovneg[2]>2) continue;
 258         if(extCovneg[5]>0.5) continue;
 259         if(extCovneg[9]>0.5) continue;
 260         if(extCovneg[14]>2) continue;
 261
 262         p4neg.SetVectM(negTrackMom, daughter2Mass);
 263
 264         p4comb=p4pos;
 265         p4comb+=p4neg;
 266         fNegKaonLikeSign->Fill(p4comb.M());
 267
 268       } //inner track loop
 269
 270   } //outer track loop
 271
 272     PostData(1, fListOfHistos);
 273 }
 274
 275 //________________________________________________________________________
 276 void AliResonanceKinkLikeSign::Terminate(Option_t *)
 277 {
 278
 279 }
 280
 281 //____________________________________________________________________//
 282
 283 Float_t AliResonanceKinkLikeSign::GetSigmaToVertex(AliESDtrack* esdTrack) const {
 284   // Calculates the number of sigma to the vertex.
 285
 286   Float_t b[2];
 287   Float_t bRes[2];
 288   Float_t bCov[3];
 289
 290     esdTrack->GetImpactParameters(b,bCov);
 291
 292   if (bCov[0]<=0 || bCov[2]<=0) {
 293     //AliDebug(1, "Estimated b resolution lower or equal zero!");
 294     bCov[0]=0; bCov[2]=0;
 295   }
 296   bRes[0] = TMath::Sqrt(bCov[0]);
 297   bRes[1] = TMath::Sqrt(bCov[2]);
 298
 299   if (bRes[0] == 0 || bRes[1] ==0) return -1;
 300
 301   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 302
 303   if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
 304
 305   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 306
 307   return d;
 308 }
 309
 310 //____________________________________________________________________//
 311
 312 const AliESDVertex* AliResonanceKinkLikeSign::GetEventVertex(const AliESDEvent* esd) const
 313
 314 {
 315   // Get the vertex
 316
 317   const AliESDVertex* vertex = esd->GetPrimaryVertex();
 318
 319   if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
 320   else
 321   {
 322      vertex = esd->GetPrimaryVertexSPD();
 323      if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
 324      else
 325      return 0;
 326   }
 327 }
 328