PWG/muon/AliAnalysisTaskDimuonCFContainerBuilder.cxx

   1 #ifndef ALIANALYSISTASKDIMUONCFCONTAINERBUILDER_CXX
   2 #define ALIANALYSISTASKDIMUONCFCONTAINERBUILDER_CXX
   3
   4 /* $Id$ */
   5
   6 #include "AliAnalysisTaskDimuonCFContainerBuilder.h"
   7 #include "AliStack.h"
   8 #include "TParticle.h"
   9 #include "TString.h"
  10 #include "TLorentzVector.h"
  11 #include "TH2D.h"
  12 #include "AliAnalysisManager.h"
  13 #include "AliESDEvent.h"
  14 #include "AliAODEvent.h"
  15 #include "AliCFManager.h"
  16 #include "AliCFContainer.h"
  17 #include "AliESDMuonTrack.h"
  18 #include "AliESDInputHandler.h"
  19 #include "AliAODInputHandler.h"
  20 #include "AliAODMCParticle.h"
  21
  22 //      Analysis task for the building of a dimuon CF container
  23 //      Also some single-muon variables are stored
  24 //      All the variable ranges and binning are set in the AddTask macro
  25 //
  26 //      author: L. Bianchi - Universita' & INFN Torino
  27
  28
  29 ClassImp(AliAnalysisTaskDimuonCFContainerBuilder)
  30
  31 //__________________________________________________________________________
  32 AliAnalysisTaskDimuonCFContainerBuilder::AliAnalysisTaskDimuonCFContainerBuilder() :
  33   fReadAODData(0),
  34   fReadMCInfo(kTRUE),
  35   fIsAccProduction(kTRUE),
  36   fCFManager(0x0),
  37   fQAHistList(0x0),
  38   fNevt(0),
  39   fBeamEnergy(3500.),
  40   fOutput(0x0),
  41   fCutOnzVtxSPD(kFALSE),
  42   fNContributors(),
  43   fCutOnNContributors(kFALSE),
  44   fTrigClassMuon(""),
  45   fTrigClassInteraction(""),
  46   fDistinguishTrigClass(kFALSE)
  47 {
  48
  49   //Default constructor
  50
  51   Double_t chilims[2]={0.,10000.};
  52   Double_t ptlims[2]={0.,100.};
  53   Double_t thetalims[2]={0.,180.};
  54   Double_t vtxlims[2]={-1000.,1000.};
  55   SetChi2Limits(chilims);
  56   SetChi2MatchLimits(chilims);
  57   SetPtSingMuLimits(ptlims);
  58   SetThetaSingMuLimits(thetalims);
  59   SetZprimVertLimits(vtxlims);
  60   SetTrigClassMuonName();
  61   SetTrigClassInteracName();
  62   TString namemuonside[4]={"CMUS1A-","CMUS1B-","CMUS1C-","CMUS1-E-"};
  63   TString nameinteractionside[4]={"CINT1A-","CINT1B-","CINT1C-","CINT1-E"};
  64   SetTrigClassMuonSideName(namemuonside);
  65   SetTrigClassInteracSideName(nameinteractionside);
  66 }
  67 //___________________________________________________________________________
  68 AliAnalysisTaskDimuonCFContainerBuilder::AliAnalysisTaskDimuonCFContainerBuilder(const Char_t* name, Bool_t readaod, Bool_t readMC, Bool_t isaccept, Double_t beamEn) :
  69   AliAnalysisTaskSE(name),
  70   fReadAODData(0),
  71   fReadMCInfo(kTRUE),
  72   fIsAccProduction(kTRUE),
  73   fCFManager(0x0),
  74   fQAHistList(0x0),
  75   fNevt(0),
  76   fBeamEnergy(3500.),
  77   fOutput(0x0),
  78   fCutOnzVtxSPD(kFALSE),
  79   fNContributors(),
  80   fCutOnNContributors(kFALSE),
  81   fTrigClassMuon(""),
  82   fTrigClassInteraction(""),
  83   fDistinguishTrigClass(kFALSE)
  84 {
  85   //
  86   // Constructor. Initialization of Inputs and Outputs
  87   //
  88   Info("AliAnalysisTaskDimuonCFContainerBuilder","Calling Constructor");
  89
  90   SetReadAODData(readaod);
  91   SetReadMCinfo(readMC);
  92   SetIsAccProd(isaccept);
  93   SetBeamEnergy(beamEn);
  94
  95   Double_t chilims[2]={0.,10000.};
  96   Double_t ptlims[2]={0.,100.};
  97   Double_t thetalims[2]={0.,180.};
  98   Double_t vtxlims[2]={-1000.,1000.};
  99   SetChi2Limits(chilims);
 100   SetChi2MatchLimits(chilims);
 101   SetPtSingMuLimits(ptlims);
 102   SetThetaSingMuLimits(thetalims);
 103   SetZprimVertLimits(vtxlims);
 104   SetTrigClassMuonName();
 105   SetTrigClassInteracName();
 106   TString namemuonside[4]={"CMUS1A-","CMUS1B-","CMUS1C-","CMUS1-E-"};
 107   TString nameinteractionside[4]={"CINT1A-","CINT1B-","CINT1C-","CINT1-E"};
 108   SetTrigClassMuonSideName(namemuonside);
 109   SetTrigClassInteracSideName(nameinteractionside);
 110
 111   DefineOutput(1,TList::Class());
 112   DefineOutput(2,AliCFContainer::Class());
 113
 114 }
 115
 116 //___________________________________________________________________________
 117 AliAnalysisTaskDimuonCFContainerBuilder& AliAnalysisTaskDimuonCFContainerBuilder::operator=(const AliAnalysisTaskDimuonCFContainerBuilder& c)
 118 {
 119   //
 120   // Assignment operator
 121   //
 122   if (this!=&c) {
 123     AliAnalysisTaskSE::operator=(c) ;
 124     fReadAODData = c.fReadAODData ;
 125     fCFManager  = c.fCFManager;
 126     fQAHistList = c.fQAHistList ;
 127     fNevt = c.fNevt ;
 128   }
 129   return *this;
 130 }
 131
 132 //___________________________________________________________________________
 133 AliAnalysisTaskDimuonCFContainerBuilder::AliAnalysisTaskDimuonCFContainerBuilder(const AliAnalysisTaskDimuonCFContainerBuilder& c) :
 134   AliAnalysisTaskSE(c),
 135   fReadAODData(c.fReadAODData),
 136   fReadMCInfo(c.fReadMCInfo),
 137   fIsAccProduction(c.fIsAccProduction),
 138   fCFManager(c.fCFManager),
 139   fQAHistList(c.fQAHistList),
 140   fNevt(c.fNevt),
 141   fBeamEnergy(c.fBeamEnergy),
 142   fOutput(c.fOutput),
 143   fCutOnzVtxSPD(c.fCutOnzVtxSPD),
 144   fCutOnNContributors(c.fCutOnNContributors),
 145   fTrigClassMuon(c.fTrigClassMuon),
 146   fTrigClassInteraction(c.fTrigClassInteraction),
 147   fDistinguishTrigClass(c.fDistinguishTrigClass)
 148 {
 149
 150   // Copy Constructor
 151
 152 }
 153
 154 //___________________________________________________________________________
 155 AliAnalysisTaskDimuonCFContainerBuilder::~AliAnalysisTaskDimuonCFContainerBuilder() {
 156   //
 157   //destructor
 158   //
 159   Info("~AliAnalysisTaskDimuonCFContainerBuilder","Calling Destructor");
 160   if (fCFManager)           delete fCFManager ;
 161   if (fQAHistList) {fQAHistList->Clear(); delete fQAHistList;}
 162 }
 163
 164 //___________________________________________________________________________
 165 void AliAnalysisTaskDimuonCFContainerBuilder::UserCreateOutputObjects(){
 166  //     UserCreateOutputObjects
 167  fOutput = new TList();
 168  fOutput->SetOwner();
 169
 170  TH1D *hnevts   = new TH1D("hnevts","hnevts",1,0,1);                                    // Stat check histos
 171  TH1D *jpsiMult = new TH1D("jpsiMult","jpsiMult",20,0,20);
 172
 173  TH1D *ptdimuREC    = new TH1D("ptdimuREC","ptdimuREC",200,0,20);                       // Dimu check histos
 174  TH1D *ydimuREC     = new TH1D("ydimuREC","ydimuREC",200,-10.,10.);
 175  TH1D *ydimuFail     = new TH1D("ydimuFail","ydimuFail",10,660,670);
 176  TH1D *costHEdimuREC= new TH1D("costHEdimuREC","costHEdimuREC",200,-1.,1.);
 177  TH1D *costCSdimuREC= new TH1D("costCSdimuREC","costCSdimuREC",200,-1.,1.);
 178  TH1D *costHEdimuFail= new TH1D("costHEdimuFail","costHEdimuFail",10,660.,670.);
 179  TH1D *costCSdimuFail= new TH1D("costCSdimuFail","costCSdimuFail",10,660.,670.);
 180  TH1D *phiHEdimuREC = new TH1D("phiHEdimuREC","phiHEdimuREC",100,0.,TMath::Pi());
 181  TH1D *phiCSdimuREC = new TH1D("phiCSdimuREC","phiCSdimuREC",100,0.,TMath::Pi());
 182  TH1D *phiHEdimuFail = new TH1D("phiHEdimuFail","phiHEdimuFail",10,660.,670.);
 183  TH1D *phiCSdimuFail = new TH1D("phiCSdimuFail","phiCSdimuFail",10,660.,670.);
 184  TH1D *imassTot     = new TH1D("imassTot","imassTot",100,0,8);
 185  TH1D *trigCond     = new TH1D("trigCond","trigCond",40,0,4);
 186
 187  TH1D *emuonREC = new TH1D("emuonREC","emuonREC",200,0.,20.);                           // Mu check histos
 188  TH1D *ptmuonREC= new TH1D("ptmuonREC","ptmuonREC",200,0.,20.);
 189  TH1D *ymuonREC = new TH1D("ymuonREC","ymuonREC",200,-10.,10.);
 190  TH1D *hdca     = new TH1D("hdca","hdca",20,0,200);
 191  TH2D *hdcay    = new TH2D("hdcay","hdcay",200,0,200,20,-5,0);
 192
 193  TH1D *zvSPD    = new TH1D("zvSPD","zvSPD",100,-50,50.);                                // Event check histos
 194  TH1D *zvSPDcut = new TH1D("zvSPDcut","zvSPDcut",100,-50,50.);
 195  TH1D *nContrib = new TH1D("nContrib","nContrib",100,0,100);
 196
 197
 198  fOutput->Add(hnevts);
 199  fOutput->Add(jpsiMult);
 200
 201  fOutput->Add(ptdimuREC);
 202  fOutput->Add(ydimuREC);
 203  fOutput->Add(ydimuFail);
 204  fOutput->Add(costHEdimuREC);
 205  fOutput->Add(costCSdimuREC);
 206  fOutput->Add(costHEdimuFail);
 207  fOutput->Add(costCSdimuFail);
 208  fOutput->Add(phiHEdimuREC);
 209  fOutput->Add(phiCSdimuREC);
 210  fOutput->Add(phiHEdimuFail);
 211  fOutput->Add(phiCSdimuFail);
 212  fOutput->Add(imassTot);
 213  fOutput->Add(trigCond);
 214
 215  fOutput->Add(emuonREC);
 216  fOutput->Add(ptmuonREC);
 217  fOutput->Add(ymuonREC);
 218  fOutput->Add(hdca);
 219  fOutput->Add(hdcay);
 220
 221  fOutput->Add(zvSPD);
 222  fOutput->Add(zvSPDcut);
 223  fOutput->Add(nContrib);
 224
 225
 226  fOutput->ls();
 227
 228
 229 }
 230
 231
 232 //_________________________________________________
 233 void AliAnalysisTaskDimuonCFContainerBuilder::UserExec(Option_t *)
 234 {
 235
 236
 237   //Info("UserExec","");
 238   fNevt++;
 239   ((TH1D*)(fOutput->FindObject("hnevts")))->Fill(0.5);
 240
 241   Double_t containerInput[15]={666,666,666,666,666,666,666,666,666,666,666,666,666,666,666};   //y, pT, costHE, phiHE, costCS, phiCS, mass, TrigCond
 242   Int_t cutAccept =1;
 243   Int_t numbJpsis =0;
 244
 245   if (!fReadAODData){     // ESD-based ANALYSIS
 246
 247     if (fReadMCInfo){
 248
 249       if (!fMCEvent) {
 250         Error("UserExec","NO MC EVENT FOUND!");
 251         //return;
 252       }
 253
 254       // MC part  ---------------------------------------------------------------------------------------------------
 255
 256       //fCFManager->SetEventInfo(fMCEvent);                                     CHANGES IN NEW VERSIONS - MANUAL CUT ON PDG
 257       AliStack* stack = fMCEvent->Stack();
 258
 259       // loop on the MC event
 260       for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) {
 261
 262         AliMCParticle *mcPart  = (AliMCParticle*) fMCEvent->GetTrack(ipart);
 263
 264         TParticle *part = mcPart->Particle();
 265
 266         // Mother kinematics
 267         Double_t e = part->Energy();
 268         Double_t pz = part->Pz();
 269         Double_t rapmc = Rap(e,pz);
 270
 271         // Selection of the resonance
 272         //if (!fCFManager->CheckParticleCuts(0,mcPart)) continue;
 273         if (part->GetPdgCode()!=443) continue;                                  // MANUAL CUT ON PDG
 274         numbJpsis++;
 275
 276         // Decays kinematics
 277         Int_t p0 = part->GetDaughter(0);
 278         TParticle *part0 = stack->Particle(p0);
 279         Int_t pdg0 = part0->GetPdgCode();
 280
 281         Int_t p1 = part->GetDaughter(1);
 282         TParticle *part1 = stack->Particle(p1);
 283         Int_t pdg1 = part1->GetPdgCode();
 284
 285         Double_t e0 = part0->Energy();
 286         Double_t pz0 = part0->Pz();
 287         Double_t py0 = part0->Py();
 288         Double_t px0 = part0->Px();
 289         Double_t charge0 = part0->GetPDG()->Charge()/3;
 290         Double_t theta0 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(px0*px0+py0*py0),pz0);
 291         Double_t pt0 = TMath::Sqrt(px0*px0+py0*py0);
 292         Double_t mom0 = part0->P();
 293
 294         Double_t e1 = part1->Energy();
 295         Double_t pz1 = part1->Pz();
 296         Double_t py1 = part1->Py();
 297         Double_t px1 = part1->Px();
 298         //Double_t charge1 = part1->GetPDG()->Charge()/3;
 299         Double_t theta1 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(px1*px1+py1*py1),pz1);
 300         Double_t pt1 = TMath::Sqrt(px1*px1+py1*py1);
 301         Double_t mom1 = part1->P();
 302
 303
 304         if(pdg0==13 || pdg1==13) {
 305
 306           Double_t ptmc = TMath::Sqrt((px0+px1)*(px0+px1)+(py0+py1)*(py0+py1));
 307           Double_t imassmc = Imass(e0,px0,py0,pz0,e1,px1,py1,pz1);
 308
 309           Double_t costCS = CostCS(px0,py0,pz0,e0,charge0,px1,py1,pz1,e1);
 310           Double_t costHE = CostHE(px0,py0,pz0,e0,charge0,px1,py1,pz1,e1);
 311           Double_t phiCS  = PhiCS(px0,py0,pz0,e0,charge0,px1,py1,pz1,e1);
 312           Double_t phiHE  = PhiHE(px0,py0,pz0,e0,charge0,px1,py1,pz1,e1);
 313
 314           containerInput[0]  = rapmc ;
 315           containerInput[1]  = ptmc;
 316           containerInput[2]  = costHE;
 317           containerInput[3]  = TMath::Abs(phiHE);
 318           containerInput[4]  = costCS;
 319           containerInput[5]  = TMath::Abs(phiCS);
 320           containerInput[6]  = imassmc;
 321           containerInput[7]  = 1.;              //for generated no trigger condition
 322           if (pt0<pt1) {
 323             containerInput[8]=pt0;
 324             containerInput[9]=pt1;
 325           } else {
 326             containerInput[8]=pt1;
 327             containerInput[9]=pt0;
 328           }
 329           if (theta0<theta1) {
 330             containerInput[10]=theta0;
 331             containerInput[11]=theta1;
 332           } else {
 333             containerInput[10]=theta1;
 334             containerInput[11]=theta0;
 335           }
 336           if (mom0<mom1) {
 337             containerInput[12]=mom0;
 338             containerInput[13]=mom1;
 339           } else {
 340             containerInput[12]=mom1;
 341             containerInput[13]=mom0;
 342           }
 343           containerInput[14]=1.;
 344
 345           // fill the container at the first step
 346           fCFManager->GetParticleContainer()->Fill(containerInput,0);
 347
 348           if(fIsAccProduction){
 349             // acceptance cuts on single mu
 350             if(theta0<fThetaSingMuCut[0]||theta0>fThetaSingMuCut[1]||theta1<fThetaSingMuCut[0]||theta1>fThetaSingMuCut[1]) cutAccept=0;
 351             if(pt0<fPtSingMuCut[0] || pt0>fPtSingMuCut[1] || pt1<fPtSingMuCut[0] || pt1>fPtSingMuCut[1]) cutAccept=0;
 352             if (cutAccept ==1) fCFManager->GetParticleContainer()->Fill(containerInput,2);
 353           }
 354         }
 355       }
 356     } //fReadMCInfo
 357
 358
 359       // ESD part  ---------------------------------------------------------------------------------------------------
 360
 361       AliESDEvent *fESD;
 362       AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 363       if ( ! esdH ) {
 364         AliError("Cannot get input event handler");
 365         return;
 366       }
 367       fESD = esdH->GetEvent();
 368       Int_t mult1 = fESD->GetNumberOfMuonTracks() ;
 369
 370       Int_t trigfired=-1;
 371       Int_t trigside=-1;
 372       if(fDistinguishTrigClass){
 373         //aod->GetHeader()->SetFiredTriggerClasses("CMU");
 374         TString trigclass = fESD->GetFiredTriggerClasses();
 375         printf("TrigClass: %s\n",trigclass.Data());
 376         if(trigclass.Contains(fTrigClassMuon)) trigfired = 1;
 377         else if (trigclass.Contains(fTrigClassInteraction)) trigfired = 0;
 378         printf("Muon container: %d\n",trigfired);
 379         for(Int_t i=0;i<4;i++){
 380           if(trigfired==1 && trigclass.Contains(fTrigClassMuonSide[i])) {trigside=i;printf("Muon Side: %d\n",trigside);}
 381           if(trigfired==0 && trigclass.Contains(fTrigClassInteractionSide[i])) {trigside=i;printf("Interaction Side: %d\n",trigside);}
 382         }
 383       }
 384
 385       ((TH1D*)(fOutput->FindObject("zvSPD")))->Fill(fESD->GetPrimaryVertexSPD()->GetZ());
 386       ((TH1D*)(fOutput->FindObject("nContrib")))->Fill(fESD->GetPrimaryVertexSPD()->GetNContributors());
 387       if (!fCutOnzVtxSPD || (fCutOnzVtxSPD && (fESD->GetPrimaryVertexSPD()->GetZ()>fzPrimVertexSPD[0] && fESD->GetPrimaryVertexSPD()->GetZ()<fzPrimVertexSPD[1]))){
 388       ((TH1D*)(fOutput->FindObject("zvSPDcut")))->Fill(fESD->GetPrimaryVertexSPD()->GetZ());
 389       if (!fCutOnNContributors || (fCutOnNContributors && (fESD->GetPrimaryVertexSPD()->GetNContributors()>fNContributors[0] && fESD->GetPrimaryVertexSPD()->GetNContributors()<fNContributors[1]))){
 390       for (Int_t j = 0; j<mult1; j++) {
 391
 392         AliESDMuonTrack* mu1 = new AliESDMuonTrack(*(fESD->GetMuonTrack(j)));
 393         if (!mu1->ContainTrackerData()) continue;
 394         if (mu1->GetChi2()<fChi2Track[0] || mu1->GetChi2()>fChi2Track[1]) continue;
 395         if (mu1->GetChi2MatchTrigger()<fChi2MatchTrig[0] || mu1->GetChi2MatchTrigger()>fChi2MatchTrig[1]) continue;
 396         Double_t zr1 = mu1->Charge();
 397         Double_t pxr1  = mu1->Px();
 398         Double_t pyr1  = mu1->Py();
 399         Double_t pzr1  = mu1->Pz();
 400         Double_t ptmu1 = TMath::Sqrt(pxr1*pxr1+pyr1*pyr1);
 401         Double_t er1 = mu1->E();
 402         Double_t mom1 = mu1->P();
 403         Double_t theta1 = (180./TMath::Pi())*mu1->Theta();
 404         Double_t rapiditymu1 = Rap(er1,pzr1);
 405         ((TH1D*)(fOutput->FindObject("emuonREC")))->Fill(er1);
 406         ((TH1D*)(fOutput->FindObject("ptmuonREC")))->Fill(ptmu1);
 407         ((TH1D*)(fOutput->FindObject("ymuonREC")))->Fill(rapiditymu1);
 408
 409         if(zr1<0){
 410
 411             for (Int_t jj = 0; jj<mult1; jj++) {
 412
 413               AliESDMuonTrack* mu2 = new AliESDMuonTrack(*(fESD->GetMuonTrack(jj)));
 414               if (!mu2->ContainTrackerData()) continue;
 415               if (mu2->GetChi2()<fChi2Track[0] || mu2->GetChi2()>fChi2Track[1]) continue;
 416               if (mu2->GetChi2MatchTrigger()<fChi2MatchTrig[0] || mu2->GetChi2MatchTrigger()>fChi2MatchTrig[1]) continue;
 417               Double_t zr2 = mu2->Charge();
 418
 419               if(zr2>0){
 420                 Double_t trigCondition=0;
 421                 if (mu1->GetMatchTrigger()>=mu2->GetMatchTrigger()) trigCondition = mu1->GetMatchTrigger()+0.1*mu2->GetMatchTrigger();
 422                 else trigCondition = mu2->GetMatchTrigger()+0.1*mu1->GetMatchTrigger();
 423                 ((TH1D*)(fOutput->FindObject("trigCond")))->Fill(trigCondition);
 424                 Double_t pxr2 = mu2->Px();
 425                 Double_t pyr2 = mu2->Py();
 426                 Double_t pzr2 = mu2->Pz();
 427                 Double_t ptmu2 = TMath::Sqrt(pxr2*pxr2+pyr2*pyr2);
 428                 Double_t er2 = mu2->E();
 429                 Double_t mom2 = mu2->P();
 430                 Double_t theta2 = (180./TMath::Pi())*mu2->Theta();
 431
 432                 Double_t ptrec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2));
 433                 ((TH1D*)(fOutput->FindObject("ptdimuREC")))->Fill(ptrec);
 434                 Double_t raprec= Rap((er1+er2),(pzr1+pzr2));
 435                 ((TH1D*)(fOutput->FindObject("ydimuREC")))->Fill(raprec);
 436                 Double_t imassrec = Imass(er1,pxr1,pyr1,pzr1,er2,pxr2,pyr2,pzr2);
 437                 ((TH1D*)(fOutput->FindObject("imassTot")))->Fill(imassrec);
 438
 439                 if (imassrec>12.) continue;
 440
 441                 Double_t costCSrec = CostCS(pxr1,pyr1,pzr1,er1,zr1,pxr2,pyr2,pzr2,er2);
 442                 ((TH1D*)(fOutput->FindObject("costCSdimuREC")))->Fill(costCSrec);
 443                 if((Int_t)costCSrec==666) ((TH1D*)(fOutput->FindObject("costCSdimuFail")))->Fill(costCSrec);
 444                 Double_t costHErec = CostHE(pxr1,pyr1,pzr1,er1,zr1,pxr2,pyr2,pzr2,er2);
 445                 ((TH1D*)(fOutput->FindObject("costHEdimuREC")))->Fill(costHErec);
 446                 if((Int_t)costHErec==666) ((TH1D*)(fOutput->FindObject("costHEdimuFail")))->Fill(costHErec);
 447                 Double_t phiCSrec  = PhiCS(pxr1,pyr1,pzr1,er1,zr1,pxr2,pyr2,pzr2,er2);
 448                 ((TH1D*)(fOutput->FindObject("phiCSdimuREC")))->Fill(phiCSrec);
 449                 if((Int_t)phiCSrec==666) ((TH1D*)(fOutput->FindObject("phiCSdimuFail")))->Fill(phiCSrec);
 450                 Double_t phiHErec  = PhiHE(pxr1,pyr1,pzr1,er1,zr1,pxr2,pyr2,pzr2,er2);
 451                 ((TH1D*)(fOutput->FindObject("phiHEdimuREC")))->Fill(phiHErec);
 452                 if((Int_t)phiHErec==666) ((TH1D*)(fOutput->FindObject("phiHEdimuFail")))->Fill(phiHErec);
 453
 454                 containerInput[0] = raprec ;
 455                 containerInput[1] = ptrec;
 456                 containerInput[2] = costHErec;
 457                 containerInput[3] = TMath::Abs(phiHErec);
 458                 containerInput[4] = costCSrec;
 459                 containerInput[5] = TMath::Abs(phiCSrec);
 460                 containerInput[6] = imassrec;
 461                 containerInput[7] = trigCondition+0.05;
 462                 if (ptmu1<ptmu2) {
 463                   containerInput[8]=ptmu1;
 464                   containerInput[9]=ptmu2;
 465                 } else {
 466                   containerInput[8]=ptmu2;
 467                   containerInput[9]=ptmu1;
 468                 }
 469                 if (theta1<theta2) {
 470                   containerInput[10]=theta1;
 471                   containerInput[11]=theta2;
 472                 } else {
 473                   containerInput[10]=theta2;
 474                   containerInput[11]=theta1;
 475                 }
 476                 if (mom1<mom2) {
 477                   containerInput[12]=mom1;
 478                   containerInput[13]=mom2;
 479                 } else {
 480                   containerInput[12]=mom2;
 481                   containerInput[13]=mom1;
 482                 }
 483                 if (fDistinguishTrigClass && trigside==0) containerInput[14]=0.;
 484                 else if (fDistinguishTrigClass && trigside==1) containerInput[14]=1.;
 485                 else if (fDistinguishTrigClass && trigside==2) containerInput[14]=2.;
 486                 else if (fDistinguishTrigClass && trigside==3) containerInput[14]=3.;
 487                 else containerInput[14]=0.;
 488
 489                 if(fDistinguishTrigClass){
 490                   if (trigfired==1) fCFManager->GetParticleContainer()->Fill(containerInput,5);
 491                   else if (trigfired==0) fCFManager->GetParticleContainer()->Fill(containerInput,1);
 492                 } else {
 493                   fCFManager->GetParticleContainer()->Fill(containerInput,1);
 494                   if(fIsAccProduction){
 495                     if(cutAccept==1) fCFManager->GetParticleContainer()->Fill(containerInput,3);
 496                   }
 497                 }
 498               }  // mu+ Selection
 499
 500             }      // second mu Loop
 501          }          // mu- Selection
 502      }
 503         }
 504       }
 505   } else {     // AOD-based ANALYSIS
 506
 507     AliAODEvent *aod;
 508     AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 509     if ( ! aodH ) {
 510       AliError("Cannot get input event handler");
 511       return;
 512     }
 513     aod = aodH->GetEvent();
 514     Int_t ntracks=aod->GetNumberOfTracks();
 515
 516     // MC part  ---------------------------------------------------------------------------------------------------
 517
 518     if (fReadMCInfo){
 519       TClonesArray *mcarray = dynamic_cast<TClonesArray*> (aod->FindListObject(AliAODMCParticle::StdBranchName()));  //array of MC particles in this event
 520       if ( ! mcarray ) {
 521         AliError("Cannot associate MC branch");
 522         return;
 523       }
 524       for(int ii=0;ii<mcarray->GetEntries();ii++){
 525         AliAODMCParticle *mctrack = (AliAODMCParticle*) mcarray->At(ii);
 526         if(mctrack->GetPdgCode()!=13) continue;
 527         Int_t numbMother = mctrack->GetMother();
 528         if (numbMother==-1) continue;
 529         AliAODMCParticle *mother = (AliAODMCParticle*) mcarray->At(numbMother);
 530         if (mother->GetPdgCode()!=443) continue;
 531         numbJpsis++;
 532         Int_t daught0 = mother->GetDaughter(0);
 533         Int_t daught1 = mother->GetDaughter(1);
 534         AliAODMCParticle *mcDaughter0 = (AliAODMCParticle*) mcarray->At(daught0);
 535         Double_t pxmc0 = mcDaughter0->Px();
 536         Double_t pymc0 = mcDaughter0->Py();
 537         Double_t pzmc0 = mcDaughter0->Pz();
 538         Double_t ptmc0 = mcDaughter0->Pt();
 539         Double_t mommc0 = mcDaughter0->P();
 540         Double_t emc0  = mcDaughter0->E();
 541         Double_t thetamc0 = (180./TMath::Pi())*mcDaughter0->Theta();
 542         Int_t charge0 = (Int_t) mcDaughter0->Charge()/3;
 543         AliAODMCParticle *mcDaughter1 = (AliAODMCParticle*) mcarray->At(daught1);
 544         Double_t pxmc1 = mcDaughter1->Px();
 545         Double_t pymc1 = mcDaughter1->Py();
 546         Double_t pzmc1 = mcDaughter1->Pz();
 547         Double_t ptmc1 = mcDaughter1->Pt();
 548         Double_t mommc1 = mcDaughter1->P();
 549         Double_t emc1  = mcDaughter1->E();
 550         Double_t thetamc1 = (180./TMath::Pi())*mcDaughter1->Theta();
 551         Int_t charge1 = (Int_t) mcDaughter1->Charge()/3;
 552
 553         if (charge0==charge1 || TMath::Abs(mcDaughter0->GetPdgCode())!=13 || TMath::Abs(mcDaughter1->GetPdgCode())!=13) continue;
 554         Double_t rapJpsi = Rap(mother->E(),mother->Pz());
 555         Double_t ptJpsi = mother->Pt();
 556         Double_t costHE = CostHE(pxmc0,pymc0,pzmc0,emc0,(Double_t)charge0,pxmc1,pymc1,pzmc1,emc1);
 557         Double_t phiHE  = PhiHE(pxmc0,pymc0,pzmc0,emc0,(Double_t)charge0,pxmc1,pymc1,pzmc1,emc1);
 558         Double_t costCS = CostCS(pxmc0,pymc0,pzmc0,emc0,(Double_t)charge0,pxmc1,pymc1,pzmc1,emc1);
 559         Double_t phiCS  = PhiCS(pxmc0,pymc0,pzmc0,emc0,(Double_t)charge0,pxmc1,pymc1,pzmc1,emc1);
 560         Double_t massJpsi = mother->M();
 561         containerInput[0] = rapJpsi;
 562         containerInput[1] = ptJpsi;
 563         containerInput[2] = costHE;
 564         containerInput[3] = TMath::Abs(phiHE);
 565         containerInput[4] = costCS;
 566         containerInput[5] = TMath::Abs(phiCS);
 567         containerInput[6] = massJpsi;
 568         containerInput[7] = 1.;         //for generated no trigger condition
 569         if (ptmc0<ptmc1) {
 570           containerInput[8]=ptmc0;
 571           containerInput[9]=ptmc1;
 572         } else {
 573           containerInput[8]=ptmc1;
 574           containerInput[9]=ptmc0;
 575         }
 576         if (thetamc0<thetamc1) {
 577           containerInput[10]=thetamc0;
 578           containerInput[11]=thetamc1;
 579         } else {
 580           containerInput[10]=thetamc1;
 581           containerInput[11]=thetamc0;
 582         }
 583         if (mommc0<mommc1) {
 584           containerInput[12]=mommc0;
 585           containerInput[13]=mommc1;
 586         } else {
 587           containerInput[12]=mommc1;
 588           containerInput[13]=mommc0;
 589         }
 590         containerInput[14]=1.;
 591
 592         if((Int_t)rapJpsi!=666 && (Int_t)costHE!=666 && (Int_t)phiHE!=666 && (Int_t)costCS!=666 && (Int_t)phiCS!=666){
 593           fCFManager->GetParticleContainer()->Fill(containerInput,0);
 594           if(fIsAccProduction){
 595             // acceptance cuts on single mu
 596             if(thetamc0<fThetaSingMuCut[0]||thetamc0>fThetaSingMuCut[1]||thetamc1<fThetaSingMuCut[0]||thetamc1>fThetaSingMuCut[1]) cutAccept=0;
 597             if(ptmc0<fPtSingMuCut[0] || ptmc0>fPtSingMuCut[1] || ptmc1<fPtSingMuCut[0] || ptmc1>fPtSingMuCut[1]) cutAccept=0;
 598             if (cutAccept ==1) fCFManager->GetParticleContainer()->Fill(containerInput,2);
 599           }
 600
 601         }
 602
 603         }
 604       }
 605
 606
 607       // AOD part  ---------------------------------------------------------------------------------------------------
 608
 609       Int_t trigfired=-1;
 610       Int_t trigside=-1;
 611       if(fDistinguishTrigClass){
 612         TString trigclass = aod->GetFiredTriggerClasses();
 613         printf("TrigClass: %s\n",trigclass.Data());
 614         if(trigclass.Contains(fTrigClassMuon)) trigfired = 1;
 615         else if (trigclass.Contains(fTrigClassInteraction)) trigfired = 0;
 616         printf("Muon container: %d\n",trigfired);
 617         for(Int_t i=0;i<4;i++){
 618           if(trigfired==1 && trigclass.Contains(fTrigClassMuonSide[i])) {trigside=i;printf("Muon Side: %d\n",trigside);}
 619           if(trigfired==0 && trigclass.Contains(fTrigClassInteractionSide[i])) {trigside=i;printf("Interaction Side: %d\n",trigside);}
 620         }
 621       }
 622
 623       ((TH1D*)(fOutput->FindObject("zvSPD")))->Fill(aod->GetPrimaryVertex()->GetZ());
 624       Int_t ncontr = aod->GetPrimaryVertex()->GetNContributors();
 625       ((TH1D*)(fOutput->FindObject("nContrib")))->Fill(ncontr);
 626       if (!fCutOnzVtxSPD || (fCutOnzVtxSPD && (aod->GetPrimaryVertex()->GetZ()>fzPrimVertexSPD[0] && aod->GetPrimaryVertex()->GetZ()<fzPrimVertexSPD[1]))){     //NOT REALLY SPD VERTEX
 627       ((TH1D*)(fOutput->FindObject("zvSPDcut")))->Fill(aod->GetPrimaryVertex()->GetZ());
 628       if (!fCutOnNContributors || (fCutOnNContributors && (aod->GetPrimaryVertex()->GetNContributors()>fNContributors[0] && aod->GetPrimaryVertex()->GetNContributors()<fNContributors[1]))){
 629       for (Int_t j = 0; j<ntracks; j++) {
 630         AliAODTrack *mu1 = aod->GetTrack(j);
 631         if(!mu1->IsMuonTrack()) continue;
 632         if (mu1->Chi2perNDF()<fChi2Track[0] || mu1->Chi2perNDF()>fChi2Track[1]) continue;
 633         if (mu1->GetChi2MatchTrigger()<fChi2MatchTrig[0] || mu1->GetChi2MatchTrigger()>fChi2MatchTrig[1]) continue;
 634         Double_t chargemu1 = mu1->Charge();
 635         Double_t pxmu1 = mu1->Px();
 636         Double_t pymu1 = mu1->Py();
 637         Double_t pzmu1 = mu1->Pz();
 638         Double_t emu1 = mu1->E();
 639         Double_t pmu1 = mu1->P();
 640         Double_t ptmu1 = mu1->Pt();
 641         Double_t rapiditymu1 = Rap(emu1,pzmu1);
 642         Double_t thetamu1 = (180./TMath::Pi())*mu1->Theta();
 643         Double_t rdcamu1 = mu1->DCA();
 644         ((TH1D*)(fOutput->FindObject("emuonREC")))->Fill(emu1);
 645         ((TH1D*)(fOutput->FindObject("ptmuonREC")))->Fill(ptmu1);
 646         ((TH1D*)(fOutput->FindObject("ymuonREC")))->Fill(rapiditymu1);
 647         if(chargemu1<0){
 648           for (Int_t jj = 0; jj<ntracks; jj++) {
 649             AliAODTrack *mu2 = aod->GetTrack(jj);
 650             if(!mu2->IsMuonTrack()) continue;
 651             if (mu2->Chi2perNDF()<fChi2Track[0] || mu2->Chi2perNDF()>fChi2Track[1]) continue;
 652             if (mu2->GetChi2MatchTrigger()<fChi2MatchTrig[0] || mu2->GetChi2MatchTrigger()>fChi2MatchTrig[1]) continue;
 653             Double_t chargemu2 = mu2->Charge();
 654             Double_t pxmu2 = mu2->Px();
 655             Double_t pymu2 = mu2->Py();
 656             Double_t pzmu2 = mu2->Pz();
 657             Double_t emu2 = mu2->E();
 658             Double_t pmu2 = mu2->P();
 659             Double_t ptmu2 = mu2->Pt();
 660             //Double_t rapiditymu2 = Rap(emu2,pzmu2);
 661             Double_t thetamu2 = (180./TMath::Pi())*mu2->Theta();
 662             Double_t rdcamu2 = mu2->DCA();
 663             if(chargemu2>0){
 664             if (mu1->GetMatchTrigger()>0) printf("Mu1: charge: %f, match: %d\n",chargemu1,1);
 665             else  printf("Mu1: charge: %f, match: %d\n",chargemu1,0);
 666             if (mu2->GetMatchTrigger()>0) printf("Mu2: charge: %f, match: %d\n",chargemu2,1);
 667             else printf("Mu2: charge: %f, match: %d\n",chargemu2,0);
 668               ((TH1D*)(fOutput->FindObject("hdca")))->Fill(rdcamu2);
 669               ((TH1D*)(fOutput->FindObject("hdca")))->Fill(rdcamu1);
 670               Double_t trigCondition=0;
 671               if (mu1->GetMatchTrigger()>=mu2->GetMatchTrigger()) trigCondition = mu1->GetMatchTrigger()+0.1*mu2->GetMatchTrigger();
 672                 else trigCondition = mu2->GetMatchTrigger()+0.1*mu1->GetMatchTrigger();
 673               containerInput[0] = (Double_t) Rap((emu1+emu2),(pzmu1+pzmu2));
 674               ((TH1D*)(fOutput->FindObject("ydimuREC")))->Fill(containerInput[0]);
 675               ((TH2D*)(fOutput->FindObject("hdcay")))->Fill(TMath::Max(rdcamu1,rdcamu2),containerInput[0]);
 676               //((TH2D*)(fOutput->FindObject("hdcay")))->Fill(Rdcamu2,containerInput[0]);
 677               containerInput[1] = (Double_t) TMath::Sqrt((pxmu1+pxmu2)*(pxmu1+pxmu2)+(pymu1+pymu2)*(pymu1+pymu2));
 678               ((TH1D*)(fOutput->FindObject("ptdimuREC")))->Fill(containerInput[1]);
 679               containerInput[2] = CostHE(pxmu1,pymu1,pzmu1,emu1,chargemu1,pxmu2,pymu2,pzmu2,emu2);
 680               ((TH1D*)(fOutput->FindObject("costHEdimuREC")))->Fill(containerInput[2]);
 681               if(containerInput[2]==666) ((TH1D*)(fOutput->FindObject("costHEdimuFail")))->Fill(containerInput[2]);
 682               containerInput[3] = TMath::Abs(PhiHE(pxmu1,pymu1,pzmu1,emu1,chargemu1,pxmu2,pymu2,pzmu2,emu2));
 683               ((TH1D*)(fOutput->FindObject("phiHEdimuREC")))->Fill(containerInput[3]);
 684               if(containerInput[3]==666) ((TH1D*)(fOutput->FindObject("phiHEdimuFail")))->Fill(containerInput[3]);
 685               containerInput[4] = CostCS(pxmu1,pymu1,pzmu1,emu1,chargemu1,pxmu2,pymu2,pzmu2,emu2);
 686               ((TH1D*)(fOutput->FindObject("costCSdimuREC")))->Fill(containerInput[4]);
 687               if(containerInput[4]==666) ((TH1D*)(fOutput->FindObject("costCSdimuFail")))->Fill(containerInput[4]);
 688               containerInput[5] = TMath::Abs(PhiCS(pxmu1,pymu1,pzmu1,emu1,chargemu1,pxmu2,pymu2,pzmu2,emu2));
 689               ((TH1D*)(fOutput->FindObject("phiCSdimuREC")))->Fill(containerInput[5]);
 690               if(containerInput[5]==666) ((TH1D*)(fOutput->FindObject("phiCSdimuFail")))->Fill(containerInput[5]);
 691               containerInput[6] = Imass(emu1,pxmu1,pymu1,pzmu1,emu2,pxmu2,pymu2,pzmu2);
 692               ((TH1D*)(fOutput->FindObject("imassTot")))->Fill(containerInput[6]);
 693               containerInput[7] = trigCondition+0.05;
 694               ((TH1D*)(fOutput->FindObject("trigCond")))->Fill(containerInput[7]);
 695               if (ptmu1<ptmu2) {
 696                 containerInput[8]=ptmu1;
 697                 containerInput[9]=ptmu2;
 698               } else {
 699                 containerInput[8]=ptmu2;
 700                 containerInput[9]=ptmu1;
 701               }
 702               if (thetamu1<thetamu2) {
 703                 containerInput[10]=thetamu1;
 704                 containerInput[11]=thetamu2;
 705               } else {
 706                 containerInput[10]=thetamu2;
 707                 containerInput[11]=thetamu1;
 708               }
 709               if (pmu1<pmu2) {
 710                 containerInput[12]=pmu1;
 711                 containerInput[13]=pmu2;
 712               } else {
 713                 containerInput[12]=pmu2;
 714                 containerInput[13]=pmu1;
 715               }
 716               if (fDistinguishTrigClass && trigside==0) containerInput[14]=0.;
 717               else if (fDistinguishTrigClass && trigside==1) containerInput[14]=1.;
 718               else if (fDistinguishTrigClass && trigside==2) containerInput[14]=2.;
 719               else if (fDistinguishTrigClass && trigside==3) containerInput[14]=3.;
 720               else containerInput[14]=0.;
 721
 722               if(containerInput[2]!=666 && containerInput[3]!=666 && containerInput[4]!=666 && containerInput[5]!=666){
 723                 if(fDistinguishTrigClass){
 724                   if (trigfired==1) fCFManager->GetParticleContainer()->Fill(containerInput,5);
 725                   else if (trigfired==0) fCFManager->GetParticleContainer()->Fill(containerInput,1);
 726                 } else {
 727                   fCFManager->GetParticleContainer()->Fill(containerInput,1);
 728                   if(fIsAccProduction){
 729                     if (cutAccept ==1) fCFManager->GetParticleContainer()->Fill(containerInput,3);
 730                   }
 731                 }
 732               } else if (containerInput[0]==666) ((TH1D*)(fOutput->FindObject("ydimuFail")))->Fill(containerInput[0]);
 733
 734           }
 735
 736         }
 737
 738       }
 739
 740     }
 741     }
 742     }
 743
 744
 745   }
 746
 747
 748 //  ----------
 749   if (fReadMCInfo) ((TH1D*)(fOutput->FindObject("jpsiMult")))->Fill(numbJpsis);
 750
 751   PostData(1,fOutput) ;
 752   PostData(2,fCFManager->GetParticleContainer()) ;
 753
 754
 755 }
 756
 757 //________________________________________________________________________
 758 Double_t AliAnalysisTaskDimuonCFContainerBuilder::Imass(Double_t e1, Double_t px1, Double_t py1, Double_t pz1,
 759                                    Double_t e2, Double_t px2, Double_t py2, Double_t pz2) const
 760 {
 761 // invariant mass calculation
 762     Double_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+
 763                                     (py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2)));
 764     return imassrec;
 765 }
 766
 767 //________________________________________________________________________
 768 Double_t AliAnalysisTaskDimuonCFContainerBuilder::Rap(Double_t e, Double_t pz) const
 769 {
 770 // calculate rapidity
 771     Double_t rap;
 772     if(e>TMath::Abs(pz)){                         // in order to avoid problems with AODs
 773         rap = 0.5*TMath::Log((e+pz)/(e-pz));
 774         return rap;
 775     }
 776     else{
 777         rap = 666.;
 778         return rap;
 779     }
 780 }
 781
 782 //________________________________________________________________________
 783 Double_t AliAnalysisTaskDimuonCFContainerBuilder::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 784 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
 785 {
 786 // Cosine of the theta decay angle (mu+) in the Collins-Soper frame
 787
 788   TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM. frame
 789   TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
 790   TVector3 beta,zaxisCS;
 791   Double_t mp=0.93827231;
 792   //
 793   // --- Fill the Lorentz vector for projectile and target in the CM frame
 794   //
 795   pProjCM.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
 796   pTargCM.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
 797   //
 798   // --- Get the muons parameters in the CM frame
 799   //
 800   pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 801   pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 802   //
 803   // --- Obtain the dimuon parameters in the CM frame
 804   //
 805   pDimuCM=pMu1CM+pMu2CM;
 806   //
 807   // --- Translate the dimuon parameters in the dimuon rest frame
 808   //
 809   beta=(-1./pDimuCM.E())*pDimuCM.Vect();
 810   if(beta.Mag()>=1) return 666.;
 811   pMu1Dimu=pMu1CM;
 812   pMu2Dimu=pMu2CM;
 813   pProjDimu=pProjCM;
 814   pTargDimu=pTargCM;
 815   pMu1Dimu.Boost(beta);
 816   pMu2Dimu.Boost(beta);
 817   pProjDimu.Boost(beta);
 818   pTargDimu.Boost(beta);
 819
 820   //Debugging part -------------------------------------
 821   Double_t debugProj[4]={0.,0.,0.,0.};
 822   Double_t debugTarg[4]={0.,0.,0.,0.};
 823   Double_t debugMu1[4]={0.,0.,0.,0.};
 824   Double_t debugMu2[4]={0.,0.,0.,0.};
 825   pMu1Dimu.GetXYZT(debugMu1);
 826   pMu2Dimu.GetXYZT(debugMu2);
 827   pProjDimu.GetXYZT(debugProj);
 828   pTargDimu.GetXYZT(debugTarg);
 829   if (debugProj[0]!=debugProj[0] ||debugProj[1]!=debugProj[1] || debugProj[2]!=debugProj[2] ||debugProj[3]!=debugProj[3]) return 666;
 830   if (debugTarg[0]!=debugTarg[0] ||debugTarg[1]!=debugTarg[1] || debugTarg[2]!=debugTarg[2] ||debugTarg[3]!=debugTarg[3]) return 666;
 831   if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
 832   if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
 833   //----------------------------------------------------
 834
 835   // --- Determine the z axis for the CS angle
 836   zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
 837
 838   // --- Determine the CS angle (angle between mu+ and the z axis defined above)
 839   Double_t cost;
 840
 841   if(charge1>0) {cost = zaxisCS.Dot((pMu1Dimu.Vect()).Unit());}
 842   else {cost = zaxisCS.Dot((pMu2Dimu.Vect()).Unit());}
 843
 844   return cost;
 845 }
 846
 847 //________________________________________________________________________
 848 Double_t AliAnalysisTaskDimuonCFContainerBuilder::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 849 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
 850 {
 851 // Cosine of the theta decay angle (mu+) in the Helicity frame
 852
 853   TLorentzVector pMu1CM, pMu2CM, pDimuCM; // In the CM frame
 854   TLorentzVector pMu1Dimu, pMu2Dimu; // In the dimuon rest frame
 855   TVector3 beta,zaxisCS;
 856   //
 857   // --- Get the muons parameters in the CM frame
 858   //
 859   pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 860   pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 861   //
 862   // --- Obtain the dimuon parameters in the CM frame
 863   //
 864   pDimuCM=pMu1CM+pMu2CM;
 865   //
 866   // --- Translate the muon parameters in the dimuon rest frame
 867   //
 868   beta=(-1./pDimuCM.E())*pDimuCM.Vect();
 869   if(beta.Mag()>=1) return 666.;
 870   pMu1Dimu=pMu1CM;
 871   pMu2Dimu=pMu2CM;
 872   pMu1Dimu.Boost(beta);
 873   pMu2Dimu.Boost(beta);
 874
 875   //Debugging part -------------------------------------
 876   Double_t debugMu1[4]={0.,0.,0.,0.};
 877   Double_t debugMu2[4]={0.,0.,0.,0.};
 878   pMu1Dimu.GetXYZT(debugMu1);
 879   pMu2Dimu.GetXYZT(debugMu2);
 880   if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
 881   if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
 882   //----------------------------------------------------
 883
 884   // --- Determine the z axis for the calculation of the polarization angle (i.e. the direction of the dimuon in the CM system)
 885   TVector3 zaxis;
 886   zaxis=(pDimuCM.Vect()).Unit();
 887
 888   // --- Calculation of the polarization angle (angle between mu+ and the z axis defined above)
 889   Double_t cost;
 890   if(charge1>0) {cost = zaxis.Dot((pMu1Dimu.Vect()).Unit());}
 891   else {cost = zaxis.Dot((pMu2Dimu.Vect()).Unit());}
 892   return cost;
 893 }
 894
 895 //________________________________________________________________________
 896 Double_t AliAnalysisTaskDimuonCFContainerBuilder::PhiCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 897 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
 898 {
 899 // Phi decay angle (mu+) in the Collins-Soper frame
 900
 901    TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM frame
 902    TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
 903    TVector3 beta,yaxisCS, xaxisCS, zaxisCS;
 904    Double_t mp=0.93827231;
 905
 906    // --- Fill the Lorentz vector for projectile and target in the CM frame
 907    pProjCM.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
 908    pTargCM.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp));
 909
 910    // --- Get the muons parameters in the CM frame
 911    pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 912    pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 913
 914    // --- Obtain the dimuon parameters in the CM frame
 915    pDimuCM=pMu1CM+pMu2CM;
 916
 917    // --- Translate the dimuon parameters in the dimuon rest frame
 918    beta=(-1./pDimuCM.E())*pDimuCM.Vect();
 919    if(beta.Mag()>=1) return 666.;
 920    pMu1Dimu=pMu1CM;
 921    pMu2Dimu=pMu2CM;
 922    pProjDimu=pProjCM;
 923    pTargDimu=pTargCM;
 924    pMu1Dimu.Boost(beta);
 925    pMu2Dimu.Boost(beta);
 926    pProjDimu.Boost(beta);
 927    pTargDimu.Boost(beta);
 928
 929    //Debugging part -------------------------------------
 930    Double_t debugProj[4]={0.,0.,0.,0.};
 931    Double_t debugTarg[4]={0.,0.,0.,0.};
 932    Double_t debugMu1[4]={0.,0.,0.,0.};
 933    Double_t debugMu2[4]={0.,0.,0.,0.};
 934    pMu1Dimu.GetXYZT(debugMu1);
 935    pMu2Dimu.GetXYZT(debugMu2);
 936    pProjDimu.GetXYZT(debugProj);
 937    pTargDimu.GetXYZT(debugTarg);
 938    if (debugProj[0]!=debugProj[0] ||debugProj[1]!=debugProj[1] || debugProj[2]!=debugProj[2] ||debugProj[3]!=debugProj[3]) return 666;
 939    if (debugTarg[0]!=debugTarg[0] ||debugTarg[1]!=debugTarg[1] || debugTarg[2]!=debugTarg[2] ||debugTarg[3]!=debugTarg[3]) return 666;
 940    if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
 941    if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
 942    //----------------------------------------------------
 943
 944    // --- Determine the z axis for the CS angle
 945    zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
 946    yaxisCS=(((pProjDimu.Vect()).Unit()).Cross((pTargDimu.Vect()).Unit())).Unit();
 947    xaxisCS=(yaxisCS.Cross(zaxisCS)).Unit();
 948
 949    Double_t phi=0.;
 950    if(charge1>0) {
 951        phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxisCS),((pMu1Dimu.Vect()).Dot(xaxisCS)));
 952    } else {
 953        phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxisCS),((pMu2Dimu.Vect()).Dot(xaxisCS)));
 954    }
 955    if (phi>TMath::Pi()) phi=phi-TMath::Pi();
 956
 957    return phi;
 958 }
 959
 960 //________________________________________________________________________
 961 Double_t AliAnalysisTaskDimuonCFContainerBuilder::PhiHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 962 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
 963 {
 964 // Phi decay angle (mu+) in the Helicity frame
 965   TLorentzVector pMu1Lab, pMu2Lab, pProjLab, pTargLab, pDimuLab; // In the lab. frame
 966   TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame
 967   TVector3 beta,xaxis, yaxis,zaxis;
 968   Double_t mp=0.93827231;
 969
 970   // --- Get the muons parameters in the LAB frame
 971   pMu1Lab.SetPxPyPzE(px1,py1,pz1,e1);
 972   pMu2Lab.SetPxPyPzE(px2,py2,pz2,e2);
 973
 974   // --- Obtain the dimuon parameters in the LAB frame
 975   pDimuLab=pMu1Lab+pMu2Lab;
 976   zaxis=(pDimuLab.Vect()).Unit();
 977
 978   // --- Translate the muon parameters in the dimuon rest frame
 979   beta=(-1./pDimuLab.E())*pDimuLab.Vect();
 980   if(beta.Mag()>=1.) return 666.;
 981
 982   pProjLab.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp)); // proiettile
 983   pTargLab.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp)); // bersaglio
 984
 985   pProjDimu=pProjLab;
 986   pTargDimu=pTargLab;
 987
 988   pProjDimu.Boost(beta);
 989   pTargDimu.Boost(beta);
 990
 991   yaxis=((pProjDimu.Vect()).Cross(pTargDimu.Vect())).Unit();
 992   xaxis=(yaxis.Cross(zaxis)).Unit();
 993
 994   pMu1Dimu=pMu1Lab;
 995   pMu2Dimu=pMu2Lab;
 996   pMu1Dimu.Boost(beta);
 997   pMu2Dimu.Boost(beta);
 998
 999   //Debugging part -------------------------------------
1000   Double_t debugProj[4]={0.,0.,0.,0.};
1001   Double_t debugTarg[4]={0.,0.,0.,0.};
1002   Double_t debugMu1[4]={0.,0.,0.,0.};
1003   Double_t debugMu2[4]={0.,0.,0.,0.};
1004   pMu1Dimu.GetXYZT(debugMu1);
1005   pMu2Dimu.GetXYZT(debugMu2);
1006   pProjDimu.GetXYZT(debugProj);
1007   pTargDimu.GetXYZT(debugTarg);
1008   if (debugProj[0]!=debugProj[0] ||debugProj[1]!=debugProj[1] || debugProj[2]!=debugProj[2] ||debugProj[3]!=debugProj[3]) return 666;
1009   if (debugTarg[0]!=debugTarg[0] ||debugTarg[1]!=debugTarg[1] || debugTarg[2]!=debugTarg[2] ||debugTarg[3]!=debugTarg[3]) return 666;
1010   if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666;
1011   if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666;
1012   //----------------------------------------------------
1013
1014   Double_t phi=0.;
1015    if(charge1 > 0) {
1016       phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxis),(pMu1Dimu.Vect()).Dot(xaxis));
1017      } else {
1018       phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxis),(pMu2Dimu.Vect()).Dot(xaxis));
1019    }
1020    return phi;
1021 }
1022
1023 //________________________________________________________________________
1024 void AliAnalysisTaskDimuonCFContainerBuilder::Terminate(Option_t *)
1025 {
1026
1027 }
1028
1029 #endif