PWG2/RESONANCES/AliRsnAnalysisPhi900GeV.cxx

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