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         hnevts(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         hnevts(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         hnevts  = new TH1D("hnevts","hnevts",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                 hnevts = c.hnevts;
 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         hnevts(c.hnevts),
 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 (hnevts)                                                             delete hnevts ;
 182   if (fQAHistList)                                      {fQAHistList->Clear(); delete fQAHistList;}
 183 }
 184 //___________________________________________________________________________
 185 void AliCFMuonResUpsilon::UserCreateOutputObjects() {
 186         // UserCreateOutputObjects
 187
 188         //TH1D *hnevts  = new TH1D("hnevts","hnevts",5,0,5);                                                    // nevent CINT1B,CMUS1B + PhysSel
 189
 190   PostData(1,hnevts) ;
 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   hnevts->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                                         hnevts->Fill(3.5);
 353                                         if(fIsPhysSelMUON) {
 354                                                 containerInput[0]=4.5;
 355                                                 hnevts->Fill(4.5);      // PhysSel
 356                                         }
 357                                 } else {        // CINT1-B
 358                                         containerInput[0]=1.5;
 359                                         hnevts->Fill(1.5);
 360                                         if(fIsPhysSelMB) {
 361                                                 containerInput[0]=2.5;
 362                                                 hnevts->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 theta_absr1 = 180.*(1.-TMath::ATan(rabsr1/505.)/TMath::Pi());  // [deg]
 390                                 Double_t theta_absr1_2 = (theta_absr1*TMath::Pi()/180.)/2.;
 391                                 Double_t tan_theta_absr1 = TMath::Tan(theta_absr1_2);
 392                                 Double_t etar1 = -TMath::Log(tan_theta_absr1);
 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 theta_absr2 = 180.*(1.-TMath::ATan(rabsr2/505.)/TMath::Pi());  // [deg]
 408                                 Double_t theta_absr2_2 = (theta_absr2*TMath::Pi()/180.)/2.;
 409                                 Double_t tan_theta_absr2 = TMath::Tan(theta_absr2_2);
 410                                 Double_t etar2 = -TMath::Log(tan_theta_absr2);
 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                                         theta_absr1 = 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                                         theta_absr2 = 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] = theta_absr1;
 451                                 containerInput[11] = theta_absr2;
 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
 483                         Int_t npart=mcArr->GetEntries();
 484
 485                         for(Int_t i=0; i<npart; i++) {
 486                                 AliAODMCParticle *mctrack = (AliAODMCParticle*) mcArr->At(i);
 487                                 // select resonance
 488                                 if(mctrack->GetPdgCode()!=fPDG) continue;
 489                                 // cuts on resonance
 490                                 if(!(mctrack->Pt()>0 && mctrack->Pt()<1000) || !(Rap(mctrack->E(),mctrack->Pz())>-4 && Rap(mctrack->E(),mctrack->Pz())<-2.4)) continue;
 491                                 Int_t daug0 = mctrack->GetDaughter(0);
 492                                 Int_t daug1 = mctrack->GetDaughter(1);
 493                                 // daughter1
 494                                 AliAODMCParticle *mcdaug0 = (AliAODMCParticle*) mcArr->At(daug0);
 495                                 Double_t pt0 = mcdaug0->Pt();
 496                                 Double_t mom0 = mcdaug0->P();
 497                                 Double_t theta0 = (180./TMath::Pi())*mcdaug0->Theta();
 498                                 Double_t eta0 = mcdaug0->Eta();
 499                                 Int_t charge0 = (Int_t) mcdaug0->Charge()/3;
 500                                 Double_t vz0 = mcdaug0->Zv();
 501                                 if(!(pt0>0.5) || !(theta0>171. && theta0<178.)) continue;
 502                                 // daughter2
 503                                 AliAODMCParticle *mcdaug1 = (AliAODMCParticle*) mcArr->At(daug1);
 504                                 Double_t pt1 = mcdaug1->Pt();
 505                                 Double_t mom1 = mcdaug1->P();
 506                                 Double_t theta1 = (180./TMath::Pi())*mcdaug1->Theta();
 507                                 Double_t eta1 = mcdaug1->Eta();
 508                                 Int_t charge1 = (Int_t) mcdaug1->Charge()/3;
 509                                 Double_t vz1 = mcdaug1->Zv();
 510                                 if(!(pt1>0.5) || !(theta1>171. && theta1<178.)) continue;
 511
 512                                 if (TMath::Abs(mcdaug0->GetPdgCode())!=13 || TMath::Abs(mcdaug1->GetPdgCode())!=13) continue;
 513
 514                                 Double_t rapmc = Rap(mctrack->E(), mctrack->Pz());
 515                                 Double_t ptmc = mctrack->Pt();
 516                                 Double_t imassmc = mctrack->M();
 517                                 if(!imassmc) continue;
 518
 519                                 containerInput[1] = rapmc ;
 520                                 containerInput[2] = ptmc ;
 521                                 containerInput[3] = imassmc ;
 522                                 containerInput[4] = 1. ;
 523                                 containerInput[5] = pt0 ;
 524                                 containerInput[5] = pt1 ;
 525                                 containerInput[6] = mom0 ;
 526                                 containerInput[6] = mom1 ;
 527                                 containerInput[7] = 1. ;
 528                                 containerInput[8] = 0. ;
 529                                 containerInput[8] = 0. ;
 530                                 containerInput[9] = charge0+charge1 ;
 531                                 containerInput[10] = eta0;
 532                                 containerInput[10] = eta1;
 533                                 containerInput[11] = theta0;
 534                                 containerInput[11] = theta1;
 535                                 containerInput[12] = vz0;
 536                                 containerInput[12] = vz1;
 537                                 containerInput[13] = 0;
 538                                 containerInput[13] = 0;
 539
 540                         fCFManager->GetParticleContainer()->Fill(containerInput,0);     // MC container
 541                                 flag=kTRUE;
 542                         }
 543                 }
 544
 545 if(fReadMCInfo && !flag) return;
 546
 547                 // RECO : AOD part ----------------------------------------------------------------------------------
 548
 549
 550                 Int_t trigfired=-1;
 551                 Int_t trigside=-1;
 552
 553                 Int_t ntrk = fAOD->GetNumberOfTracks();
 554
 555                 if(!fReadMCInfo) {
 556                         fIsPhysSelMB=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMB));
 557                         fIsPhysSelMUON=(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMUON));
 558
 559                         if(fDistinguishTrigClass) {
 560                                 TString trigclass = fAOD->GetFiredTriggerClasses();
 561                                 if(trigclass.Contains(fTrigClassMuon)) trigfired = 1;
 562                                 else if(trigclass.Contains(fTrigClassInteraction)) trigfired = 0;
 563                                 for(Int_t i=0;i<3;i++) {
 564                                         if(trigfired==1 && trigclass.Contains(fTrigClassSide[i].Data())) {trigside=i;}
 565                                         if(trigfired==0 && trigclass.Contains(fTrigClassSide[i].Data())) {trigside=i;}
 566                                 }
 567                         }
 568
 569                 if(trigside==1) {       // Beam-Beam event
 570                                 if(trigfired) { // CMUS1-B
 571                                         containerInput[0]=3.5;
 572                                         hnevts->Fill(3.5);
 573                                         if(fIsPhysSelMUON) {
 574                                                 containerInput[0]=4.5;
 575                                                 hnevts->Fill(4.5);      // PhysSel
 576                                         }
 577                                 } else {        // CINT1-B
 578                                         containerInput[0]=1.5;
 579                                         hnevts->Fill(1.5);
 580                                         if(fIsPhysSelMB) {
 581                                                 containerInput[0]=2.5;
 582                                                 hnevts->Fill(2.5);      // PhysSel
 583                                         }
 584                                 }
 585                         }
 586                 }
 587                 else trigside = 0;      // for MC
 588
 589
 590                 for(Int_t j=0; j<ntrk; j++) {
 591                         for(Int_t jj=j+1; jj<ntrk; jj++) {
 592                                 AliAODTrack *mu1 = fAOD->GetTrack(j);
 593                                 if(!mu1->IsMuonTrack() || !(mu1->Y()>-4 && mu1->Y()<-2.5)) continue;
 594                                 AliAODTrack *mu2 = fAOD->GetTrack(jj);
 595                                 if(!mu2->IsMuonTrack() || !(mu2->Y()>-4 && mu2->Y()<-2.5)) continue;
 596
 597                 Double_t trigCondition=0;
 598                 if (mu1->GetMatchTrigger()>=mu2->GetMatchTrigger()) trigCondition = mu1->GetMatchTrigger()+0.1*mu2->GetMatchTrigger();
 599                                 else trigCondition = mu2->GetMatchTrigger()+0.1*mu1->GetMatchTrigger();
 600
 601                                 // mu1
 602                                 Double_t zr1 = mu1->Charge();
 603                                 Double_t pxr1 = mu1->Px();
 604                                 Double_t pyr1 = mu1->Py();
 605                                 Double_t pzr1 = mu1->Pz();
 606                                 Double_t ptr1 = TMath::Sqrt(pxr1*pxr1+pyr1*pyr1);
 607                                 Double_t er1 = mu1->E();
 608                                 Double_t momr1 = mu1->P();
 609                                 Double_t vzr1 = mu1->Zv();
 610                                 Double_t dcar1 = mu1->DCA();
 611                                 Double_t rabsr1 = mu1->GetRAtAbsorberEnd();
 612                                 Double_t theta_absr1 = 180.*(1.-TMath::ATan(rabsr1/505.)/TMath::Pi());  // [deg]
 613                                 Double_t theta_absr1_2 = (theta_absr1*TMath::Pi()/180.)/2.;
 614                                 Double_t tan_theta_absr1 = TMath::Tan(theta_absr1_2);
 615                                 Double_t etar1 = -TMath::Log(tan_theta_absr1);
 616
 617                                 // mu2
 618                                 Double_t zr2 = mu2->Charge();
 619                                 Double_t pxr2 = mu2->Px();
 620                                 Double_t pyr2 = mu2->Py();
 621                                 Double_t pzr2 = mu2->Pz();
 622                                 Double_t ptr2 = TMath::Sqrt(pxr2*pxr2+pyr2*pyr2);
 623                                 Double_t er2 = mu2->E();
 624                                 Double_t momr2 = mu2->P();
 625                                 Double_t vzr2 = mu2->Zv();
 626                                 Double_t dcar2 = mu2->DCA();
 627                                 Double_t rabsr2 = mu1->GetRAtAbsorberEnd();
 628                                 Double_t theta_absr2 = 180.*(1.-TMath::ATan(rabsr2/505.)/TMath::Pi());  // [deg]
 629                                 Double_t theta_absr2_2 = (theta_absr2*TMath::Pi()/180.)/2.;
 630                                 Double_t tan_theta_absr2 = TMath::Tan(theta_absr2_2);
 631                                 Double_t etar2 = -TMath::Log(tan_theta_absr2);
 632
 633                                 // for MC
 634                                 if(fReadMCInfo) {
 635                                         //mu1
 636                                         rabsr1 = 0;
 637                                         Double_t thetar1 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(pxr1*pxr1+pyr1*pyr1),pzr1);
 638                                         theta_absr1 = thetar1;
 639                                         etar1 = mu1->Eta();
 640                                         //mu2
 641                                         rabsr2 = 0;
 642                                         Double_t thetar2 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(pxr2*pxr2+pyr2*pyr2),pzr2);
 643                                         theta_absr2 = thetar2;
 644                                         etar2 = mu2->Eta();
 645                                 }
 646
 647                                 if(TMath::Abs(etar1) > 8 || TMath::Abs(etar2) > 8) continue;
 648
 649                                 // dimuon
 650                                 Double_t ptrec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2));
 651                                 Double_t raprec= Rap((er1+er2),(pzr1+pzr2));
 652                                 Double_t imassrec = Imass(er1,pxr1,pyr1,pzr1,er2,pxr2,pyr2,pzr2);
 653                                 if(!imassrec) continue;
 654
 655                                 containerInput[1] = raprec ;
 656                                 containerInput[2] = ptrec ;
 657                                 containerInput[3] = imassrec ;
 658                                 containerInput[4] = trigCondition+0.05 ;
 659                                 containerInput[5] = ptr1 ;
 660                                 containerInput[5] = ptr2 ;
 661                                 containerInput[6] = momr1;
 662                                 containerInput[6] = momr2;
 663                                 containerInput[7] = trigside;
 664                                 containerInput[8] = rabsr1;
 665                                 containerInput[8] = rabsr2;
 666                                 containerInput[9] = zr1 + zr2;
 667                                 containerInput[10] = etar1;
 668                                 containerInput[10] = etar2;
 669                                 containerInput[11] = theta_absr1;
 670                                 containerInput[11] = theta_absr2;
 671                                 containerInput[12] = vzr1;
 672                                 containerInput[12] = vzr2;
 673                                 containerInput[13] = dcar1;
 674                                 containerInput[13] = dcar2;
 675
 676                         if (trigfired==1 && trigside==1) fCFManager->GetParticleContainer()->Fill(containerInput,4);
 677                         if (trigfired==0 && trigside==1) fCFManager->GetParticleContainer()->Fill(containerInput,1);
 678                                 if(fReadMCInfo) fCFManager->GetParticleContainer()->Fill(containerInput,1); // Rec container
 679
 680                         } // second mu loop
 681                 } // first mu loop
 682         } // end AOD-based ANALYSIS
 683
 684         // end User analysis loop
 685 }
 686 //________________________________________________________________________
 687 Double_t AliCFMuonResUpsilon::Imass(Float_t e1, Float_t px1, Float_t py1, Float_t pz1,
 688                                    Float_t e2, Float_t px2, Float_t py2, Float_t pz2) const
 689 {
 690 // invariant mass calculation
 691     Float_t imass = (e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+(py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2));
 692                 if(imass<0) return 0;
 693     Float_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+(py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2)));
 694     return imassrec;
 695 }
 696 //________________________________________________________________________
 697 Double_t AliCFMuonResUpsilon::Rap(Float_t e, Float_t pz) const
 698 {
 699 // calculate rapidity
 700     Float_t rap;
 701     if(e!=pz){
 702         rap = 0.5*TMath::Log((e+pz)/(e-pz));
 703         return rap;
 704     }
 705     else{
 706         rap = -200;
 707         return rap;
 708     }
 709 }
 710 //________________________________________________________________________
 711 Double_t AliCFMuonResUpsilon::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 712 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2,
 713 Double_t Energy)
 714 {
 715   TLorentzVector PMu1CM, PMu2CM, PProjCM, PTargCM, PDimuCM; // In the CM. frame
 716   TLorentzVector PMu1Dimu, PMu2Dimu, PProjDimu, PTargDimu; // In the dimuon rest frame
 717   TVector3 beta,zaxisCS;
 718   Double_t mp=0.93827231;
 719   //
 720   // --- Fill the Lorentz vector for projectile and target in the CM frame
 721   //
 722   PProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
 723   PTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
 724   //
 725   // --- Get the muons parameters in the CM frame
 726   //
 727   PMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 728   PMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 729   //
 730   // --- Obtain the dimuon parameters in the CM frame
 731   //
 732   PDimuCM=PMu1CM+PMu2CM;
 733   //
 734   // --- Translate the dimuon parameters in the dimuon rest frame
 735   //
 736   beta=(-1./PDimuCM.E())*PDimuCM.Vect();
 737   PMu1Dimu=PMu1CM;
 738   PMu2Dimu=PMu2CM;
 739   PProjDimu=PProjCM;
 740   PTargDimu=PTargCM;
 741   PMu1Dimu.Boost(beta);
 742   PMu2Dimu.Boost(beta);
 743   PProjDimu.Boost(beta);
 744   PTargDimu.Boost(beta);
 745   //
 746   // --- Determine the z axis for the CS angle
 747   //
 748   zaxisCS=(((PProjDimu.Vect()).Unit())-((PTargDimu.Vect()).Unit())).Unit();
 749   //
 750   // --- Determine the CS angle (angle between mu+ and the z axis defined above)
 751   Double_t cost;
 752
 753   if(charge1>0) cost = zaxisCS.Dot((PMu1Dimu.Vect()).Unit());
 754   else cost = zaxisCS.Dot((PMu2Dimu.Vect()).Unit());
 755   return cost;
 756 }
 757 //________________________________________________________________________
 758
 759 //________________________________________________________________________
 760 Double_t AliCFMuonResUpsilon::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 761 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
 762 {
 763   TLorentzVector PMu1CM, PMu2CM, PDimuCM; // In the CM frame
 764   TLorentzVector PMu1Dimu, PMu2Dimu; // In the dimuon rest frame
 765   TVector3 beta,zaxisCS;
 766   //
 767   // --- Get the muons parameters in the CM frame
 768   //
 769   PMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 770   PMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 771   //
 772   // --- Obtain the dimuon parameters in the CM frame
 773   //
 774   PDimuCM=PMu1CM+PMu2CM;
 775   //
 776   // --- Translate the muon parameters in the dimuon rest frame
 777   //
 778   beta=(-1./PDimuCM.E())*PDimuCM.Vect();
 779   PMu1Dimu=PMu1CM;
 780   PMu2Dimu=PMu2CM;
 781   PMu1Dimu.Boost(beta);
 782   PMu2Dimu.Boost(beta);
 783   //
 784   // --- Determine the z axis for the calculation of the polarization angle (i.e. the direction of the dimuon in the CM system)
 785   //
 786   TVector3 zaxis;
 787   zaxis=(PDimuCM.Vect()).Unit();
 788
 789   // --- Calculation of the polarization angle (angle between mu+ and the z axis defined above)
 790   Double_t cost;
 791   if(charge1>0) cost = zaxis.Dot((PMu1Dimu.Vect()).Unit());
 792   else cost = zaxis.Dot((PMu2Dimu.Vect()).Unit());
 793
 794   return cost;
 795 }
 796 //________________________________________________________________________
 797
 798 //________________________________________________________________________
 799 Double_t AliCFMuonResUpsilon::PhiCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 800 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2,
 801 Double_t Energy)
 802 {
 803    TLorentzVector PMu1CM, PMu2CM, PProjCM, PTargCM, PDimuCM; // In the CM frame
 804    TLorentzVector PMu1Dimu, PMu2Dimu, PProjDimu, PTargDimu; // In the dimuon rest frame
 805    TVector3 beta,yaxisCS, xaxisCS, zaxisCS;
 806    Double_t mp=0.93827231;
 807    //
 808    // --- Fill the Lorentz vector for projectile and target in the CM frame
 809    //
 810    PProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
 811    PTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
 812    //
 813    // --- Get the muons parameters in the CM frame
 814    //
 815    PMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 816    PMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 817    //
 818    // --- Obtain the dimuon parameters in the CM frame
 819    //
 820    PDimuCM=PMu1CM+PMu2CM;
 821    //
 822    // --- Translate the dimuon parameters in the dimuon rest frame
 823    //
 824    beta=(-1./PDimuCM.E())*PDimuCM.Vect();
 825    PMu1Dimu=PMu1CM;
 826    PMu2Dimu=PMu2CM;
 827    PProjDimu=PProjCM;
 828    PTargDimu=PTargCM;
 829    PMu1Dimu.Boost(beta);
 830    PMu2Dimu.Boost(beta);
 831    PProjDimu.Boost(beta);
 832    PTargDimu.Boost(beta);
 833    //
 834    // --- Determine the z axis for the CS angle
 835    //
 836    zaxisCS=(((PProjDimu.Vect()).Unit())-((PTargDimu.Vect()).Unit())).Unit();
 837    yaxisCS=(((PProjDimu.Vect()).Unit()).Cross((PTargDimu.Vect()).Unit())).Unit();
 838    xaxisCS=(yaxisCS.Cross(zaxisCS)).Unit();
 839
 840    Double_t phi;
 841    if(charge1>0) phi = TMath::ATan2((PMu1Dimu.Vect()).Dot(yaxisCS),((PMu1Dimu.Vect()).Dot(xaxisCS)));
 842    else phi = TMath::ATan2((PMu2Dimu.Vect()).Dot(yaxisCS),((PMu2Dimu.Vect()).Dot(xaxisCS)));
 843
 844    return phi;
 845 }
 846 //________________________________________________________________________
 847
 848 //________________________________________________________________________
 849 Double_t AliCFMuonResUpsilon::PhiHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 850 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2, Double_t Energy)
 851 {
 852    TLorentzVector PMu1CM, PMu2CM, PProjCM, PTargCM, PDimuCM; // In the CM frame
 853    TLorentzVector PMu1Dimu, PMu2Dimu, PProjDimu, PTargDimu; // In the dimuon rest frame
 854    TVector3 beta,xaxis,yaxis,zaxis;
 855    Double_t mp=0.93827231;
 856
 857    //
 858    // --- Get the muons parameters in the CM frame
 859    //
 860    PMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 861    PMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 862    //
 863    // --- Obtain the dimuon parameters in the CM frame
 864    //
 865    PDimuCM=PMu1CM+PMu2CM;
 866    //
 867    // --- Translate the muon parameters in the dimuon rest frame
 868    //
 869    zaxis=(PDimuCM.Vect()).Unit();
 870
 871    beta=(-1./PDimuCM.E())*PDimuCM.Vect();
 872
 873    PProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
 874    PTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
 875
 876    PProjDimu=PProjCM;
 877    PTargDimu=PTargCM;
 878
 879    PProjDimu.Boost(beta);
 880    PTargDimu.Boost(beta);
 881
 882    yaxis=((PProjDimu.Vect()).Cross(PTargDimu.Vect())).Unit();
 883    xaxis=(yaxis.Cross(zaxis)).Unit();
 884
 885    PMu1Dimu=PMu1CM;
 886    PMu2Dimu=PMu2CM;
 887    PMu1Dimu.Boost(beta);
 888    PMu2Dimu.Boost(beta);
 889
 890    Double_t phi;
 891    if(charge1>0) phi = TMath::ATan2((PMu1Dimu.Vect()).Dot(yaxis),(PMu1Dimu.Vect()).Dot(xaxis));
 892    else phi = TMath::ATan2((PMu2Dimu.Vect()).Dot(yaxis),(PMu2Dimu.Vect()).Dot(xaxis));
 893
 894    return phi;
 895 }
 896
 897 //________________________________________________________________________
 898 void AliCFMuonResUpsilon::Terminate(Option_t *)
 899 {
 900   // draw result of the Invariant mass MC and ESD
 901
 902 /*
 903                 TH1D *h1 = dynamic_cast<TH1D*>(GetOutputData(1));
 904     AliCFContainer *cont = dynamic_cast<AliCFContainer*> (GetOutputData(2));
 905
 906     TH1D *kmass = cont->ShowProjection(3,0);
 907     TH1D *rmass = cont->ShowProjection(3,1);
 908                 TH1D *mmass = cont->ShowProjection(3,4);
 909
 910     TCanvas *c1 = new TCanvas("AliCFMuonResUpsilon","UPSILON Container",0,0,800,800);
 911     c1->Divide(2,2);
 912     c1->cd(1);
 913     kmass->Draw("HIST");
 914     c1->cd(2);
 915     rmass->Draw("HIST");
 916                 c1->cd(3);
 917                 mmass->Draw("HIST");
 918                 c1->cd(4);
 919                 h1->Draw("HIST");
 920                 */
 921 }
 922 //________________________________________________________________________
 923
 924 #endif