PWG2/RESONANCES/AliRsnAnalysisPhi900GeV.cxx

   1 #include "TTree.h"
   2 #include "TParticle.h"
   3 #include "TRandom.h"
   4 #include "TLorentzVector.h"
   5
   6 #include "AliLog.h"
   7 #include "AliESDEvent.h"
   8 #include "AliESDVertex.h"
   9 #include "AliESDtrack.h"
  10 #include "AliStack.h"
  11 #include "AliMCEvent.h"
  12
  13 #include "AliRsnAnalysisPhi900GeV.h"
  14
  15
  16 AliRsnAnalysisPhi900GeV::AliRsnAnalysisPhi900GeV(const char *name) :
  17   AliAnalysisTaskSE(name),
  18   fUseMC(kFALSE),
  19   fPDG(0),
  20   fIM(0.0),
  21   fPt(0.0),
  22   fY(0.0),
  23   fEta(0.0),
  24   fDCAr(1E6),
  25   fDCAz(1E6),
  26   fChi2(1E6),
  27   fNTPC(0),
  28   fMinTPC(-1E6),
  29   fMaxTPC( 1E6),
  30   fOutTree(0x0),
  31   fTOFESD(kFALSE),
  32   fTOFSigma(210.0),
  33   fTOFmaker(0x0),
  34   fTOFSettings(AliRsnTOFT0maker::kNone)
  35 {
  36 //
  37 // Constructor
  38 //
  39
  40   DefineOutput(1, TTree::Class());
  41 }
  42
  43
  44 AliRsnAnalysisPhi900GeV::AliRsnAnalysisPhi900GeV(const AliRsnAnalysisPhi900GeV& copy) :
  45   AliAnalysisTaskSE(copy),
  46   fUseMC(copy.fUseMC),
  47   fPDG(0),
  48   fIM(0.0),
  49   fPt(0.0),
  50   fY(0.0),
  51   fEta(0.0),
  52   fDCAr(copy.fDCAr),
  53   fDCAz(copy.fDCAz),
  54   fChi2(copy.fChi2),
  55   fNTPC(copy.fNTPC),
  56   fMinTPC(copy.fMinTPC),
  57   fMaxTPC(copy.fMaxTPC),
  58   fOutTree(0x0),
  59   fTOFESD(copy.fTOFESD),
  60   fTOFSigma(copy.fTOFSigma),
  61   fTOFmaker(0x0),
  62   fTOFSettings(copy.fTOFSettings)
  63 {
  64 //
  65 // Copy constructor
  66 //
  67 }
  68
  69
  70 AliRsnAnalysisPhi900GeV& AliRsnAnalysisPhi900GeV::operator=(const AliRsnAnalysisPhi900GeV& copy)
  71 {
  72 //
  73 // Assignment operator
  74 //
  75
  76   fUseMC = copy.fUseMC;
  77
  78   fDCAr = copy.fDCAr;
  79   fDCAz = copy.fDCAz;
  80   fChi2 = copy.fChi2;
  81   fNTPC = copy.fNTPC;
  82
  83   fMinTPC = copy.fMinTPC;
  84   fMaxTPC = copy.fMaxTPC;
  85
  86   fTOFESD      = copy.fTOFESD;
  87   fTOFSigma    = copy.fTOFSigma;
  88   fTOFSettings = copy.fTOFSettings;
  89
  90   return (*this);
  91 }
  92
  93
  94 AliRsnAnalysisPhi900GeV::~AliRsnAnalysisPhi900GeV()
  95 {
  96 //
  97 // Destructor
  98 //
  99
 100   if (fOutTree) delete fOutTree;
 101 }
 102
 103
 104 void AliRsnAnalysisPhi900GeV::UserCreateOutputObjects()
 105 {
 106 //
 107 // Create the output data container
 108 //
 109
 110   // setup TOF maker
 111   fTOFmaker = new AliRsnTOFT0maker;
 112   fTOFmaker->SetTimeResolution(fTOFSigma * 1E-12);
 113   fTOFmaker->SetESDdata(fTOFESD);
 114   fTOFmaker->fSettings = fTOFSettings;
 115   AliInfo(Form("TOF sigma    = %f", fTOFSigma));
 116   AliInfo(Form("TOF ESD      = %s", (fTOFESD ? "YES" : "NO")));
 117   AliInfo(Form("TOF settings = %s", fTOFmaker->Settings().Data()));
 118
 119   // load dead channel map
 120   fTOFmaker->LoadChannelMap("tofmap.root");
 121   fTOFmaker->SetMaskOffChannel();
 122
 123   // initialize random
 124   gRandom->SetSeed(0);
 125
 126   // create output trees
 127   OpenFile(1);
 128   fOutTree = new TTree("rsnTree", "Pairs");
 129
 130   fOutTree->Branch("pdg", &fPDG, "pdg/S");
 131   fOutTree->Branch("im" , &fIM , "im/F" );
 132   fOutTree->Branch("y"  , &fY  , "y/F"  );
 133   fOutTree->Branch("pt" , &fPt , "pt/F" );
 134   fOutTree->Branch("eta", &fEta, "eta/F");
 135 }
 136
 137
 138 void AliRsnAnalysisPhi900GeV::UserExec(Option_t *)
 139 {
 140 //
 141 // Main execution function
 142 //
 143
 144   // retrieve ESD event and related stack (if available)
 145   AliESDEvent *esd   = dynamic_cast<AliESDEvent*>(fInputEvent);
 146   AliStack    *stack = (fMCEvent ? fMCEvent->Stack() : 0x0);
 147
 148   // get the best primary vertex:
 149   // first try the one with tracks
 150   Int_t type = 0;
 151   const AliESDVertex *v = esd->GetPrimaryVertexTracks();
 152   if(v->GetNContributors() < 1)
 153   {
 154     // if not good, try SPD vertex
 155     type = 1;
 156     v = esd->GetPrimaryVertexSPD();
 157
 158     // if this is not good skip this event
 159     if (v->GetNContributors() < 1) return;
 160   }
 161
 162   // if the Z position is larger than 10, skip this event
 163   if (TMath::Abs(v->GetZv()) > 10.0) return;
 164
 165   // smear TOF times in case of MC
 166   if (stack) RemakeTOFtimeMC(esd);
 167
 168   // get time zero for TOF
 169   Double_t *tof = fTOFmaker->RemakePID(esd);
 170
 171   if (!fUseMC) ProcessESD(esd, v, tof[0], stack);
 172   else ProcessMC(stack);
 173 }
 174
 175
 176 void AliRsnAnalysisPhi900GeV::Terminate(Option_t *)
 177 {
 178 //
 179 // Terminate
 180 //
 181 }
 182
 183
 184 void AliRsnAnalysisPhi900GeV::ProcessESD
 185 (AliESDEvent *esd, const AliESDVertex *v, Double_t time0, AliStack *stack)
 186 {
 187 //
 188 // This function works with the ESD object
 189 //
 190
 191   // prepare to look on all tracks to select the ones
 192   // which pass all the cuts
 193   Int_t   ntracks = esd->GetNumberOfTracks();
 194   TArrayI pos(ntracks);
 195   TArrayI neg(ntracks);
 196
 197   // define fixed functions for TOF compatibility range
 198   Double_t a1 = 0.01, a2 = -0.03;
 199   Double_t b1 = 0.25, b2 =  0.25;
 200   Double_t c1 = 0.05, c2 = -0.03;
 201   Double_t ymax, ymin;
 202
 203   // loop on all tracks
 204   Int_t    i, charge, nSPD, npos = 0, nneg = 0;
 205   Float_t  chi2, b[2], bCov[3];
 206   Double_t tpc, bb, mom, tofTime, tofRef, tofRel, times[10];
 207   Bool_t   okTOF;
 208   for (i = 0; i < ntracks; i++)
 209   {
 210     AliESDtrack *track = esd->GetTrack(i);
 211     if (!track) continue;
 212
 213     // skip if it has not the required flags
 214     if (!track->IsOn(AliESDtrack::kTPCin)) continue;
 215     if (!track->IsOn(AliESDtrack::kTPCrefit)) continue;
 216     if (!track->IsOn(AliESDtrack::kITSrefit)) continue;
 217
 218     // skip if it has not the TPC inner wall projection
 219     if (!track->GetInnerParam()) continue;
 220
 221     // skip kink daughters
 222     if ((Int_t)track->GetKinkIndex(0) > 0) continue;
 223
 224     // check clusters in TPC
 225     if (track->GetTPCclusters(0) < fNTPC) continue;
 226
 227     // check chi2
 228     chi2  = (Float_t)track->GetTPCchi2();
 229     chi2 /= (Float_t)track->GetTPCclusters(0);
 230     if (chi2 > fChi2) continue;
 231
 232     // check that has at least 1 SPD cluster
 233     nSPD = 0;
 234     if (track->HasPointOnITSLayer(0)) nSPD++;
 235     if (track->HasPointOnITSLayer(1)) nSPD++;
 236     if (nSPD < 1) continue;
 237
 238     // check primary by reverting to vertex
 239     // and checking DCA
 240     if (!track->RelateToVertex(v, esd->GetMagneticField(), kVeryBig)) continue;
 241     track->GetImpactParameters(b, bCov);
 242     if (b[0] > fDCAr) continue;
 243     if (b[1] > fDCAz) continue;
 244
 245     // check TPC dE/dx
 246     AliExternalTrackParam trackIn(*track->GetInnerParam());
 247     mom = trackIn.P();
 248     tpc = (Double_t)track->GetTPCsignal();
 249     bb  = AlephBB(mom);
 250     tpc = (tpc - bb) / bb;
 251     if (tpc < fMinTPC || tpc > fMaxTPC) continue;
 252
 253     // if possible, check TOF
 254     okTOF = kTRUE;
 255     if (track->IsOn(AliESDtrack::kTOFpid))
 256     {
 257       mom = track->P();
 258       if (mom <= 0.26)
 259         okTOF = kTRUE;
 260       else
 261       {
 262         track->GetIntegratedTimes(times);
 263         tofTime = (Double_t)track->GetTOFsignal() - time0;
 264         tofRef  = times[AliPID::kKaon];
 265         tofRel  = (tofTime - tofRef) / tofRef;
 266         ymax    = a1 / (mom - b1) + c1;
 267         ymin    = a2 / (mom - b2) + c2;
 268         okTOF   = (tofRel >= ymin && tofRel <= ymax);
 269       }
 270     }
 271     if (!okTOF) continue;
 272
 273     // if we arrive here, all cuts were passed
 274     // and we add the track to one array depending on charge
 275     charge = (Int_t)track->Charge();
 276     if (charge > 0)
 277       pos[npos++] = i;
 278     else if (charge < 0)
 279       neg[nneg++] = i;
 280   }
 281
 282   // resize arrays accordingly
 283   pos.Set(npos);
 284   neg.Set(nneg);
 285
 286   // loop to compute invariant mass
 287   Int_t           ip, in, lp, ln;
 288   AliPID          pid;
 289   Double_t        kmass = pid.ParticleMass(AliPID::kKaon);
 290   Double_t        phimass = 1.019455;
 291   TParticle      *partp = 0x0, *partn = 0x0;
 292   AliESDtrack    *tp = 0x0, *tn = 0x0;
 293   TLorentzVector  vp, vn, vsum, vref;
 294   for (ip = 0; ip < npos; ip++)
 295   {
 296     tp = esd->GetTrack(pos[ip]);
 297     lp = TMath::Abs(tp->GetLabel());
 298     if (stack) partp = stack->Particle(lp);
 299
 300     for (in = 0; in < nneg; in++)
 301     {
 302       if (pos[ip] == neg[in]) continue;
 303       tn = esd->GetTrack(neg[in]);
 304       ln = TMath::Abs(tn->GetLabel());
 305       if (stack) partn = stack->Particle(ln);
 306
 307       fPDG = 0;
 308       if (partp && partn)
 309       {
 310         if (partp->GetFirstMother() == partn->GetFirstMother())
 311         {
 312           if (partp->GetFirstMother() > 0)
 313           {
 314             TParticle *mum = stack->Particle(partp->GetFirstMother());
 315             fPDG = mum->GetPdgCode();
 316           }
 317         }
 318       }
 319       fPDG = TMath::Abs(fPDG);
 320
 321       vp.SetXYZM(tp->Px(), tp->Py(), tp->Pz(), kmass);
 322       vn.SetXYZM(tn->Px(), tn->Py(), tn->Pz(), kmass);
 323       vsum = vp + vn;
 324       vref.SetXYZM(vsum.X(), vsum.Y(), vsum.Z(), phimass);
 325
 326       fIM  = (Float_t)vsum.M();
 327       fPt  = (Float_t)vsum.Perp();
 328       fEta = (Float_t)vsum.Eta();
 329       fY   = (Float_t)vref.Rapidity();
 330
 331       fOutTree->Fill();
 332     }
 333   }
 334 }
 335
 336
 337 void AliRsnAnalysisPhi900GeV::ProcessMC(AliStack *stack)
 338 {
 339 //
 340 // Function to process stack only
 341 //
 342
 343   if (!stack) return;
 344   Int_t nPart = stack->GetNtrack();
 345
 346   // loop to compute invariant mass
 347   Int_t           ip, in;
 348   AliPID          pid;
 349   Double_t        kmass = pid.ParticleMass(AliPID::kKaon);
 350   Double_t        phimass = 1.019455;
 351   TParticle      *partp = 0x0, *partn = 0x0;
 352   TLorentzVector  vp, vn, vsum, vref;
 353
 354   for (ip = 0; ip < nPart; ip++)
 355   {
 356     partp = stack->Particle(ip);
 357     if (partp->GetPdgCode() != 321) continue;
 358
 359     for (in = 0; in < nPart; in++)
 360     {
 361       partn = stack->Particle(in);
 362       if (partn->GetPdgCode() != -321) continue;
 363
 364       fPDG = 0;
 365       if (partp->GetFirstMother() == partn->GetFirstMother())
 366       {
 367         if (partp->GetFirstMother() > 0)
 368         {
 369           TParticle *mum = stack->Particle(partp->GetFirstMother());
 370           fPDG = mum->GetPdgCode();
 371         }
 372       }
 373       fPDG = TMath::Abs(fPDG);
 374       if (fPDG != 333) continue;
 375
 376       vp.SetXYZM(partp->Px(), partp->Py(), partp->Pz(), kmass);
 377       vn.SetXYZM(partn->Px(), partn->Py(), partn->Pz(), kmass);
 378       vsum = vp + vn;
 379       vref.SetXYZM(vsum.X(), vsum.Y(), vsum.Z(), phimass);
 380
 381       fIM  = (Float_t)vsum.M();
 382       fPt  = (Float_t)vsum.Perp();
 383       fEta = (Float_t)vsum.Eta();
 384       fY   = (Float_t)vref.Rapidity();
 385
 386       fOutTree->Fill();
 387     }
 388   }
 389 }
 390
 391
 392 Double_t AliRsnAnalysisPhi900GeV::AlephBB(Double_t p, Double_t mass)
 393 {
 394 //
 395 // Compute expected Bethe-Bloch for that momentum and mass
 396 //
 397
 398   if (mass < 1E-6) return 0.0;
 399
 400   Double_t aa, bb, out, beta, bg;
 401
 402   bg   = p / mass;
 403   beta = bg / TMath::Sqrt(1.0 + bg * bg);
 404   aa   = TMath::Power(beta, fTPCpar[3]);
 405   bb   = TMath::Power(1./bg, fTPCpar[4]);
 406   bb   = TMath::Log(fTPCpar[2] + bb);
 407   out  = (fTPCpar[1] - aa - bb) * fTPCpar[0]/aa;
 408
 409   return out;
 410 }
 411
 412
 413 Double_t AliRsnAnalysisPhi900GeV::RemakeTOFtimeMC(AliESDEvent *& esd)
 414 {
 415 //
 416 // Smear initial time for TOF in order to reproduce data resolution
 417 //
 418
 419   Double_t t0 = gRandom->Gaus(0,135.); //spread in ps
 420   Int_t ntracks = esd->GetNumberOfTracks();
 421   Double_t sigmaStandard = 80.; // 80 ps from TOF
 422
 423   while (ntracks--)
 424   {
 425     AliESDtrack *t = esd->GetTrack(ntracks);
 426     if ((t->GetStatus()&AliESDtrack::kTOFout) == 0 || (t->GetStatus()&AliESDtrack::kTIME)==0) continue;
 427     Double_t time = t->GetTOFsignal();
 428     if(fTOFSigma > sigmaStandard)
 429     {
 430       Double_t sigmaAdded = TMath::Sqrt(fTOFSigma*fTOFSigma - sigmaStandard*sigmaStandard);
 431       Double_t timerandomtrack = gRandom->Gaus(0, sigmaAdded); //spread in ps
 432       time += timerandomtrack;
 433    }
 434    time += t0;
 435    t->SetTOFsignal(time);
 436  }
 437  return t0;
 438 }