PWG2/RESONANCES/AliRsnAnalysisPhi900GeV.cxx

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