PWG/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 = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(j));
 596                                 if(!mu1) AliFatal("Not a standard AOD");
 597                                 if(!mu1->IsMuonTrack() || !(mu1->Y()>-4 && mu1->Y()<-2.5)) continue;
 598                                 AliAODTrack *mu2 = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(jj));
 599                                 if(!mu2) AliFatal("Not a standard AOD");
 600                                 if(!mu2->IsMuonTrack() || !(mu2->Y()>-4 && mu2->Y()<-2.5)) continue;
 601
 602                 Double_t trigCondition=0;
 603                 if (mu1->GetMatchTrigger()>=mu2->GetMatchTrigger()) trigCondition = mu1->GetMatchTrigger()+0.1*mu2->GetMatchTrigger();
 604                                 else trigCondition = mu2->GetMatchTrigger()+0.1*mu1->GetMatchTrigger();
 605
 606                                 // mu1
 607                                 Double_t zr1 = mu1->Charge();
 608                                 Double_t pxr1 = mu1->Px();
 609                                 Double_t pyr1 = mu1->Py();
 610                                 Double_t pzr1 = mu1->Pz();
 611                                 Double_t ptr1 = TMath::Sqrt(pxr1*pxr1+pyr1*pyr1);
 612                                 Double_t er1 = mu1->E();
 613                                 Double_t momr1 = mu1->P();
 614                                 Double_t vzr1 = mu1->Zv();
 615                                 Double_t dcar1 = mu1->DCA();
 616                                 Double_t rabsr1 = mu1->GetRAtAbsorberEnd();
 617                                 Double_t thetaabsr1 = 180.*(1.-TMath::ATan(rabsr1/505.)/TMath::Pi());   // [deg]
 618                                 Double_t thetaabsr12 = (thetaabsr1*TMath::Pi()/180.)/2.;
 619                                 Double_t tanthetaabsr1 = TMath::Tan(thetaabsr12);
 620                                 Double_t etar1 = -TMath::Log(tanthetaabsr1);
 621
 622                                 // mu2
 623                                 Double_t zr2 = mu2->Charge();
 624                                 Double_t pxr2 = mu2->Px();
 625                                 Double_t pyr2 = mu2->Py();
 626                                 Double_t pzr2 = mu2->Pz();
 627                                 Double_t ptr2 = TMath::Sqrt(pxr2*pxr2+pyr2*pyr2);
 628                                 Double_t er2 = mu2->E();
 629                                 Double_t momr2 = mu2->P();
 630                                 Double_t vzr2 = mu2->Zv();
 631                                 Double_t dcar2 = mu2->DCA();
 632                                 Double_t rabsr2 = mu1->GetRAtAbsorberEnd();
 633                                 Double_t thetaabsr2 = 180.*(1.-TMath::ATan(rabsr2/505.)/TMath::Pi());   // [deg]
 634                                 Double_t thetaabsr22 = (thetaabsr2*TMath::Pi()/180.)/2.;
 635                                 Double_t tanthetaabsr2 = TMath::Tan(thetaabsr22);
 636                                 Double_t etar2 = -TMath::Log(tanthetaabsr2);
 637
 638                                 // for MC
 639                                 if(fReadMCInfo) {
 640                                         //mu1
 641                                         rabsr1 = 0;
 642                                         Double_t thetar1 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(pxr1*pxr1+pyr1*pyr1),pzr1);
 643                                         thetaabsr1 = thetar1;
 644                                         etar1 = mu1->Eta();
 645                                         //mu2
 646                                         rabsr2 = 0;
 647                                         Double_t thetar2 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(pxr2*pxr2+pyr2*pyr2),pzr2);
 648                                         thetaabsr2 = thetar2;
 649                                         etar2 = mu2->Eta();
 650                                 }
 651
 652                                 if(TMath::Abs(etar1) > 8 || TMath::Abs(etar2) > 8) continue;
 653
 654                                 // dimuon
 655                                 Double_t ptrec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2));
 656                                 Double_t raprec= Rap((er1+er2),(pzr1+pzr2));
 657                                 Double_t imassrec = Imass(er1,pxr1,pyr1,pzr1,er2,pxr2,pyr2,pzr2);
 658                                 if(!imassrec) continue;
 659
 660                                 containerInput[1] = raprec ;
 661                                 containerInput[2] = ptrec ;
 662                                 containerInput[3] = imassrec ;
 663                                 containerInput[4] = trigCondition+0.05 ;
 664                                 containerInput[5] = ptr1 ;
 665                                 containerInput[5] = ptr2 ;
 666                                 containerInput[6] = momr1;
 667                                 containerInput[6] = momr2;
 668                                 containerInput[7] = trigside;
 669                                 containerInput[8] = rabsr1;
 670                                 containerInput[8] = rabsr2;
 671                                 containerInput[9] = zr1 + zr2;
 672                                 containerInput[10] = etar1;
 673                                 containerInput[10] = etar2;
 674                                 containerInput[11] = thetaabsr1;
 675                                 containerInput[11] = thetaabsr2;
 676                                 containerInput[12] = vzr1;
 677                                 containerInput[12] = vzr2;
 678                                 containerInput[13] = dcar1;
 679                                 containerInput[13] = dcar2;
 680
 681                         if (trigfired==1 && trigside==1) fCFManager->GetParticleContainer()->Fill(containerInput,4);
 682                         if (trigfired==0 && trigside==1) fCFManager->GetParticleContainer()->Fill(containerInput,1);
 683                                 if(fReadMCInfo) fCFManager->GetParticleContainer()->Fill(containerInput,1); // Rec container
 684
 685                         } // second mu loop
 686                 } // first mu loop
 687         } // end AOD-based ANALYSIS
 688
 689         // end User analysis loop
 690 }
 691 //________________________________________________________________________
 692 Double_t AliCFMuonResUpsilon::Imass(Float_t e1, Float_t px1, Float_t py1, Float_t pz1,
 693                                    Float_t e2, Float_t px2, Float_t py2, Float_t pz2) const
 694 {
 695 // invariant mass calculation
 696     Float_t imass = (e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+(py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2));
 697                 if(imass<0) return 0;
 698     Float_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+(py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2)));
 699     return imassrec;
 700 }
 701 //________________________________________________________________________
 702 Double_t AliCFMuonResUpsilon::Rap(Float_t e, Float_t pz) const
 703 {
 704 // calculate rapidity
 705     Float_t rap;
 706     if(e!=pz){
 707         rap = 0.5*TMath::Log((e+pz)/(e-pz));
 708         return rap;
 709     }
 710     else{
 711         rap = -200;
 712         return rap;
 713     }
 714 }
 715 //________________________________________________________________________
 716 Double_t AliCFMuonResUpsilon::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 717 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2,
 718 Double_t Energy)
 719 {
 720         // CS angle
 721   TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM. frame
 722   TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
 723   TVector3 beta,zaxisCS;
 724   Double_t mp=0.93827231;
 725   //
 726   // --- Fill the Lorentz vector for projectile and target in the CM frame
 727   //
 728   pProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
 729   pTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
 730   //
 731   // --- Get the muons parameters in the CM frame
 732   //
 733   pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 734   pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 735   //
 736   // --- Obtain the dimuon parameters in the CM frame
 737   //
 738   pDimuCM=pMu1CM+pMu2CM;
 739   //
 740   // --- Translate the dimuon parameters in the dimuon rest frame
 741   //
 742   beta=(-1./pDimuCM.E())*pDimuCM.Vect();
 743   pMu1Dimu=pMu1CM;
 744   pMu2Dimu=pMu2CM;
 745   pProjDimu=pProjCM;
 746   pTargDimu=pTargCM;
 747   pMu1Dimu.Boost(beta);
 748   pMu2Dimu.Boost(beta);
 749   pProjDimu.Boost(beta);
 750   pTargDimu.Boost(beta);
 751   //
 752   // --- Determine the z axis for the CS angle
 753   //
 754   zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
 755   //
 756   // --- Determine the CS angle (angle between mu+ and the z axis defined above)
 757   Double_t cost;
 758
 759   if(charge1>0) cost = zaxisCS.Dot((pMu1Dimu.Vect()).Unit());
 760   else cost = zaxisCS.Dot((pMu2Dimu.Vect()).Unit());
 761   return cost;
 762 }
 763 //________________________________________________________________________
 764
 765 //________________________________________________________________________
 766 Double_t AliCFMuonResUpsilon::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 767 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
 768 {
 769         // Helicity
 770   TLorentzVector pMu1CM, pMu2CM, pDimuCM; // In the CM frame
 771   TLorentzVector pMu1Dimu, pMu2Dimu; // In the dimuon rest frame
 772   TVector3 beta,zaxisCS;
 773   //
 774   // --- Get the muons parameters in the CM frame
 775   //
 776   pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 777   pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 778   //
 779   // --- Obtain the dimuon parameters in the CM frame
 780   //
 781   pDimuCM=pMu1CM+pMu2CM;
 782   //
 783   // --- Translate the muon parameters in the dimuon rest frame
 784   //
 785   beta=(-1./pDimuCM.E())*pDimuCM.Vect();
 786   pMu1Dimu=pMu1CM;
 787   pMu2Dimu=pMu2CM;
 788   pMu1Dimu.Boost(beta);
 789   pMu2Dimu.Boost(beta);
 790   //
 791   // --- Determine the z axis for the calculation of the polarization angle (i.e. the direction of the dimuon in the CM system)
 792   //
 793   TVector3 zaxis;
 794   zaxis=(pDimuCM.Vect()).Unit();
 795
 796   // --- Calculation of the polarization angle (angle between mu+ and the z axis defined above)
 797   Double_t cost;
 798   if(charge1>0) cost = zaxis.Dot((pMu1Dimu.Vect()).Unit());
 799   else cost = zaxis.Dot((pMu2Dimu.Vect()).Unit());
 800
 801   return cost;
 802 }
 803 //________________________________________________________________________
 804
 805 //________________________________________________________________________
 806 Double_t AliCFMuonResUpsilon::PhiCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 807 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2,
 808 Double_t Energy)
 809 {
 810         // CS phi
 811    TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM frame
 812    TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
 813    TVector3 beta,yaxisCS, xaxisCS, zaxisCS;
 814    Double_t mp=0.93827231;
 815    //
 816    // --- Fill the Lorentz vector for projectile and target in the CM frame
 817    //
 818    pProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
 819    pTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
 820    //
 821    // --- Get the muons parameters in the CM frame
 822    //
 823    pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 824    pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 825    //
 826    // --- Obtain the dimuon parameters in the CM frame
 827    //
 828    pDimuCM=pMu1CM+pMu2CM;
 829    //
 830    // --- Translate the dimuon parameters in the dimuon rest frame
 831    //
 832    beta=(-1./pDimuCM.E())*pDimuCM.Vect();
 833    pMu1Dimu=pMu1CM;
 834    pMu2Dimu=pMu2CM;
 835    pProjDimu=pProjCM;
 836    pTargDimu=pTargCM;
 837    pMu1Dimu.Boost(beta);
 838    pMu2Dimu.Boost(beta);
 839    pProjDimu.Boost(beta);
 840    pTargDimu.Boost(beta);
 841    //
 842    // --- Determine the z axis for the CS angle
 843    //
 844    zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
 845    yaxisCS=(((pProjDimu.Vect()).Unit()).Cross((pTargDimu.Vect()).Unit())).Unit();
 846    xaxisCS=(yaxisCS.Cross(zaxisCS)).Unit();
 847
 848    Double_t phi;
 849    if(charge1>0) phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxisCS),((pMu1Dimu.Vect()).Dot(xaxisCS)));
 850    else phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxisCS),((pMu2Dimu.Vect()).Dot(xaxisCS)));
 851
 852    return phi;
 853 }
 854 //________________________________________________________________________
 855
 856 //________________________________________________________________________
 857 Double_t AliCFMuonResUpsilon::PhiHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 858 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2, Double_t Energy)
 859 {
 860         // Helicity phi
 861    TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM frame
 862    TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
 863    TVector3 beta,xaxis,yaxis,zaxis;
 864    Double_t mp=0.93827231;
 865
 866    //
 867    // --- Get the muons parameters in the CM frame
 868    //
 869    pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 870    pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 871    //
 872    // --- Obtain the dimuon parameters in the CM frame
 873    //
 874    pDimuCM=pMu1CM+pMu2CM;
 875    //
 876    // --- Translate the muon parameters in the dimuon rest frame
 877    //
 878    zaxis=(pDimuCM.Vect()).Unit();
 879
 880    beta=(-1./pDimuCM.E())*pDimuCM.Vect();
 881
 882    pProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
 883    pTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
 884
 885    pProjDimu=pProjCM;
 886    pTargDimu=pTargCM;
 887
 888    pProjDimu.Boost(beta);
 889    pTargDimu.Boost(beta);
 890
 891    yaxis=((pProjDimu.Vect()).Cross(pTargDimu.Vect())).Unit();
 892    xaxis=(yaxis.Cross(zaxis)).Unit();
 893
 894    pMu1Dimu=pMu1CM;
 895    pMu2Dimu=pMu2CM;
 896    pMu1Dimu.Boost(beta);
 897    pMu2Dimu.Boost(beta);
 898
 899    Double_t phi;
 900    if(charge1>0) phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxis),(pMu1Dimu.Vect()).Dot(xaxis));
 901    else phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxis),(pMu2Dimu.Vect()).Dot(xaxis));
 902
 903    return phi;
 904 }
 905
 906 //________________________________________________________________________
 907 void AliCFMuonResUpsilon::Terminate(Option_t *)
 908 {
 909   // draw result of the Invariant mass MC and ESD
 910
 911 /*
 912                 TH1D *h1 = dynamic_cast<TH1D*>(GetOutputData(1));
 913     AliCFContainer *cont = dynamic_cast<AliCFContainer*> (GetOutputData(2));
 914
 915     TH1D *kmass = cont->ShowProjection(3,0);
 916     TH1D *rmass = cont->ShowProjection(3,1);
 917                 TH1D *mmass = cont->ShowProjection(3,4);
 918
 919     TCanvas *c1 = new TCanvas("AliCFMuonResUpsilon","UPSILON Container",0,0,800,800);
 920     c1->Divide(2,2);
 921     c1->cd(1);
 922     kmass->Draw("HIST");
 923     c1->cd(2);
 924     rmass->Draw("HIST");
 925                 c1->cd(3);
 926                 mmass->Draw("HIST");
 927                 c1->cd(4);
 928                 h1->Draw("HIST");
 929                 */
 930 }
 931 //________________________________________________________________________
 932
 933 #endif