PWG3/muon/AliCFMuonResUpsilon.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //-----------------------------------------------------------------------
  17 // Example of task (running locally, on AliEn and CAF),
  18 // which provides standard way of calculating acceptance and efficiency
  19 // between different steps of the procedure.
  20 // The ouptut of the task is a AliCFContainer from which the efficiencies
  21 // can be calculated
  22 //-----------------------------------------------------------------------
  23 // Author : R. Vernet, Consorzio Cometa - Catania (it)
  24 //-----------------------------------------------------------------------
  25 // Modification done by X. Lopez - LPC Clermont (fr)
  26 //-----------------------------------------------------------------------
  27 // Modification done by S. Ahn - LPC Clermont (fr), Konkuk University (kr)
  28 //-----------------------------------------------------------------------
  29
  30
  31
  32 #ifndef ALICFMUONRESUPSILON_CXX
  33 #define ALICFMUONRESUPSILON_CXX
  34
  35 #include "TH1.h"
  36 #include "TParticle.h"
  37 #include "TChain.h"
  38 #include "TLorentzVector.h"
  39 #include "TCanvas.h"
  40
  41 #include "AliHeader.h"
  42 #include "AliStack.h"
  43 #include "AliMCEvent.h"
  44 #include "AliAnalysisManager.h"
  45 #include "AliLog.h"
  46
  47 #include "AliVEvent.h"
  48 #include "AliPhysicsSelection.h"
  49
  50 #include "AliCFMuonResUpsilon.h"
  51 #include "AliCFManager.h"
  52 #include "AliCFCutBase.h"
  53 #include "AliCFContainer.h"
  54
  55 #include "AliESDEvent.h"
  56 #include "AliESDHeader.h"
  57 #include "AliESDtrack.h"
  58 #include "AliESDMuonTrack.h"
  59 #include "AliESDtrack.h"
  60 #include "AliESDInputHandler.h"
  61
  62 #include "AliAODEvent.h"
  63 #include "AliAODInputHandler.h"
  64 #include "AliAODMCParticle.h"
  65
  66 ClassImp(AliCFMuonResUpsilon)
  67
  68 //__________________________________________________________________________
  69 AliCFMuonResUpsilon::AliCFMuonResUpsilon() :
  70         AliAnalysisTaskSE(""),
  71         fReadAODData(kFALSE),
  72         fReadMCInfo(kFALSE),
  73   fCFManager(0x0),
  74         fnevts(0x0),
  75         fIsPhysSelMB(kFALSE),
  76         fIsPhysSelMUON(kFALSE),
  77   fQAHistList(0x0),
  78         fPDG(0),
  79         fPtMin(0),
  80         fPtMax(0),
  81         fYMin(0),
  82         fYMax(0),
  83         fTrigClassMuon(""),
  84         fTrigClassInteraction(""),
  85         fDistinguishTrigClass(kTRUE)
  86 {
  87
  88   //Default ctor
  89
  90 }
  91 //___________________________________________________________________________
  92 AliCFMuonResUpsilon::AliCFMuonResUpsilon(const Char_t* name) :
  93   AliAnalysisTaskSE(name),
  94         fReadAODData(kFALSE),
  95         fReadMCInfo(kFALSE),
  96   fCFManager(0x0),
  97         fnevts(0x0),
  98         fIsPhysSelMB(kFALSE),
  99         fIsPhysSelMUON(kFALSE),
 100   fQAHistList(0x0),
 101         fPDG(0),
 102         fPtMin(0),
 103         fPtMax(0),
 104         fYMin(0),
 105         fYMax(0),
 106         fTrigClassMuon(""),
 107         fTrigClassInteraction(""),
 108         fDistinguishTrigClass(kTRUE)
 109 {
 110
 111   // Constructor. Initialization of Inputs and Outputs
 112
 113   Info("AliCFMuonResUpsilon","Calling Constructor");
 114
 115         fnevts  = new TH1D("fnevts","fnevts",5,0,5);                                                    // nevent CINT1B,CMUS1B + PhysSel
 116
 117         TString nameside[3]={"AC","B","E"};
 118
 119   SetTrigClassMuonName();
 120   SetTrigClassInteracName();
 121   SetTrigClassSideName(nameside);
 122
 123         DefineInput(0, TChain::Class());
 124         DefineOutput(1,TH1D::Class());
 125   DefineOutput(2,AliCFContainer::Class());
 126 }
 127
 128 //___________________________________________________________________________
 129 AliCFMuonResUpsilon& AliCFMuonResUpsilon::operator=(const AliCFMuonResUpsilon& c)
 130 {
 131
 132   // Assignment operator
 133
 134   if (this!=&c) {
 135     AliAnalysisTaskSE::operator=(c);
 136     fReadAODData = c.fReadAODData;
 137                 fReadMCInfo = c.fReadMCInfo;
 138     fCFManager  = c.fCFManager;
 139     fQAHistList = c.fQAHistList;
 140                 fnevts = c.fnevts;
 141                 fPDG = c.fPDG;
 142                 fPtMin = c.fPtMin;
 143                 fPtMax = c.fPtMax;
 144                 fYMin = c.fYMin;
 145                 fYMax = c.fYMax;
 146   }
 147   return *this;
 148 }
 149
 150 //___________________________________________________________________________
 151 AliCFMuonResUpsilon::AliCFMuonResUpsilon(const AliCFMuonResUpsilon& c) :
 152   AliAnalysisTaskSE(c),
 153   fReadAODData(c.fReadAODData),
 154         fReadMCInfo(c.fReadMCInfo),
 155   fCFManager(c.fCFManager),
 156         fnevts(c.fnevts),
 157         fIsPhysSelMB(kFALSE),
 158         fIsPhysSelMUON(kFALSE),
 159   fQAHistList(c.fQAHistList),
 160         fPDG(c.fPDG),
 161         fPtMin(c.fPtMin),
 162         fPtMax(c.fPtMax),
 163         fYMin(c.fYMin),
 164         fYMax(c.fYMax),
 165   fTrigClassMuon(c.fTrigClassMuon),
 166   fTrigClassInteraction(c.fTrigClassInteraction),
 167   fDistinguishTrigClass(c.fDistinguishTrigClass)
 168 {
 169
 170   // Copy Constructor
 171
 172 }
 173
 174 //___________________________________________________________________________
 175 AliCFMuonResUpsilon::~AliCFMuonResUpsilon() {
 176
 177   //destructor
 178
 179   Info("~AliCFMuonResUpsilon","Calling Destructor");
 180   if (fCFManager)           delete fCFManager ;
 181         if (fnevts)                                                             delete fnevts ;
 182   if (fQAHistList)                                      {fQAHistList->Clear(); delete fQAHistList;}
 183 }
 184 //___________________________________________________________________________
 185 void AliCFMuonResUpsilon::UserCreateOutputObjects() {
 186         // UserCreateOutputObjects
 187
 188         //TH1D *fnevts  = new TH1D("fnevts","fnevts",5,0,5);                                                    // nevent CINT1B,CMUS1B + PhysSel
 189
 190   PostData(1,fnevts) ;
 191   PostData(2,fCFManager->GetParticleContainer()) ;
 192 }
 193
 194 //_________________________________________________
 195 void AliCFMuonResUpsilon::UserExec(Option_t *)
 196 {
 197
 198   // Main loop function
 199
 200         // variables
 201         fIsPhysSelMB=kFALSE;                    // physics selection : MB
 202         fIsPhysSelMUON=kFALSE;          // physics selection : MUON
 203
 204   Double_t containerInput[14] = {0,} ;
 205
 206   fnevts->Fill(0.5);
 207         containerInput[0] = 0.5;
 208
 209         if(!fReadAODData) {     // ESD-based ANALYSIS
 210
 211                 Bool_t flag=kFALSE;
 212
 213                 if(fReadMCInfo) {
 214                 if (!fMCEvent) {
 215                 Error("UserExec","NO MC EVENT FOUND!");
 216                 return;
 217                 }
 218
 219                         // MC part ----------------------------------------------------------------------------------
 220                         fCFManager->SetMCEventInfo(fMCEvent);
 221                 AliStack *stack = fMCEvent->Stack();
 222                         Int_t npart=fMCEvent->GetNumberOfTracks();
 223
 224                         //printf("ESD: npart=%d\n", npart);
 225
 226                 for (Int_t ipart=0; ipart<npart; ipart++) {
 227
 228                                 AliMCParticle *mcPart  = (AliMCParticle*)fMCEvent->GetTrack(ipart);
 229
 230                 // Mother kinematics
 231                         TParticle *part = mcPart->Particle();
 232                 Double_t e = part->Energy();
 233                 Double_t pz = part->Pz();
 234                 Double_t rapmc = Rap(e,pz);
 235
 236                                 // Selection of the resonance
 237                 if(!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
 238                                 if(part->Pt()<0 || part->Pt()>1000) continue;
 239                                 if(rapmc<-4 || rapmc>-2.4) continue;
 240
 241                                 // Decays kinematics
 242                 Int_t p0 = part->GetDaughter(0);
 243                 TParticle *part0 = stack->Particle(p0);
 244                                 AliMCParticle *mcpart0 = new AliMCParticle(part0);
 245                 Int_t pdg0 = part0->GetPdgCode();
 246
 247                 Int_t p1 = part->GetDaughter(1);
 248                 TParticle *part1 = stack->Particle(p1);
 249                                 AliMCParticle *mcpart1 = new AliMCParticle(part1);
 250                 Int_t pdg1 = part1->GetPdgCode();
 251
 252                                 Double_t e0 = part0->Energy();
 253                                 Double_t pz0 = part0->Pz();
 254                                 Double_t py0 = part0->Py();
 255                                 Double_t px0 = part0->Px();
 256                                 Double_t charge0 = part0->GetPDG()->Charge()/3;
 257                                 Double_t pt0 = TMath::Sqrt(px0*px0+py0*py0);
 258                                 Double_t theta0 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(px0*px0+py0*py0),pz0);
 259                                 Double_t eta0 = part0->Eta();
 260                                 Double_t vz0 = part0->Vz();
 261                                 Double_t mom0 = part0->P();
 262
 263                 if(!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcpart0)) continue;
 264                                 if(pt0<0.5) continue;
 265                                 if(theta0<171. || theta0>178.) continue;
 266
 267                                 Double_t e1 = part1->Energy();
 268                                 Double_t pz1 = part1->Pz();
 269                                 Double_t py1 = part1->Py();
 270                                 Double_t px1 = part1->Px();
 271                                 Double_t charge1 = part1->GetPDG()->Charge()/3;
 272                                 Double_t pt1 = TMath::Sqrt(px1*px1+py1*py1);
 273                                 Double_t theta1 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(px1*px1+py1*py1),pz1);
 274                                 Double_t eta1 = part1->Eta();
 275                                 Double_t vz1 = part1->Vz();
 276                                 Double_t mom1 = part1->P();
 277
 278                 if(!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcpart1)) continue;
 279                                 if(pt1<0.5) continue;
 280                                 if(theta1<171. || theta1>178.) continue;
 281
 282                     if(pdg0==-13 || pdg1==-13) {
 283
 284                                         Double_t ptmc = TMath::Sqrt((px0+px1)*(px0+px1)+(py0+py1)*(py0+py1));
 285                                         Double_t imassmc = Imass(e0,px0,py0,pz0,e1,px1,py1,pz1);
 286                                         if(!imassmc) continue;
 287
 288                                         containerInput[1] = rapmc ;
 289                                         containerInput[2] = ptmc ;
 290                                         containerInput[3] = imassmc ;
 291                                         containerInput[4] = 1. ;
 292                                         containerInput[5] = pt0 ;
 293                                         containerInput[5] = pt1 ;
 294                                         containerInput[6] = mom0 ;
 295                                         containerInput[6] = mom1 ;
 296                                         containerInput[7] = 1. ;
 297                                         containerInput[8] = 0. ;
 298                                         containerInput[8] = 0. ;
 299                                         containerInput[9] = charge0+charge1 ;
 300                                         containerInput[10] = eta0;
 301                                         containerInput[10] = eta1;
 302                                         containerInput[11] = theta0;
 303                                         containerInput[11] = theta1;
 304                                         containerInput[12] = vz0;
 305                                         containerInput[12] = vz1;
 306                                         containerInput[13] = 0;
 307                                         containerInput[13] = 0;
 308
 309                                         // fill the container at the first step
 310                                         fCFManager->GetParticleContainer()->Fill(containerInput,0);             // MC container
 311                                         flag=kTRUE;
 312                 }       // second mu loop
 313                         }       // first mu loop
 314                 } // end fReadMCInfo
 315
 316                 if(fReadMCInfo && !flag) return;
 317                 // RECO : ESD part ----------------------------------------------------------------------------------
 318
 319                 // ESD handler
 320         AliESDEvent *fESD = 0x0;
 321
 322                 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 323                 if( ! esdH) {
 324                         AliError("Cannot get input event handler");
 325                         return;
 326                 }
 327         fESD = esdH->GetEvent();
 328
 329
 330                 Int_t ntrk=fESD->GetNumberOfMuonTracks();
 331
 332                 Int_t trigfired=-1;
 333                 Int_t trigside=-1;
 334
 335                 if(!fReadMCInfo) {
 336                         fIsPhysSelMB=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMB));
 337                         fIsPhysSelMUON=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMUON));
 338
 339                         if(fDistinguishTrigClass) {
 340                                 TString trigclass = fESD->GetFiredTriggerClasses();
 341                                 if(trigclass.Contains(fTrigClassMuon)) trigfired = 1;
 342                                 else if(trigclass.Contains(fTrigClassInteraction)) trigfired = 0;
 343                                 for(Int_t i=0;i<3;i++) {
 344                                         if(trigfired==1 && trigclass.Contains(fTrigClassSide[i])) {trigside=i;}
 345                                         if(trigfired==0 && trigclass.Contains(fTrigClassSide[i])) {trigside=i;}
 346                                 }
 347                         }
 348
 349                 if(trigside==1) {       // Beam-Beam event
 350                                 if(trigfired) { // CMUS1-B
 351                                         containerInput[0]=3.5;
 352                                         fnevts->Fill(3.5);
 353                                         if(fIsPhysSelMUON) {
 354                                                 containerInput[0]=4.5;
 355                                                 fnevts->Fill(4.5);      // PhysSel
 356                                         }
 357                                 } else {        // CINT1-B
 358                                         containerInput[0]=1.5;
 359                                         fnevts->Fill(1.5);
 360                                         if(fIsPhysSelMB) {
 361                                                 containerInput[0]=2.5;
 362                                                 fnevts->Fill(2.5);      // PhysSel
 363                                         }
 364                                 }
 365                         }
 366                 }
 367                 else trigside = 0;      // for MC
 368
 369                 for(Int_t j=0; j<ntrk-1; j++) {
 370                         for(Int_t jj=j+1; jj<ntrk; jj++) {
 371                                 AliESDMuonTrack *mu1 = new AliESDMuonTrack(*(fESD->GetMuonTrack(j)));
 372                                 AliESDMuonTrack *mu2 = new AliESDMuonTrack(*(fESD->GetMuonTrack(jj)));
 373
 374                                 Double_t trigCondition = 0;
 375                                 if (mu1->GetMatchTrigger()>=mu2->GetMatchTrigger()) trigCondition = mu1->GetMatchTrigger()+0.1*mu2->GetMatchTrigger();
 376                                 else trigCondition = mu2->GetMatchTrigger()+0.1*mu1->GetMatchTrigger();
 377
 378                                 // mu1
 379                                 Double_t zr1 = mu1->Charge();
 380                                 Double_t pxr1 = mu1->Px();
 381                                 Double_t pyr1 = mu1->Py();
 382                                 Double_t pzr1 = mu1->Pz();
 383                 Double_t ptr1 = TMath::Sqrt(pxr1*pxr1+pyr1*pyr1);
 384                                 Double_t er1 = mu1->E();
 385                                 Double_t momr1 = mu1->P();
 386                                 Double_t vzr1 = mu1->GetZ();
 387                                 Double_t dcar1 = mu1->GetDCA();
 388                                 Double_t rabsr1 = mu1->GetRAtAbsorberEnd();
 389                                 Double_t thetaabsr1 = 180.*(1.-TMath::ATan(rabsr1/505.)/TMath::Pi());  // [deg]
 390                                 Double_t thetaabsr12 = (thetaabsr1*TMath::Pi()/180.)/2.;
 391                                 Double_t tanthetaabsr1 = TMath::Tan(thetaabsr12);
 392                                 Double_t etar1 = -TMath::Log(tanthetaabsr1);
 393
 394                                 if(Rap(er1,pzr1)<-4 || Rap(er1,pzr1)>-2.5) continue;
 395
 396                                 // mu2
 397                                 Double_t zr2 = mu2->Charge();
 398                                 Double_t pxr2 = mu2->Px();
 399                                 Double_t pyr2 = mu2->Py();
 400                                 Double_t pzr2 = mu2->Pz();
 401                 Double_t ptr2 = TMath::Sqrt(pxr2*pxr2+pyr2*pyr2);
 402                                 Double_t er2 = mu2->E();
 403                                 Double_t momr2 = mu2->P();
 404                                 Double_t vzr2 = mu2->GetZ();
 405                                 Double_t dcar2 = mu2->GetDCA();
 406                                 Double_t rabsr2 = mu2->GetRAtAbsorberEnd();
 407                                 Double_t thetaabsr2 = 180.*(1.-TMath::ATan(rabsr2/505.)/TMath::Pi());  // [deg]
 408                                 Double_t thetaabsr22 = (thetaabsr2*TMath::Pi()/180.)/2.;
 409                                 Double_t tanthetaabsr2 = TMath::Tan(thetaabsr22);
 410                                 Double_t etar2 = -TMath::Log(tanthetaabsr2);
 411
 412                                 if(Rap(er2,pzr2)<-4 || Rap(er2,pzr2)>-2.5) continue;
 413
 414                                 // for MC
 415                                 if(fReadMCInfo) {
 416                                         //mu1
 417                                         rabsr1 = 0;
 418                                         Double_t thetar1 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(pxr1*pxr1+pyr1*pyr1),pzr1);
 419                                         thetaabsr1 = thetar1;
 420                                         etar1 = mu1->Eta();
 421                                         //mu2
 422                                         rabsr2 = 0;
 423                                         Double_t thetar2 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(pxr2*pxr2+pyr2*pyr2),pzr2);
 424                                         thetaabsr2 = thetar2;
 425                                         etar2 = mu2->Eta();
 426                                 }
 427
 428                                 if(TMath::Abs(etar1) > 8 || TMath::Abs(etar2) > 8) continue;
 429
 430                                 // dimuon
 431                                 Double_t ptrec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2));
 432                                 Double_t raprec= Rap((er1+er2),(pzr1+pzr2));
 433                                 Double_t imassrec = Imass(er1,pxr1,pyr1,pzr1,er2,pxr2,pyr2,pzr2);
 434                                 if(!imassrec) continue;
 435
 436                                 containerInput[1] = raprec ;
 437                                 containerInput[2] = ptrec ;
 438                                 containerInput[3] = imassrec ;
 439                                 containerInput[4] = trigCondition+0.05 ;
 440                                 containerInput[5] = ptr1 ;
 441                                 containerInput[5] = ptr2 ;
 442                                 containerInput[6] = momr1;
 443                                 containerInput[6] = momr2;
 444                                 containerInput[7] = trigside;
 445                                 containerInput[8] = rabsr1;
 446                                 containerInput[8] = rabsr2;
 447                                 containerInput[9] = zr1 + zr2;
 448                                 containerInput[10] = etar1;
 449                                 containerInput[10] = etar2;
 450                                 containerInput[11] = thetaabsr1;
 451                                 containerInput[11] = thetaabsr2;
 452                                 containerInput[12] = vzr1;
 453                                 containerInput[12] = vzr2;
 454                                 containerInput[13] = dcar1;
 455                                 containerInput[13] = dcar2;
 456
 457                         if (trigfired==1 && trigside==1) fCFManager->GetParticleContainer()->Fill(containerInput,4);
 458                         if (trigfired==0 && trigside==1) fCFManager->GetParticleContainer()->Fill(containerInput,1);
 459                                 if(fReadMCInfo) fCFManager->GetParticleContainer()->Fill(containerInput,1); // Rec container
 460                         }       // second mu loop
 461                 }       // first mu loop
 462         } // end ESD-based ANALYSIS
 463
 464         else {  // AOD-based ANALYSIS
 465
 466                 Bool_t flag=kTRUE;
 467
 468                 // AOD handler
 469                 AliAODEvent *fAOD = 0x0;
 470
 471         AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 472                 if( ! aodH) {
 473                         AliError("Cannot get input event handler");
 474                         return;
 475                 }
 476                 fAOD = aodH->GetEvent();
 477
 478                 // MC part -----------------------------------------------------------------------------------
 479
 480                 if(fReadMCInfo) {
 481                         TClonesArray *mcArr = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
 482                         if( ! mcArr) {
 483                           AliError("Cannot get MC innf in AOD MC branch");
 484                           return;
 485                         }
 486                         Int_t npart=mcArr->GetEntries();
 487
 488                         for(Int_t i=0; i<npart; i++) {
 489                                 AliAODMCParticle *mctrack = (AliAODMCParticle*) mcArr->At(i);
 490                                 // select resonance
 491                                 if(mctrack->GetPdgCode()!=fPDG) continue;
 492                                 // cuts on resonance
 493                                 if(!(mctrack->Pt()>0 && mctrack->Pt()<1000) || !(Rap(mctrack->E(),mctrack->Pz())>-4 && Rap(mctrack->E(),mctrack->Pz())<-2.4)) continue;
 494                                 Int_t daug0 = mctrack->GetDaughter(0);
 495                                 Int_t daug1 = mctrack->GetDaughter(1);
 496                                 // daughter1
 497                                 AliAODMCParticle *mcdaug0 = (AliAODMCParticle*) mcArr->At(daug0);
 498                                 Double_t pt0 = mcdaug0->Pt();
 499                                 Double_t mom0 = mcdaug0->P();
 500                                 Double_t theta0 = (180./TMath::Pi())*mcdaug0->Theta();
 501                                 Double_t eta0 = mcdaug0->Eta();
 502                                 Int_t charge0 = (Int_t) mcdaug0->Charge()/3;
 503                                 Double_t vz0 = mcdaug0->Zv();
 504                                 if(!(pt0>0.5) || !(theta0>171. && theta0<178.)) continue;
 505                                 // daughter2
 506                                 AliAODMCParticle *mcdaug1 = (AliAODMCParticle*) mcArr->At(daug1);
 507                                 Double_t pt1 = mcdaug1->Pt();
 508                                 Double_t mom1 = mcdaug1->P();
 509                                 Double_t theta1 = (180./TMath::Pi())*mcdaug1->Theta();
 510                                 Double_t eta1 = mcdaug1->Eta();
 511                                 Int_t charge1 = (Int_t) mcdaug1->Charge()/3;
 512                                 Double_t vz1 = mcdaug1->Zv();
 513                                 if(!(pt1>0.5) || !(theta1>171. && theta1<178.)) continue;
 514
 515                                 if (TMath::Abs(mcdaug0->GetPdgCode())!=13 || TMath::Abs(mcdaug1->GetPdgCode())!=13) continue;
 516
 517                                 Double_t rapmc = Rap(mctrack->E(), mctrack->Pz());
 518                                 Double_t ptmc = mctrack->Pt();
 519                                 Double_t imassmc = mctrack->M();
 520                                 if(!imassmc) continue;
 521
 522                                 containerInput[1] = rapmc ;
 523                                 containerInput[2] = ptmc ;
 524                                 containerInput[3] = imassmc ;
 525                                 containerInput[4] = 1. ;
 526                                 containerInput[5] = pt0 ;
 527                                 containerInput[5] = pt1 ;
 528                                 containerInput[6] = mom0 ;
 529                                 containerInput[6] = mom1 ;
 530                                 containerInput[7] = 1. ;
 531                                 containerInput[8] = 0. ;
 532                                 containerInput[8] = 0. ;
 533                                 containerInput[9] = charge0+charge1 ;
 534                                 containerInput[10] = eta0;
 535                                 containerInput[10] = eta1;
 536                                 containerInput[11] = theta0;
 537                                 containerInput[11] = theta1;
 538                                 containerInput[12] = vz0;
 539                                 containerInput[12] = vz1;
 540                                 containerInput[13] = 0;
 541                                 containerInput[13] = 0;
 542
 543                         fCFManager->GetParticleContainer()->Fill(containerInput,0);     // MC container
 544                                 flag=kTRUE;
 545                         }
 546                 }
 547
 548 if(fReadMCInfo && !flag) return;
 549
 550                 // RECO : AOD part ----------------------------------------------------------------------------------
 551
 552
 553                 Int_t trigfired=-1;
 554                 Int_t trigside=-1;
 555
 556                 Int_t ntrk = fAOD->GetNumberOfTracks();
 557
 558                 if(!fReadMCInfo) {
 559                         fIsPhysSelMB=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMB));
 560                         fIsPhysSelMUON=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMUON));
 561
 562                         if(fDistinguishTrigClass) {
 563                                 TString trigclass = fAOD->GetFiredTriggerClasses();
 564                                 if(trigclass.Contains(fTrigClassMuon)) trigfired = 1;
 565                                 else if(trigclass.Contains(fTrigClassInteraction)) trigfired = 0;
 566                                 for(Int_t i=0;i<3;i++) {
 567                                         if(trigfired==1 && trigclass.Contains(fTrigClassSide[i].Data())) {trigside=i;}
 568                                         if(trigfired==0 && trigclass.Contains(fTrigClassSide[i].Data())) {trigside=i;}
 569                                 }
 570                         }
 571
 572                 if(trigside==1) {       // Beam-Beam event
 573                                 if(trigfired) { // CMUS1-B
 574                                         containerInput[0]=3.5;
 575                                         fnevts->Fill(3.5);
 576                                         if(fIsPhysSelMUON) {
 577                                                 containerInput[0]=4.5;
 578                                                 fnevts->Fill(4.5);      // PhysSel
 579                                         }
 580                                 } else {        // CINT1-B
 581                                         containerInput[0]=1.5;
 582                                         fnevts->Fill(1.5);
 583                                         if(fIsPhysSelMB) {
 584                                                 containerInput[0]=2.5;
 585                                                 fnevts->Fill(2.5);      // PhysSel
 586                                         }
 587                                 }
 588                         }
 589                 }
 590                 else trigside = 0;      // for MC
 591
 592
 593                 for(Int_t j=0; j<ntrk; j++) {
 594                         for(Int_t jj=j+1; jj<ntrk; jj++) {
 595                                 AliAODTrack *mu1 = fAOD->GetTrack(j);
 596                                 if(!mu1->IsMuonTrack() || !(mu1->Y()>-4 && mu1->Y()<-2.5)) continue;
 597                                 AliAODTrack *mu2 = fAOD->GetTrack(jj);
 598                                 if(!mu2->IsMuonTrack() || !(mu2->Y()>-4 && mu2->Y()<-2.5)) continue;
 599
 600                 Double_t trigCondition=0;
 601                 if (mu1->GetMatchTrigger()>=mu2->GetMatchTrigger()) trigCondition = mu1->GetMatchTrigger()+0.1*mu2->GetMatchTrigger();
 602                                 else trigCondition = mu2->GetMatchTrigger()+0.1*mu1->GetMatchTrigger();
 603
 604                                 // mu1
 605                                 Double_t zr1 = mu1->Charge();
 606                                 Double_t pxr1 = mu1->Px();
 607                                 Double_t pyr1 = mu1->Py();
 608                                 Double_t pzr1 = mu1->Pz();
 609                                 Double_t ptr1 = TMath::Sqrt(pxr1*pxr1+pyr1*pyr1);
 610                                 Double_t er1 = mu1->E();
 611                                 Double_t momr1 = mu1->P();
 612                                 Double_t vzr1 = mu1->Zv();
 613                                 Double_t dcar1 = mu1->DCA();
 614                                 Double_t rabsr1 = mu1->GetRAtAbsorberEnd();
 615                                 Double_t thetaabsr1 = 180.*(1.-TMath::ATan(rabsr1/505.)/TMath::Pi());   // [deg]
 616                                 Double_t thetaabsr12 = (thetaabsr1*TMath::Pi()/180.)/2.;
 617                                 Double_t tanthetaabsr1 = TMath::Tan(thetaabsr12);
 618                                 Double_t etar1 = -TMath::Log(tanthetaabsr1);
 619
 620                                 // mu2
 621                                 Double_t zr2 = mu2->Charge();
 622                                 Double_t pxr2 = mu2->Px();
 623                                 Double_t pyr2 = mu2->Py();
 624                                 Double_t pzr2 = mu2->Pz();
 625                                 Double_t ptr2 = TMath::Sqrt(pxr2*pxr2+pyr2*pyr2);
 626                                 Double_t er2 = mu2->E();
 627                                 Double_t momr2 = mu2->P();
 628                                 Double_t vzr2 = mu2->Zv();
 629                                 Double_t dcar2 = mu2->DCA();
 630                                 Double_t rabsr2 = mu1->GetRAtAbsorberEnd();
 631                                 Double_t thetaabsr2 = 180.*(1.-TMath::ATan(rabsr2/505.)/TMath::Pi());   // [deg]
 632                                 Double_t thetaabsr22 = (thetaabsr2*TMath::Pi()/180.)/2.;
 633                                 Double_t tanthetaabsr2 = TMath::Tan(thetaabsr22);
 634                                 Double_t etar2 = -TMath::Log(tanthetaabsr2);
 635
 636                                 // for MC
 637                                 if(fReadMCInfo) {
 638                                         //mu1
 639                                         rabsr1 = 0;
 640                                         Double_t thetar1 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(pxr1*pxr1+pyr1*pyr1),pzr1);
 641                                         thetaabsr1 = thetar1;
 642                                         etar1 = mu1->Eta();
 643                                         //mu2
 644                                         rabsr2 = 0;
 645                                         Double_t thetar2 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(pxr2*pxr2+pyr2*pyr2),pzr2);
 646                                         thetaabsr2 = thetar2;
 647                                         etar2 = mu2->Eta();
 648                                 }
 649
 650                                 if(TMath::Abs(etar1) > 8 || TMath::Abs(etar2) > 8) continue;
 651
 652                                 // dimuon
 653                                 Double_t ptrec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2));
 654                                 Double_t raprec= Rap((er1+er2),(pzr1+pzr2));
 655                                 Double_t imassrec = Imass(er1,pxr1,pyr1,pzr1,er2,pxr2,pyr2,pzr2);
 656                                 if(!imassrec) continue;
 657
 658                                 containerInput[1] = raprec ;
 659                                 containerInput[2] = ptrec ;
 660                                 containerInput[3] = imassrec ;
 661                                 containerInput[4] = trigCondition+0.05 ;
 662                                 containerInput[5] = ptr1 ;
 663                                 containerInput[5] = ptr2 ;
 664                                 containerInput[6] = momr1;
 665                                 containerInput[6] = momr2;
 666                                 containerInput[7] = trigside;
 667                                 containerInput[8] = rabsr1;
 668                                 containerInput[8] = rabsr2;
 669                                 containerInput[9] = zr1 + zr2;
 670                                 containerInput[10] = etar1;
 671                                 containerInput[10] = etar2;
 672                                 containerInput[11] = thetaabsr1;
 673                                 containerInput[11] = thetaabsr2;
 674                                 containerInput[12] = vzr1;
 675                                 containerInput[12] = vzr2;
 676                                 containerInput[13] = dcar1;
 677                                 containerInput[13] = dcar2;
 678
 679                         if (trigfired==1 && trigside==1) fCFManager->GetParticleContainer()->Fill(containerInput,4);
 680                         if (trigfired==0 && trigside==1) fCFManager->GetParticleContainer()->Fill(containerInput,1);
 681                                 if(fReadMCInfo) fCFManager->GetParticleContainer()->Fill(containerInput,1); // Rec container
 682
 683                         } // second mu loop
 684                 } // first mu loop
 685         } // end AOD-based ANALYSIS
 686
 687         // end User analysis loop
 688 }
 689 //________________________________________________________________________
 690 Double_t AliCFMuonResUpsilon::Imass(Float_t e1, Float_t px1, Float_t py1, Float_t pz1,
 691                                    Float_t e2, Float_t px2, Float_t py2, Float_t pz2) const
 692 {
 693 // invariant mass calculation
 694     Float_t imass = (e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+(py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2));
 695                 if(imass<0) return 0;
 696     Float_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+(py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2)));
 697     return imassrec;
 698 }
 699 //________________________________________________________________________
 700 Double_t AliCFMuonResUpsilon::Rap(Float_t e, Float_t pz) const
 701 {
 702 // calculate rapidity
 703     Float_t rap;
 704     if(e!=pz){
 705         rap = 0.5*TMath::Log((e+pz)/(e-pz));
 706         return rap;
 707     }
 708     else{
 709         rap = -200;
 710         return rap;
 711     }
 712 }
 713 //________________________________________________________________________
 714 Double_t AliCFMuonResUpsilon::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 715 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2,
 716 Double_t Energy)
 717 {
 718         // CS angle
 719   TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM. frame
 720   TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
 721   TVector3 beta,zaxisCS;
 722   Double_t mp=0.93827231;
 723   //
 724   // --- Fill the Lorentz vector for projectile and target in the CM frame
 725   //
 726   pProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
 727   pTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
 728   //
 729   // --- Get the muons parameters in the CM frame
 730   //
 731   pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 732   pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 733   //
 734   // --- Obtain the dimuon parameters in the CM frame
 735   //
 736   pDimuCM=pMu1CM+pMu2CM;
 737   //
 738   // --- Translate the dimuon parameters in the dimuon rest frame
 739   //
 740   beta=(-1./pDimuCM.E())*pDimuCM.Vect();
 741   pMu1Dimu=pMu1CM;
 742   pMu2Dimu=pMu2CM;
 743   pProjDimu=pProjCM;
 744   pTargDimu=pTargCM;
 745   pMu1Dimu.Boost(beta);
 746   pMu2Dimu.Boost(beta);
 747   pProjDimu.Boost(beta);
 748   pTargDimu.Boost(beta);
 749   //
 750   // --- Determine the z axis for the CS angle
 751   //
 752   zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
 753   //
 754   // --- Determine the CS angle (angle between mu+ and the z axis defined above)
 755   Double_t cost;
 756
 757   if(charge1>0) cost = zaxisCS.Dot((pMu1Dimu.Vect()).Unit());
 758   else cost = zaxisCS.Dot((pMu2Dimu.Vect()).Unit());
 759   return cost;
 760 }
 761 //________________________________________________________________________
 762
 763 //________________________________________________________________________
 764 Double_t AliCFMuonResUpsilon::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 765 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
 766 {
 767         // Helicity
 768   TLorentzVector pMu1CM, pMu2CM, pDimuCM; // In the CM frame
 769   TLorentzVector pMu1Dimu, pMu2Dimu; // In the dimuon rest frame
 770   TVector3 beta,zaxisCS;
 771   //
 772   // --- Get the muons parameters in the CM frame
 773   //
 774   pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 775   pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 776   //
 777   // --- Obtain the dimuon parameters in the CM frame
 778   //
 779   pDimuCM=pMu1CM+pMu2CM;
 780   //
 781   // --- Translate the muon parameters in the dimuon rest frame
 782   //
 783   beta=(-1./pDimuCM.E())*pDimuCM.Vect();
 784   pMu1Dimu=pMu1CM;
 785   pMu2Dimu=pMu2CM;
 786   pMu1Dimu.Boost(beta);
 787   pMu2Dimu.Boost(beta);
 788   //
 789   // --- Determine the z axis for the calculation of the polarization angle (i.e. the direction of the dimuon in the CM system)
 790   //
 791   TVector3 zaxis;
 792   zaxis=(pDimuCM.Vect()).Unit();
 793
 794   // --- Calculation of the polarization angle (angle between mu+ and the z axis defined above)
 795   Double_t cost;
 796   if(charge1>0) cost = zaxis.Dot((pMu1Dimu.Vect()).Unit());
 797   else cost = zaxis.Dot((pMu2Dimu.Vect()).Unit());
 798
 799   return cost;
 800 }
 801 //________________________________________________________________________
 802
 803 //________________________________________________________________________
 804 Double_t AliCFMuonResUpsilon::PhiCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 805 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2,
 806 Double_t Energy)
 807 {
 808         // CS phi
 809    TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM frame
 810    TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
 811    TVector3 beta,yaxisCS, xaxisCS, zaxisCS;
 812    Double_t mp=0.93827231;
 813    //
 814    // --- Fill the Lorentz vector for projectile and target in the CM frame
 815    //
 816    pProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
 817    pTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
 818    //
 819    // --- Get the muons parameters in the CM frame
 820    //
 821    pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 822    pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 823    //
 824    // --- Obtain the dimuon parameters in the CM frame
 825    //
 826    pDimuCM=pMu1CM+pMu2CM;
 827    //
 828    // --- Translate the dimuon parameters in the dimuon rest frame
 829    //
 830    beta=(-1./pDimuCM.E())*pDimuCM.Vect();
 831    pMu1Dimu=pMu1CM;
 832    pMu2Dimu=pMu2CM;
 833    pProjDimu=pProjCM;
 834    pTargDimu=pTargCM;
 835    pMu1Dimu.Boost(beta);
 836    pMu2Dimu.Boost(beta);
 837    pProjDimu.Boost(beta);
 838    pTargDimu.Boost(beta);
 839    //
 840    // --- Determine the z axis for the CS angle
 841    //
 842    zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
 843    yaxisCS=(((pProjDimu.Vect()).Unit()).Cross((pTargDimu.Vect()).Unit())).Unit();
 844    xaxisCS=(yaxisCS.Cross(zaxisCS)).Unit();
 845
 846    Double_t phi;
 847    if(charge1>0) phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxisCS),((pMu1Dimu.Vect()).Dot(xaxisCS)));
 848    else phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxisCS),((pMu2Dimu.Vect()).Dot(xaxisCS)));
 849
 850    return phi;
 851 }
 852 //________________________________________________________________________
 853
 854 //________________________________________________________________________
 855 Double_t AliCFMuonResUpsilon::PhiHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 856 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2, Double_t Energy)
 857 {
 858         // Helicity phi
 859    TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM frame
 860    TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
 861    TVector3 beta,xaxis,yaxis,zaxis;
 862    Double_t mp=0.93827231;
 863
 864    //
 865    // --- Get the muons parameters in the CM frame
 866    //
 867    pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 868    pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 869    //
 870    // --- Obtain the dimuon parameters in the CM frame
 871    //
 872    pDimuCM=pMu1CM+pMu2CM;
 873    //
 874    // --- Translate the muon parameters in the dimuon rest frame
 875    //
 876    zaxis=(pDimuCM.Vect()).Unit();
 877
 878    beta=(-1./pDimuCM.E())*pDimuCM.Vect();
 879
 880    pProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
 881    pTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
 882
 883    pProjDimu=pProjCM;
 884    pTargDimu=pTargCM;
 885
 886    pProjDimu.Boost(beta);
 887    pTargDimu.Boost(beta);
 888
 889    yaxis=((pProjDimu.Vect()).Cross(pTargDimu.Vect())).Unit();
 890    xaxis=(yaxis.Cross(zaxis)).Unit();
 891
 892    pMu1Dimu=pMu1CM;
 893    pMu2Dimu=pMu2CM;
 894    pMu1Dimu.Boost(beta);
 895    pMu2Dimu.Boost(beta);
 896
 897    Double_t phi;
 898    if(charge1>0) phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxis),(pMu1Dimu.Vect()).Dot(xaxis));
 899    else phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxis),(pMu2Dimu.Vect()).Dot(xaxis));
 900
 901    return phi;
 902 }
 903
 904 //________________________________________________________________________
 905 void AliCFMuonResUpsilon::Terminate(Option_t *)
 906 {
 907   // draw result of the Invariant mass MC and ESD
 908
 909 /*
 910                 TH1D *h1 = dynamic_cast<TH1D*>(GetOutputData(1));
 911     AliCFContainer *cont = dynamic_cast<AliCFContainer*> (GetOutputData(2));
 912
 913     TH1D *kmass = cont->ShowProjection(3,0);
 914     TH1D *rmass = cont->ShowProjection(3,1);
 915                 TH1D *mmass = cont->ShowProjection(3,4);
 916
 917     TCanvas *c1 = new TCanvas("AliCFMuonResUpsilon","UPSILON Container",0,0,800,800);
 918     c1->Divide(2,2);
 919     c1->cd(1);
 920     kmass->Draw("HIST");
 921     c1->cd(2);
 922     rmass->Draw("HIST");
 923                 c1->cd(3);
 924                 mmass->Draw("HIST");
 925                 c1->cd(4);
 926                 h1->Draw("HIST");
 927                 */
 928 }
 929 //________________________________________________________________________
 930
 931 #endif