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