PWG2/RESONANCES/AliRsnAnalysisPhi7TeV.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 #include <iomanip>
   9
  10 #include "TH1.h"
  11 #include "TTree.h"
  12 #include "TParticle.h"
  13 #include "TRandom.h"
  14 #include "TLorentzVector.h"
  15
  16 #include "AliLog.h"
  17 #include "AliESDpid.h"
  18 #include "AliESDEvent.h"
  19 #include "AliESDVertex.h"
  20 #include "AliESDtrack.h"
  21 #include "AliStack.h"
  22 #include "AliMCEvent.h"
  23 #include "AliTOFT0maker.h"
  24 #include "AliTOFcalib.h"
  25 #include "AliCDBManager.h"
  26 #include "AliITSPIDResponse.h"
  27
  28 #include "AliRsnAnalysisPhi7TeV.h"
  29
  30 //__________________________________________________________________________________________________
  31 AliRsnAnalysisPhi7TeV::AliRsnAnalysisPhi7TeV(const char *name) :
  32   AliAnalysisTaskSE(name),
  33   fUseMC(kFALSE),
  34   fCheckITS(kTRUE),
  35   fCheckTPC(kTRUE),
  36   fCheckTOF(kTRUE),
  37   fPDG(0),
  38   fCh(0),
  39   fIM(0.0),
  40   fPt(0.0),
  41   fY(0.0),
  42   fEta(0.0),
  43   fMaxVz(1E6),
  44   fMaxITSband(1E6),
  45   fTPCpLimit(0.35),
  46   fMinTPCband(-1E6),
  47   fMaxTPCband( 1E6),
  48   fRsnTreeComp(0x0),
  49   fRsnTreeTrue(0x0),
  50   fOutList(0x0),
  51   fHEvents(0x0),
  52   fESDtrackCutsTPC(),
  53   fESDtrackCutsITS(),
  54   fESDpid(0x0),
  55   fTOFmaker(0x0),
  56   fTOFcalib(0x0),
  57   fTOFcalibrateESD(kFALSE),
  58   fTOFcorrectTExp(kFALSE),
  59   fTOFuseT0(kFALSE),
  60   fTOFtuneMC(kFALSE),
  61   fTOFresolution(0.0)
  62
  63 {
  64 //
  65 // Constructor
  66 //
  67
  68   DefineOutput(1, TTree::Class());
  69   DefineOutput(2, TTree::Class());
  70   DefineOutput(3, TList::Class());
  71 }
  72
  73 //__________________________________________________________________________________________________
  74 AliRsnAnalysisPhi7TeV::AliRsnAnalysisPhi7TeV(const AliRsnAnalysisPhi7TeV& copy) :
  75   AliAnalysisTaskSE(copy),
  76   fUseMC(copy.fUseMC),
  77   fCheckITS(copy.fCheckITS),
  78   fCheckTPC(copy.fCheckTPC),
  79   fCheckTOF(copy.fCheckTOF),
  80   fPDG(0),
  81   fCh(0),
  82   fIM(0.0),
  83   fPt(0.0),
  84   fY(0.0),
  85   fEta(0.0),
  86   fMaxVz(copy.fMaxVz),
  87   fMaxITSband(copy.fMaxITSband),
  88   fTPCpLimit(copy.fTPCpLimit),
  89   fMinTPCband(copy.fMinTPCband),
  90   fMaxTPCband(copy.fMaxTPCband),
  91   fRsnTreeComp(0x0),
  92   fRsnTreeTrue(0x0),
  93   fOutList(0x0),
  94   fHEvents(0x0),
  95   fESDtrackCutsTPC(copy.fESDtrackCutsTPC),
  96   fESDtrackCutsITS(copy.fESDtrackCutsITS),
  97   fESDpid(0x0),
  98   fTOFmaker(0x0),
  99   fTOFcalib(0x0),
 100   fTOFcalibrateESD(kFALSE),
 101   fTOFcorrectTExp(kFALSE),
 102   fTOFuseT0(kFALSE),
 103   fTOFtuneMC(kFALSE),
 104   fTOFresolution(0.0)
 105 {
 106 //
 107 // Copy constructor
 108 //
 109 }
 110
 111 //__________________________________________________________________________________________________
 112 AliRsnAnalysisPhi7TeV& AliRsnAnalysisPhi7TeV::operator=(const AliRsnAnalysisPhi7TeV& copy)
 113 {
 114 //
 115 // Assignment operator
 116 //
 117
 118   fUseMC = copy.fUseMC;
 119   fCheckITS = copy.fCheckITS;
 120   fCheckTPC = copy.fCheckTPC;
 121   fCheckTOF = copy.fCheckTOF;
 122
 123   fMaxVz   = copy.fMaxVz;
 124   fMaxITSband = copy.fMaxITSband;
 125
 126   fTPCpLimit  = copy.fTPCpLimit;
 127   fMinTPCband = copy.fMinTPCband;
 128   fMaxTPCband = copy.fMaxTPCband;
 129
 130   fESDtrackCutsTPC = copy.fESDtrackCutsTPC;
 131   fESDtrackCutsITS = copy.fESDtrackCutsITS;
 132
 133   fTOFcalibrateESD = copy.fTOFcalibrateESD;
 134   fTOFcorrectTExp = copy.fTOFcorrectTExp;
 135   fTOFuseT0 = copy.fTOFuseT0;
 136   fTOFtuneMC = copy.fTOFtuneMC;
 137   fTOFresolution = copy.fTOFresolution;
 138
 139   return (*this);
 140 }
 141
 142 //__________________________________________________________________________________________________
 143 AliRsnAnalysisPhi7TeV::~AliRsnAnalysisPhi7TeV()
 144 {
 145 //
 146 // Destructor
 147 //
 148
 149   if (fRsnTreeComp) delete fRsnTreeComp;
 150   if (fRsnTreeTrue) delete fRsnTreeTrue;
 151   if (fHEvents)     delete fHEvents;
 152   if (fESDpid)      delete fESDpid;
 153 }
 154
 155 //__________________________________________________________________________________________________
 156 void AliRsnAnalysisPhi7TeV::UserCreateOutputObjects()
 157 {
 158 //
 159 // Create the output data container
 160 //
 161
 162   // setup TPC response
 163   fESDpid = new AliESDpid;
 164   fESDpid->GetTPCResponse().SetBetheBlochParameters(fTPCpar[0],fTPCpar[1],fTPCpar[2],fTPCpar[3],fTPCpar[4]);
 165
 166   // setup TOF maker & calibration
 167   fTOFcalib = new AliTOFcalib;
 168   fTOFmaker = new AliTOFT0maker(fESDpid, fTOFcalib);
 169   fTOFmaker->SetTimeResolution(fTOFresolution);
 170
 171   // initialize random
 172   gRandom->SetSeed(0);
 173
 174   // create output trees
 175   OpenFile(1);
 176   fRsnTreeComp = new TTree("rsnTree", "Pairs");
 177
 178   fRsnTreeComp->Branch("pdg", &fPDG, "pdg/S"   );
 179   fRsnTreeComp->Branch("ch" , &fCh , "ch/S"    );
 180   fRsnTreeComp->Branch("im" , &fIM , "im/F"    );
 181   fRsnTreeComp->Branch("y"  , &fY  , "y/F"     );
 182   fRsnTreeComp->Branch("pt" , &fPt , "pt/F"    );
 183   fRsnTreeComp->Branch("eta", &fEta, "eta/F"   );
 184   fRsnTreeComp->Branch("its", &fITS, "its[2]/S");
 185
 186   OpenFile(2);
 187   fRsnTreeTrue = new TTree("rsnTrue", "True pairs");
 188
 189   fRsnTreeTrue->Branch("im" , &fIM , "im/F" );
 190   fRsnTreeTrue->Branch("y"  , &fY  , "y/F"  );
 191   fRsnTreeTrue->Branch("pt" , &fPt , "pt/F" );
 192   fRsnTreeTrue->Branch("eta", &fEta, "eta/F");
 193
 194   OpenFile(3);
 195   fOutList    = new TList;
 196   fHEvents    = new TH1I("hEvents", "Event details", 5, 0, 5);
 197   fVertexX[0] = new TH1F("hVertexTracksX", "X position of primary vertex (tracks)", 200,  -2,  2);
 198   fVertexY[0] = new TH1F("hVertexTracksY", "Y position of primary vertex (tracks)", 200,  -2,  2);
 199   fVertexZ[0] = new TH1F("hVertexTracksZ", "Z position of primary vertex (tracks)", 400, -40, 40);
 200   fVertexX[1] = new TH1F("hVertexSPDX", "X position of primary vertex (SPD)", 200,  -2,  2);
 201   fVertexY[1] = new TH1F("hVertexSPDY", "Y position of primary vertex (SPD)", 200,  -2,  2);
 202   fVertexZ[1] = new TH1F("hVertexSPDZ", "Z position of primary vertex (SPD)", 400, -40, 40);
 203
 204   fHEvents->GetXaxis()->SetBinLabel(1, "Good vertex with tracks");
 205   fHEvents->GetXaxis()->SetBinLabel(2, "Good vertex with SPD");
 206   fHEvents->GetXaxis()->SetBinLabel(3, "Far vertex with tracks");
 207   fHEvents->GetXaxis()->SetBinLabel(4, "Far vertex with SPD");
 208   fHEvents->GetXaxis()->SetBinLabel(5, "No good vertex");
 209
 210   fOutList->Add(fHEvents);
 211   fOutList->Add(fVertexX[0]);
 212   fOutList->Add(fVertexY[0]);
 213   fOutList->Add(fVertexZ[0]);
 214   fOutList->Add(fVertexX[1]);
 215   fOutList->Add(fVertexY[1]);
 216   fOutList->Add(fVertexZ[1]);
 217 }
 218
 219 //__________________________________________________________________________________________________
 220 void AliRsnAnalysisPhi7TeV::UserExec(Option_t *)
 221 {
 222 //
 223 // Main execution function.
 224 // Fills the fHEvents data member with the following legenda:
 225 // 0 -- event OK, prim vertex with tracks
 226 // 1 -- event OK, prim vertex with SPD
 227 // 2 -- event OK but vz large
 228 // 3 -- event bad
 229 //
 230
 231   static Int_t evNum = 0;
 232   evNum++;
 233
 234   // retrieve ESD event and related stack (if available)
 235   AliESDEvent *esd   = dynamic_cast<AliESDEvent*>(fInputEvent);
 236   AliStack    *stack = (fMCEvent ? fMCEvent->Stack() : 0x0);
 237
 238   // check the event
 239   Int_t eval = EventEval(esd);
 240   fHEvents->Fill(eval);
 241
 242   // if the event is good for analysis, process it
 243   if (eval == kGoodTracksPrimaryVertex || eval == kGoodSPDPrimaryVertex)
 244   {
 245     ProcessESD(esd, stack);
 246     ProcessMC(stack);
 247   }
 248
 249   // update histogram container
 250   PostData(3, fOutList);
 251 }
 252
 253 //__________________________________________________________________________________________________
 254 void AliRsnAnalysisPhi7TeV::Terminate(Option_t *)
 255 {
 256 //
 257 // Terminate
 258 //
 259 }
 260
 261 //__________________________________________________________________________________________________
 262 Int_t AliRsnAnalysisPhi7TeV::EventEval(AliESDEvent *esd)
 263 {
 264 //
 265 // Checks if the event is good for analysis.
 266 // Returns one of the flag values defined in the header
 267 //
 268
 269   static Int_t evNum = 0;
 270   evNum++;
 271
 272   // get the best primary vertex:
 273   // first try the one with tracks
 274   const AliESDVertex *vTrk  = esd->GetPrimaryVertexTracks();
 275   const AliESDVertex *vSPD  = esd->GetPrimaryVertexSPD();
 276   Double_t            vzTrk = 1000.0;
 277   Double_t            vzSPD = 1000.0;
 278   Int_t               ncTrk = -1;
 279   Int_t               ncSPD = -1;
 280   if (vTrk) ncTrk = (Int_t)vTrk->GetNContributors();
 281   if (vSPD) ncSPD = (Int_t)vSPD->GetNContributors();
 282   if (vTrk) vzTrk = TMath::Abs(vTrk->GetZv());
 283   if (vSPD) vzSPD = TMath::Abs(vSPD->GetZv());
 284   if(vTrk && ncTrk > 0)
 285   {
 286     // fill the histograms
 287     fVertexX[0]->Fill(vTrk->GetXv());
 288     fVertexY[0]->Fill(vTrk->GetYv());
 289     fVertexZ[0]->Fill(vTrk->GetZv());
 290
 291     // check VZ position
 292     if (vzTrk <= fMaxVz)
 293       return kGoodTracksPrimaryVertex;
 294     else
 295       return kFarTracksPrimaryVertex;
 296   }
 297   else if (vSPD && ncSPD > 0)
 298   {
 299     // fill the histograms
 300     fVertexX[1]->Fill(vSPD->GetXv());
 301     fVertexY[1]->Fill(vSPD->GetYv());
 302     fVertexZ[1]->Fill(vSPD->GetZv());
 303
 304     // check VZ position
 305     if (vzSPD <= fMaxVz)
 306       return kGoodSPDPrimaryVertex;
 307     else
 308       return kFarSPDPrimaryVertex;
 309   }
 310   else
 311     return kNoGoodPrimaryVertex;
 312 }
 313
 314 //__________________________________________________________________________________________________
 315 void AliRsnAnalysisPhi7TeV::ProcessESD
 316 (AliESDEvent *esd, AliStack *stack)
 317 {
 318 //
 319 // This function works with the ESD object
 320 //
 321
 322   static Int_t lastRun = -1;
 323
 324   // ITS stuff #1 create the response function
 325   Bool_t isMC = (stack != 0x0);
 326   AliITSPIDResponse itsrsp(isMC);
 327
 328   // TOF stuff #1: init OCDB
 329   Int_t run = esd->GetRunNumber();
 330   if (run != lastRun)
 331   {
 332     cout << "Run = " << run << " -- LAST = " << lastRun << endl;
 333     lastRun = run;
 334     AliCDBManager *cdb = AliCDBManager::Instance();
 335     cdb->SetDefaultStorage("raw://");
 336     cdb->SetRun(run);
 337     fTOFcalib->SetCorrectTExp(fTOFcorrectTExp);
 338     fTOFcalib->Init();
 339   }
 340   //AliCDBManager *cdb = AliCDBManager::Instance();
 341   //cdb->SetDefaultStorage("raw://");
 342   //cdb->SetRun(run);
 343   // TOF stuff #2: init calibration
 344   //fTOFcalib->SetCorrectTExp(fTOFcorrectTExp);
 345   //fTOFcalib->Init();
 346   // TOF stuff #3: calibrate
 347   if (fTOFcalibrateESD) fTOFcalib->CalibrateESD(esd);
 348   if (fTOFtuneMC) fTOFmaker->TuneForMC(esd);
 349   if (fTOFuseT0)
 350   {
 351     fTOFmaker->ComputeT0TOF(esd);
 352     fTOFmaker->ApplyT0TOF(esd);
 353     fESDpid->MakePID(esd, kFALSE, 0.);
 354   }
 355   // TOF stuff #4: define fixed function for compatibility range
 356   Double_t a1 = 0.01, a2 = -0.03;
 357   Double_t b1 = 0.25, b2 =  0.25;
 358   Double_t c1 = 0.05, c2 = -0.03;
 359   Double_t ymax, ymin;
 360
 361   // prepare to look on all tracks to select the ones
 362   // which pass all the cuts
 363   Int_t   ntracks = esd->GetNumberOfTracks();
 364   TArrayI pos(ntracks);
 365   TArrayI neg(ntracks);
 366   TArrayI itspos(ntracks);
 367   TArrayI itsneg(ntracks);
 368
 369   // loop on all tracks
 370   ULong_t  status;
 371   Int_t    i, k, charge, npos = 0, nneg = 0, nITS;
 372   Double_t times[10], tpcNSigma, tpcMaxNSigma, itsSignal, itsNSigma, mom, tofTime, tofSigma, tofRef, tofDiff, tofRel;
 373   Bool_t   okQuality, okTOF, okTPC, okITS, okTrack, isTPC, isITSSA;
 374   UChar_t  itsCluMap;
 375   for (i = 0; i < ntracks; i++)
 376   {
 377     AliESDtrack *track = esd->GetTrack(i);
 378     if (!track) continue;
 379
 380     // get commonly used variables
 381     status  = (ULong_t)track->GetStatus();
 382     mom     = track->P();
 383     isTPC   = ((status & AliESDtrack::kTPCin)  != 0);
 384     isITSSA = ((status & AliESDtrack::kTPCin)  == 0 && (status & AliESDtrack::kITSrefit) != 0 && (status & AliESDtrack::kITSpureSA) == 0 && (status & AliESDtrack::kITSpid) != 0);
 385     //cout << isTPC << ' ' << isITSSA << endl;
 386
 387     // accept only tracks which are TPC+ITS or ITS standalone
 388     if (!isTPC && !isITSSA) continue;
 389
 390     // define selection properties depending on track type
 391     // it track is standard TPC+ITS+TOF track, check standard cuts and TOF
 392     // if the track is an ITS standalone, check its specific cuts only
 393     okTrack = kTRUE;
 394
 395     if (isTPC)
 396     {
 397       // check standard ESD cuts
 398       okQuality = fESDtrackCutsTPC.IsSelected(track);
 399       //cout << "GLOBAL -- quality = " << (okQuality ? "GOOD" : "BAD") << endl;
 400       if (!okQuality) continue;
 401
 402       // check TPC dE/dx
 403       if (fCheckTPC)
 404       {
 405         tpcNSigma = TMath::Abs(fESDpid->NumberOfSigmasTPC(track, AliPID::kKaon));
 406         if (track->GetInnerParam()->P() > fTPCpLimit) tpcMaxNSigma = fMinTPCband; else tpcMaxNSigma = fMaxTPCband;
 407         okTPC = (tpcNSigma <= tpcMaxNSigma);
 408         //cout << "ALT -- TPC    -- nsigma = " << tpcNSigma << ", max = " << tpcMaxNSigma << " --> " << (okTPC ? "OK" : "FAILED") << endl;
 409         //cout << "ALTTPC -- " << fTPCpar[0] << ' ' << fTPCpar[1] << ' ' << fTPCpar[2] << ' ' << fTPCpar[3] << ' ' << fTPCpar[4] << endl;
 410       }
 411       else
 412       {
 413         okTPC = kTRUE;
 414       }
 415
 416       // check TOF (only if momentum is large than function asymptote and flags are OK)
 417       if (fCheckTOF)
 418       {
 419         if (((status & AliESDtrack::kTOFout) != 0) && ((status & AliESDtrack::kTIME) != 0) && mom > TMath::Max(b1, b2))
 420         {
 421           track->GetIntegratedTimes(times);
 422           tofTime  = (Double_t)track->GetTOFsignal();
 423           tofSigma = fTOFmaker->GetExpectedSigma(mom, times[AliPID::kKaon], AliPID::ParticleMass(AliPID::kKaon));
 424           tofRef   = times[AliPID::kKaon];
 425           if (tofRef > 0.0)
 426           {
 427             tofDiff  = (tofTime - tofRef);
 428             tofRel   = (tofTime - tofRef) / tofRef;
 429             ymax     = a1 / (mom - b1) + c1;
 430             ymin     = a2 / (mom - b2) + c2;
 431             okTOF    = (tofRel >= ymin && tofRel <= ymax);
 432             //cout << "TOF    -- diff = " << tofDiff << ", rel diff = " << tofRel << ", range = " << ymin << " to " << ymax << ", sigma = " << tofSigma << " --> " << (okTOF ? "OK" : "FAILED") << endl;
 433           }
 434           else
 435           {
 436             okTOF = kTRUE;
 437             //cout << "TOF    -- not checked due to ZERO reference time" << endl;
 438           }
 439         }
 440         else
 441         {
 442           okTOF = kTRUE;
 443           //cout << "TOF    -- not checked because TOF pid absent" << endl;
 444         }
 445       }
 446       else
 447       {
 448         okTOF = kTRUE;
 449       }
 450
 451       okTrack = okQuality && okTPC && okTOF;
 452       //cout << "GLOBAL -- overall = " << (okTrack ? "ACCEPTED" : "REJECTED") << endl;
 453     }
 454     else
 455     {
 456       // check standard ESD cuts
 457       okQuality = fESDtrackCutsITS.IsSelected(track);
 458       //cout << "ITSSA  -- quality = " << (okQuality ? "GOOD" : "BAD") << endl;
 459       if (!okQuality) continue;
 460
 461       // check dE/dx
 462       if (fCheckITS)
 463       {
 464         itsSignal = track->GetITSsignal();
 465         itsCluMap = track->GetITSClusterMap();
 466         nITS      = 0;
 467         for(k = 2; k < 6; k++) if(itsCluMap & (1 << k)) ++nITS;
 468         if (nITS < 3)
 469         {
 470           okITS = kFALSE;
 471           //cout << "ITS    -- not checked due to too few PID clusters" << endl;
 472         }
 473         else
 474         {
 475           itsNSigma = itsrsp.GetNumberOfSigmas(mom, itsSignal, AliPID::kKaon, nITS, kTRUE);
 476           okITS = (TMath::Abs(itsNSigma) <= fMaxITSband);
 477           //cout << "ITS    -- nsigma = " << itsNSigma << ", max = " << fMaxITSband << " --> " << (okITS ? "OK" : "FAILED") << endl;
 478         }
 479       }
 480       else
 481       {
 482         okITS = kTRUE;
 483       }
 484
 485       okTrack = okQuality && okITS;
 486       //cout << "ITSSA  -- overall = " << (okTrack ? "ACCEPTED" : "REJECTED") << endl;
 487     }
 488
 489     // skip tracks not passing cuts
 490     if (!okTrack) continue;
 491
 492     // if all checks are passed, add the track index in one of the
 493     // charged tracks arrays
 494     charge = (Int_t)track->Charge();
 495     if (charge > 0)
 496     {
 497       pos[npos] = i;
 498       if (isITSSA) itspos[npos] = 1; else itspos[npos] = 0;
 499       npos++;
 500     }
 501     else if (charge < 0)
 502     {
 503       neg[nneg] = i;
 504       if (isITSSA) itsneg[nneg] = 1; else itsneg[nneg] = 0;
 505       nneg++;
 506     }
 507   }
 508
 509   // resize arrays accordingly
 510   pos.Set(npos);
 511   neg.Set(nneg);
 512   itspos.Set(npos);
 513   itsneg.Set(nneg);
 514
 515   /*
 516   // fill unlike-sign
 517   Int_t i1, i2;
 518   for (i1 = 0; i1 < npos; i1++)
 519     for (i2 = 0; i2 < nneg; i2++)
 520       AddEntryFromESD(esd, i1, i2, itspos[i1], itsneg[i2], 0, stack);
 521
 522   // fill ++ like-sign
 523   for (i1 = 0; i1 < npos; i1++)
 524     for (i2 = i1 + 1; i2 < npos; i2++)
 525       AddEntryFromESD(esd, i1, i2, itspos[i1], itspos[i2], 1, 0x0);
 526
 527   // fill -- like-sign
 528   for (i1 = 0; i1 < nneg; i1++)
 529     for (i2 = i1 + 1; i2 < nneg; i2++)
 530       AddEntryFromESD(esd, i1, i2, itsneg[i1], itsneg[i2], -1, 0x0);
 531       */
 532
 533   // loop on unlike-sign pairs to compute invariant mass signal
 534   Int_t           ip, in, lp, ln;
 535   AliPID          pid;
 536   Double_t        kmass = pid.ParticleMass(AliPID::kKaon);
 537   Double_t        phimass = 1.019455;
 538   TParticle      *partp = 0x0, *partn = 0x0;
 539   AliESDtrack    *tp = 0x0, *tn = 0x0;
 540   TLorentzVector  vp, vn, vsum, vref;
 541   for (ip = 0; ip < npos; ip++)
 542   {
 543     tp = esd->GetTrack(pos[ip]);
 544     lp = TMath::Abs(tp->GetLabel());
 545     if (stack) partp = stack->Particle(lp);
 546
 547     for (in = 0; in < nneg; in++)
 548     {
 549       if (pos[ip] == neg[in])
 550       {
 551         AliError("POS = NEG");
 552         continue;
 553       }
 554       tn = esd->GetTrack(neg[in]);
 555       ln = TMath::Abs(tn->GetLabel());
 556       if (stack) partn = stack->Particle(ln);
 557
 558       fPDG = 0;
 559       if (partp && partn)
 560       {
 561         if (partp->GetFirstMother() == partn->GetFirstMother())
 562         {
 563           if (partp->GetFirstMother() > 0)
 564           {
 565             TParticle *mum = stack->Particle(partp->GetFirstMother());
 566             fPDG = mum->GetPdgCode();
 567           }
 568         }
 569       }
 570       fPDG = TMath::Abs(fPDG);
 571
 572       vp.SetXYZM(tp->Px(), tp->Py(), tp->Pz(), kmass);
 573       vn.SetXYZM(tn->Px(), tn->Py(), tn->Pz(), kmass);
 574       vsum = vp + vn;
 575       vref.SetXYZM(vsum.X(), vsum.Y(), vsum.Z(), phimass);
 576
 577       fCh     = 0;
 578       fIM     = (Float_t)vsum.M();
 579       fPt     = (Float_t)vsum.Perp();
 580       fEta    = (Float_t)vsum.Eta();
 581       fY      = (Float_t)vref.Rapidity();
 582       fITS[0] = itspos[ip];
 583       fITS[1] = itsneg[in];
 584
 585       //if (fIM < 0.9 || fIM >  1.4) continue;
 586       //if (fPt < 0.0 || fPt > 20.0) continue;
 587
 588       fRsnTreeComp->Fill();
 589     }
 590   }
 591
 592   // loop on like-sign pairs to compute invariant mass background
 593   Int_t           i1, i2;
 594   AliESDtrack    *t1 = 0x0, *t2 = 0x0;
 595   TLorentzVector  v1, v2;
 596
 597   // pos-pos
 598   for (i1 = 0; i1 < npos; i1++)
 599   {
 600     t1 = esd->GetTrack(pos[i1]);
 601
 602     for (i2 = i1+1; i2 < npos; i2++)
 603     {
 604       t2 = esd->GetTrack(pos[i2]);
 605
 606       v1.SetXYZM(t1->Px(), t1->Py(), t1->Pz(), kmass);
 607       v2.SetXYZM(t2->Px(), t2->Py(), t2->Pz(), kmass);
 608       vsum = v1 + v2;
 609       vref.SetXYZM(vsum.X(), vsum.Y(), vsum.Z(), phimass);
 610
 611       fPDG    = 0;
 612       fCh     = 1;
 613       fIM     = (Float_t)vsum.M();
 614       fPt     = (Float_t)vsum.Perp();
 615       fEta    = (Float_t)vsum.Eta();
 616       fY      = (Float_t)vref.Rapidity();
 617       fITS[0] = itspos[i1];
 618       fITS[1] = itspos[i2];
 619
 620       //if (fIM < 0.9 || fIM >  1.4) continue;
 621       //if (fPt < 0.0 || fPt > 20.0) continue;
 622
 623       fRsnTreeComp->Fill();
 624     }
 625   }
 626   // neg-neg
 627   for (i1 = 0; i1 < nneg; i1++)
 628   {
 629     t1 = esd->GetTrack(neg[i1]);
 630
 631     for (i2 = i1+1; i2 < nneg; i2++)
 632     {
 633       t2 = esd->GetTrack(neg[i2]);
 634
 635       v1.SetXYZM(t1->Px(), t1->Py(), t1->Pz(), kmass);
 636       v2.SetXYZM(t2->Px(), t2->Py(), t2->Pz(), kmass);
 637       vsum = v1 + v2;
 638       vref.SetXYZM(vsum.X(), vsum.Y(), vsum.Z(), phimass);
 639
 640       fPDG    = 0;
 641       fCh     = -1;
 642       fIM     = (Float_t)vsum.M();
 643       fPt     = (Float_t)vsum.Perp();
 644       fEta    = (Float_t)vsum.Eta();
 645       fY      = (Float_t)vref.Rapidity();
 646       fITS[0] = itsneg[i1];
 647       fITS[1] = itsneg[i2];
 648
 649       //if (fIM < 0.9 || fIM >  1.4) continue;
 650       //if (fPt < 0.0 || fPt > 20.0) continue;
 651
 652       fRsnTreeComp->Fill();
 653     }
 654   }
 655
 656   PostData(1, fRsnTreeComp);
 657 }
 658
 659 //__________________________________________________________________________________________________
 660 void AliRsnAnalysisPhi7TeV::AddEntryFromESD(AliESDEvent *esd, Int_t i1, Int_t i2, Int_t its1, Int_t its2, Short_t charge, AliStack *stack)
 661 {
 662 //
 663 // Add an entry to the output TTree computed from two tracks in the ESD
 664 //
 665
 666   // loop on unlike-sign pairs to compute invariant mass signal
 667   AliPID          pid;
 668   Double_t        kmass = pid.ParticleMass(AliPID::kKaon);
 669   Double_t        phimass = 1.019455;
 670   TLorentzVector  v1, v2, vsum, vref;
 671
 672   AliESDtrack *t1 = esd->GetTrack(i1);
 673   AliESDtrack *t2 = esd->GetTrack(i2);
 674
 675   v1.SetXYZM(t1->Px(), t1->Py(), t1->Pz(), kmass);
 676   v2.SetXYZM(t2->Px(), t2->Py(), t2->Pz(), kmass);
 677   vsum = v1 + v2;
 678   vref.SetXYZM(vsum.X(), vsum.Y(), vsum.Z(), phimass);
 679
 680   // if stack is present, search for true pair
 681   fPDG = 0;
 682   Int_t l1 = TMath::Abs(t1->GetLabel());
 683   Int_t l2 = TMath::Abs(t2->GetLabel());
 684   if (stack)
 685   {
 686     TParticle *part1 = stack->Particle(l1);
 687     TParticle *part2 = stack->Particle(l2);
 688     if (part1 && part2)
 689     {
 690       if (part1->GetFirstMother() == part2->GetFirstMother())
 691       {
 692         if (part1->GetFirstMother() > 0)
 693         {
 694           TParticle *mum = stack->Particle(part1->GetFirstMother());
 695           fPDG = mum->GetPdgCode();
 696         }
 697       }
 698     }
 699   }
 700
 701   fPDG = TMath::Abs(fPDG);
 702   fCh     = charge;
 703   fIM     = (Float_t)vsum.M();
 704   fPt     = (Float_t)vsum.Perp();
 705   fEta    = (Float_t)vsum.Eta();
 706   fY      = (Float_t)vref.Rapidity();
 707   fITS[0] = its1;
 708   fITS[1] = its2;
 709
 710   //if (fIM < 0.9 || fIM >  1.4) return;
 711   //if (fPt < 0.0 || fPt > 20.0) return;
 712
 713   fRsnTreeComp->Fill();
 714 }
 715
 716 //__________________________________________________________________________________________________
 717 void AliRsnAnalysisPhi7TeV::ProcessMC(AliStack *stack)
 718 {
 719 //
 720 // Function to process stack only
 721 //
 722
 723   if (!stack) return;
 724   Int_t nPart = stack->GetNtrack();
 725
 726   // loop to compute invariant mass
 727   Int_t           ip, in;
 728   AliPID          pid;
 729   Double_t        kmass = pid.ParticleMass(AliPID::kKaon);
 730   Double_t        phimass = 1.019455;
 731   TParticle      *partp = 0x0, *partn = 0x0;
 732   TLorentzVector  vp, vn, vsum, vref;
 733
 734   for (ip = 0; ip < nPart; ip++)
 735   {
 736     partp = stack->Particle(ip);
 737     if (partp->GetPdgCode() != 321) continue;
 738
 739     for (in = 0; in < nPart; in++)
 740     {
 741       partn = stack->Particle(in);
 742       if (partn->GetPdgCode() != -321) continue;
 743
 744       fPDG = 0;
 745       if (partp->GetFirstMother() == partn->GetFirstMother())
 746       {
 747         if (partp->GetFirstMother() > 0)
 748         {
 749           TParticle *mum = stack->Particle(partp->GetFirstMother());
 750           fPDG = mum->GetPdgCode();
 751         }
 752       }
 753       fPDG = TMath::Abs(fPDG);
 754       if (fPDG != 333) continue;
 755
 756       vp.SetXYZM(partp->Px(), partp->Py(), partp->Pz(), kmass);
 757       vn.SetXYZM(partn->Px(), partn->Py(), partn->Pz(), kmass);
 758       vsum = vp + vn;
 759       vref.SetXYZM(vsum.X(), vsum.Y(), vsum.Z(), phimass);
 760
 761       fIM  = (Float_t)vsum.M();
 762       fPt  = (Float_t)vsum.Perp();
 763       fEta = (Float_t)vsum.Eta();
 764       fY   = (Float_t)vref.Rapidity();
 765
 766       fRsnTreeTrue->Fill();
 767     }
 768   }
 769
 770   PostData(2, fRsnTreeTrue);
 771 }