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