PWG/muon/AliAnalysisMuMuMinv.cxx

   1 #include "AliAnalysisMuMuMinv.h"
   2
   3 /**
   4  *
   5  * \ingroup pwg-muon-mumu
   6  *
   7  * \class AliAnalysisMuMuMinv
   8  *
   9  * Analysis which fills a bunch of histograms for invariant mass analysis of J/psi
  10  *
  11  * Can be used on real data and/or simulated (for instance to get Acc x Eff)
  12  *
  13  * Can optionally use as input an already computed Acc x Eff matrix that will be applied
  14  * when filling the invariant mass histograms.
  15  *
  16  */
  17
  18 #include "TH2F.h"
  19 #include "AliLog.h"
  20 #include "TObjArray.h"
  21 #include "AliAnalysisMuMuBinning.h"
  22 #include "TString.h"
  23 #include "TLorentzVector.h"
  24 #include "TString.h"
  25 #include "AliMCEvent.h"
  26 #include "AliMergeableCollection.h"
  27 #include "AliAnalysisMuonUtility.h"
  28 #include "TParameter.h"
  29 #include <cassert>
  30
  31 ClassImp(AliAnalysisMuMuMinv)
  32
  33 //_____________________________________________________________________________
  34 AliAnalysisMuMuMinv::AliAnalysisMuMuMinv(TH2* accEffHisto, Bool_t computeMeanPt, Int_t systLevel)
  35 : AliAnalysisMuMuBase(),
  36 fcomputeMeanPt(computeMeanPt),
  37 fAccEffHisto(0x0),
  38 fsystLevel(systLevel)
  39 {
  40   // FIXME : find the AccxEff histogram from HistogramCollection()->Histo("/EXCHANGE/JpsiAccEff")
  41
  42   if ( accEffHisto )
  43   {
  44     fAccEffHisto = static_cast<TH2F*>(accEffHisto->Clone());
  45     fAccEffHisto->SetDirectory(0);
  46   }
  47 }
  48
  49 //_____________________________________________________________________________
  50 AliAnalysisMuMuMinv::~AliAnalysisMuMuMinv()
  51 {
  52   /// dtor
  53   delete fAccEffHisto;
  54 }
  55
  56 //_____________________________________________________________________________
  57 void
  58 AliAnalysisMuMuMinv::DefineHistogramCollection(const char* eventSelection,
  59                                                const char* triggerClassName,
  60                                                const char* centrality)
  61 {
  62   /// Define the histograms this analysis will use
  63
  64   if ( Histo(eventSelection,triggerClassName,centrality,"AliAnalysisMuMuMinv") )
  65   {
  66     return;
  67   }
  68
  69   // dummy histogram to signal that we already defined all our histograms (see above)
  70   CreateEventHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"AliAnalysisMuMuMinv","Dummy semaphore",1,0,1);
  71
  72   /// Create invariant mass histograms
  73
  74   Double_t minvMin = 0;
  75   Double_t minvMax = 16;
  76   Int_t nMinvBins = GetNbins(minvMin,minvMax,0.025);
  77
  78   Int_t nMCMinvBins = GetNbins(minvMin,minvMax,0.1);
  79
  80   Double_t rapidityMin = -5;
  81   Double_t rapidityMax = -2;
  82   Int_t nbinsRapidity = GetNbins(rapidityMin,rapidityMax,0.05);
  83
  84   Double_t etaMin = -5;
  85   Double_t etaMax = -2;
  86   Int_t nbinsEta = GetNbins(etaMin,etaMax,0.05);
  87
  88   TObjArray* bins = Binning()->CreateBinObjArray("psi","integrated,ptvsy,yvspt,pt,y,phi,nch,dnchdeta,ntrcorr,ntrcorrpt,ntrcorry,relntrcorr","");//We may include ,v0a,v0acent
  89
  90   CreatePairHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"Pt","#mu+#mu- Pt distribution",
  91                   200,0,20,-2);
  92
  93   CreatePairHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"Y","#mu+#mu- Y distribution",
  94                    nbinsRapidity,rapidityMin,rapidityMax,-2);
  95
  96   CreatePairHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,"Eta","#mu+#mu- Eta distribution",
  97                    nbinsEta,etaMin,etaMax);
  98
  99
 100   //___Histos for pure MC
 101   CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"Pt","MCINPUT #mu+#mu- Pt distribution",
 102                    200,0,20,-2);
 103
 104   CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"Y","MCINPUT #mu+#mu- Y distribution",
 105                    nbinsRapidity,rapidityMin,rapidityMax,-2);
 106
 107   CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,"Eta","MCINPUT #mu+#mu- Eta distribution",
 108                    nbinsEta,etaMin,etaMax);
 109
 110   CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,Form("%s/INYRANGE",centrality),"Pt","MCINPUT #mu+#mu- Pt distribution",
 111                     200,0,20,-2);
 112
 113   CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,Form("%s/INYRANGE",centrality),"Y","MCINPUT #mu+#mu- Y distribution",
 114                     nbinsRapidity,rapidityMin,rapidityMax,-2);
 115
 116   CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,Form("%s/INYRANGE",centrality),"Eta","MCINPUT #mu+#mu- Eta distribution",
 117                     nbinsEta,etaMin,etaMax);
 118   //____
 119
 120   //  CreatePairHistos(eventSelection,triggerClassName,centrality,"BinFlowPt","#mu+#mu- BinFlowPt distribution",
 121   //                  200,0,20);
 122
 123   CreatePairHistos(kHistoForData,eventSelection,triggerClassName,centrality,"PtRecVsSim","#mu+#mu- Pt distribution rec vs sim",
 124                   200,0,20,200,0,20);
 125
 126   TIter next(bins);
 127   AliAnalysisMuMuBinning::Range* r;
 128   Int_t nb(0);
 129
 130   while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(next()) ) )
 131   {
 132     TString minvName(GetMinvHistoName(*r,kFALSE));
 133
 134     ++nb;
 135
 136     if ( !IsHistogramDisabled(minvName.Data()) )
 137     {
 138
 139       AliDebug(1,Form("bin %d %s histoname = %s",nb,r->AsString().Data(),minvName.Data()));
 140
 141       CreatePairHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,minvName.Data(),
 142                        Form("#mu+#mu- inv. mass %s;M_{#mu^{+}#mu^{-}} (GeV/c^2);Counts",
 143                             r->AsString().Data()),
 144                        nMinvBins,minvMin,minvMax,-2);
 145
 146       CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,minvName.Data(),
 147                        Form("MCINPUT #mu+#mu- inv. mass %s;M_{#mu^{+}#mu^{-}} (GeV/c^2);Counts",
 148                             r->AsString().Data()),
 149                        nMCMinvBins,minvMin,minvMax,-2); // Pure MC histo
 150
 151       CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,Form("%s/INYRANGE",centrality),minvName.Data(),
 152                         Form("MCINPUT #mu+#mu- inv. mass %s;M_{#mu^{+}#mu^{-}} (GeV/c^2);Counts",
 153                              r->AsString().Data()),
 154                         nMCMinvBins,minvMin,minvMax,-2); // Pure MC histo
 155
 156
 157       if ( fcomputeMeanPt )
 158       {
 159         TString mPtName(Form("MeanPtVs%s",minvName.Data()));
 160         CreatePairHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,mPtName.Data(),
 161                          Form("#mu+#mu- mean p_{T} %s;M_{#mu^{+}#mu^{-}} (GeV/c^2);<p_{T}^{#mu^{+}#mu^{-} (GeV/c^2)}>",
 162                               r->AsString().Data()),
 163                          nMinvBins,minvMin,minvMax,0);
 164
 165         CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,mPtName.Data(),
 166                          Form("#mu+#mu- mean p_{T} %s;M_{#mu^{+}#mu^{-}} (GeV/c^2);<p_{T}^{#mu^{+}#mu^{-} (GeV/c^2)}>",
 167                               r->AsString().Data()),
 168                          nMinvBins,minvMin,minvMax,0); //Pure MC Histo
 169
 170         CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,Form("%s/INYRANGE",centrality),mPtName.Data(),
 171                           Form("#mu+#mu- mean p_{T} %s;M_{#mu^{+}#mu^{-}} (GeV/c^2);<p_{T}^{#mu^{+}#mu^{-} (GeV/c^2)}>",
 172                                r->AsString().Data()),
 173                           nMinvBins,minvMin,minvMax,0); //Pure MC Histo
 174       }
 175
 176 //      if ( HasMC() )
 177 //      {
 178 //        TH1* h = new TH1F(minvName.Data(),Form("MC #mu+#mu- inv. mass %s",r->AsString().Data()),
 179 //                          nMCMinvBins,minvMin,minvMax);
 180 //
 181 //        HistogramCollection()->Adopt(Form("/%s/ALL",MCInputPrefix()),h);
 182 //
 183 //        HistogramCollection()->Adopt(Form("/%s/INYRANGE",MCInputPrefix()),static_cast<TH1*>(h->Clone()));
 184 //      }
 185     }
 186
 187     if ( ShouldCorrectDimuonForAccEff() )
 188     {
 189       minvName = GetMinvHistoName(*r,kTRUE);
 190
 191       if ( !IsHistogramDisabled(minvName.Data()) )
 192       {
 193
 194         AliDebug(1,Form("bin %d %s histoname = %s",nb,r->AsString().Data(),minvName.Data()));
 195
 196         CreatePairHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,minvName.Data(),
 197                          Form("#mu+#mu- inv. mass %s (Acc #times Eff Corrected);M_{#mu^{+}#mu^{-}}(GeV/c^2);Counts",
 198                               r->AsString().Data()),
 199                          nMinvBins,minvMin,minvMax,-2);
 200
 201         CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,centrality,minvName.Data(),
 202                          Form("#mu+#mu- inv. mass %s (Acc #times Eff Corrected);M_{#mu^{+}#mu^{-}} (GeV/c^2);Counts",
 203                               r->AsString().Data()),
 204                          nMCMinvBins,minvMin,minvMax,-2); // Pure MC histo
 205
 206         CreateEventHistos(kHistoForMCInput,eventSelection,triggerClassName,Form("%s/INYRANGE",centrality),minvName.Data(),
 207                           Form("#mu+#mu- inv. mass %s (Acc #times Eff Corrected);M_{#mu^{+}#mu^{-}} (GeV/c^2);Counts",
 208                                r->AsString().Data()),
 209                           nMCMinvBins,minvMin,minvMax,-2); // Pure MC histo
 210
 211         if ( fcomputeMeanPt )
 212         {
 213           TString mPtName(Form("MeanPtVs%s",minvName.Data()));
 214           CreatePairHistos(kHistoForData | kHistoForMCInput,eventSelection,triggerClassName,centrality,mPtName.Data(),
 215                            Form("#mu+#mu- mean p_{T} %s (Acc #times Eff Corrected);M_{#mu^{+}#mu^{-}} (GeV/c^2);<p_{T}^{#mu^{+}#mu^{-}}>",
 216                                 r->AsString().Data()),
 217                            nMinvBins,minvMin,minvMax,0);
 218         }
 219
 220 //        if ( HasMC() )
 221 //        {
 222 //          TH1*  h = new TH1F(minvName.Data(),Form("MC #mu+#mu- inv. mass %s",r->AsString().Data()),
 223 //                             nMCMinvBins,minvMin,minvMax);
 224 //
 225 //          h->Sumw2();
 226 //
 227 //          HistogramCollection()->Adopt(Form("/%s/ALL",MCInputPrefix()),h);
 228 //
 229 //          HistogramCollection()->Adopt(Form("/%s/INYRANGE",MCInputPrefix()),static_cast<TH1*>(h->Clone()));
 230 //        }
 231       }
 232     }
 233   }
 234
 235   delete bins;
 236 }
 237
 238 //_____________________________________________________________________________
 239 void AliAnalysisMuMuMinv::FillHistosForPair(const char* eventSelection, const char* triggerClassName,
 240                                             const char* centrality, const char* pairCutName,
 241                                             const AliVParticle& tracki,
 242                                             const AliVParticle& trackj)
 243 {
 244
 245   TLorentzVector pi(tracki.Px(),tracki.Py(),tracki.Pz(),
 246                     TMath::Sqrt(AliAnalysisMuonUtility::MuonMass2()+tracki.P()*tracki.P()));
 247
 248   if (!AliAnalysisMuonUtility::IsMuonTrack(&tracki) ) return;
 249   if (!AliAnalysisMuonUtility::IsMuonTrack(&trackj) ) return;
 250
 251   TLorentzVector pair4Momentum(trackj.Px(),trackj.Py(),trackj.Pz(),
 252                                TMath::Sqrt(AliAnalysisMuonUtility::MuonMass2()+trackj.P()*trackj.P()));
 253
 254   pair4Momentum += pi;
 255
 256
 257 ////  if (!IsHistogramDisabled("Chi12"))
 258 ////  {
 259 ////    Histo(eventSelection,triggerClassName,centrality,pairCutName,"Chi12")
 260 ////    ->Fill(
 261 ////           AliAnalysisMuonUtility::GetChi2perNDFtracker(&tracki),
 262 ////           AliAnalysisMuonUtility::GetChi2perNDFtracker(&trackj));
 263 ////  }
 264 ////
 265 ////  if (!IsHistogramDisabled("Rabs12"))
 266 ////  {
 267 ////    Histo(eventSelection,triggerClassName,centrality,pairCutName,"Rabs12")
 268 ////    ->Fill(AliAnalysisMuonUtility::GetRabs(&tracki),
 269 ////           AliAnalysisMuonUtility::GetRabs(&trackj));
 270 ////  }
 271
 272   if ( ( tracki.Charge() != trackj.Charge() ) )
 273   {
 274     Double_t inputWeight = WeightDistribution(pair4Momentum.Pt(),pair4Momentum.Rapidity());
 275
 276     if ( !IsHistogramDisabled("Pt") )
 277     {
 278       Histo(eventSelection,triggerClassName,centrality,pairCutName,"Pt")->Fill(pair4Momentum.Pt(),inputWeight);
 279     }
 280     if ( !IsHistogramDisabled("Y") )
 281     {
 282       Histo(eventSelection,triggerClassName,centrality,pairCutName,"Y")->Fill(pair4Momentum.Rapidity(),inputWeight);
 283     }
 284     if ( !IsHistogramDisabled("Eta") )
 285     {
 286       Histo(eventSelection,triggerClassName,centrality,pairCutName,"Eta")->Fill(pair4Momentum.Eta());
 287     }
 288
 289     TLorentzVector* pair4MomentumMC(0x0);
 290     Double_t inputWeightMC(1.);
 291     if ( HasMC() )
 292     {
 293       Int_t labeli = tracki.GetLabel();
 294       Int_t labelj = trackj.GetLabel();
 295
 296       if ( labeli < 0 || labelj < 0 )
 297       {
 298         AliError("Got negative labels!");
 299       }
 300       else
 301       {
 302         AliVParticle* mcTracki = MCEvent()->GetTrack(labeli);
 303         AliVParticle* mcTrackj = MCEvent()->GetTrack(labelj);
 304
 305         TLorentzVector mcpi(mcTracki->Px(),mcTracki->Py(),mcTracki->Pz(),
 306                             TMath::Sqrt(AliAnalysisMuonUtility::MuonMass2()+mcTracki->P()*mcTracki->P()));
 307         TLorentzVector mcpj(mcTrackj->Px(),mcTrackj->Py(),mcTrackj->Pz(),
 308                             TMath::Sqrt(AliAnalysisMuonUtility::MuonMass2()+mcTrackj->P()*mcTrackj->P()));
 309
 310         mcpj += mcpi;
 311
 312         inputWeightMC = WeightDistribution(mcpj.Pt(),mcpj.Rapidity());
 313
 314         Histo(eventSelection,triggerClassName,centrality,pairCutName,"PtRecVsSim")->Fill(mcpj.Pt(),pair4Momentum.Pt());
 315
 316         if ( !IsHistogramDisabled("Pt") )
 317         {
 318           MCHisto(eventSelection,triggerClassName,centrality,pairCutName,"Pt")->Fill(mcpj.Pt(),inputWeightMC);
 319         }
 320         if ( !IsHistogramDisabled("Y") )
 321         {
 322           MCHisto(eventSelection,triggerClassName,centrality,pairCutName,"Y")->Fill(mcpj.Rapidity(),inputWeightMC);
 323         }
 324         if ( !IsHistogramDisabled("Eta") )
 325         {
 326           MCHisto(eventSelection,triggerClassName,centrality,pairCutName,"Eta")->Fill(mcpj.Eta());
 327         }
 328         pair4MomentumMC = &mcpj;
 329
 330       }
 331     }
 332
 333
 334     TObjArray* bins = Binning()->CreateBinObjArray("psi","integrated,ptvsy,yvspt,pt,y,phi,nch,dnchdeta,ntrcorr,ntrcorrpt,ntrcorry,relntrcorr",""); // We may include: ,v0a,v0acent
 335     TIter nextBin(bins);
 336     AliAnalysisMuMuBinning::Range* r;
 337
 338     while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) )
 339     {
 340       Bool_t ok(kFALSE);
 341       Bool_t okMC(kFALSE);
 342
 343       if ( r->IsIntegrated() )
 344       {
 345         ok = kTRUE;
 346         if ( pair4MomentumMC ) okMC = kTRUE;
 347       }
 348       else if ( r->Is2D() )
 349       {
 350         if ( r->AsString().BeginsWith("PTVSY") )
 351         {
 352           ok = r->IsInRange(pair4Momentum.Rapidity(),pair4Momentum.Pt());
 353           if ( pair4MomentumMC ) okMC = r->IsInRange(pair4MomentumMC->Rapidity(),pair4MomentumMC->Pt());
 354         }
 355         else if ( r->AsString().BeginsWith("YVSPT") )
 356         {
 357           ok = r->IsInRange(pair4Momentum.Pt(),pair4Momentum.Rapidity());
 358           if ( pair4MomentumMC ) okMC = r->IsInRange(pair4MomentumMC->Pt(),pair4MomentumMC->Rapidity());
 359         }
 360         else if ( r->Quantity() == "NTRCORRPT" )
 361         {
 362           TList* list = static_cast<TList*>(Event()->FindListObject("NCH"));
 363           if (list)
 364           {
 365             Int_t i(-1);
 366             Bool_t parFound(kFALSE);
 367             while ( i < list->GetEntries() - 1 && !parFound )
 368             {
 369               i++;
 370               while ( list->At(i)->IsA() != TParameter<Double_t>::Class() && i < list->GetEntries() - 1 ) // In case there is a diferent object, just to skip it
 371               {
 372                 i++;
 373               }
 374
 375               TParameter<Double_t>* p = static_cast<TParameter<Double_t>*>(list->At(i));
 376
 377               if ( TString(p->GetName()).Contains("NtrCorr") )
 378               {
 379                 parFound = kTRUE;
 380                 ok = r->IsInRange(p->GetVal(),pair4Momentum.Pt());
 381               }
 382             }
 383           }
 384
 385         }
 386         else if ( r->Quantity() == "NTRCORRY" )
 387         {
 388           TList* list = static_cast<TList*>(Event()->FindListObject("NCH"));
 389           if (list)
 390           {
 391             Int_t i(-1);
 392             Bool_t parFound(kFALSE);
 393             while ( i < list->GetEntries() - 1 && !parFound )
 394             {
 395               i++;
 396               while ( list->At(i)->IsA() != TParameter<Double_t>::Class() && i < list->GetEntries() - 1 ) // In case there is a diferent object, just to skip it
 397               {
 398                 i++;
 399               }
 400
 401               TParameter<Double_t>* p = static_cast<TParameter<Double_t>*>(list->At(i));
 402
 403               if ( TString(p->GetName()).Contains("NtrCorr") )
 404               {
 405                 parFound = kTRUE;
 406                 ok = r->IsInRange(p->GetVal(),pair4Momentum.Rapidity());
 407               }
 408             }
 409           }
 410
 411         }
 412
 413         else
 414         {
 415           AliError(Form("Don't know how to deal with 2D bin %s",r->AsString().Data()));
 416         }
 417       }
 418       else
 419       {
 420         if ( r->Quantity() == "PT" )
 421         {
 422           ok = r->IsInRange(pair4Momentum.Pt());
 423           if ( pair4MomentumMC ) okMC = r->IsInRange(pair4MomentumMC->Pt());
 424         }
 425         else if ( r->Quantity() == "Y" )
 426         {
 427           ok = r->IsInRange(pair4Momentum.Rapidity());
 428           if ( pair4MomentumMC ) okMC = r->IsInRange(pair4MomentumMC->Rapidity());
 429         }
 430         else if ( r->Quantity() == "PHI" )
 431         {
 432           ok = r->IsInRange(pair4Momentum.Phi());
 433           if ( pair4MomentumMC ) okMC = r->IsInRange(pair4MomentumMC->Phi());
 434         }
 435         else if ( r->Quantity() == "DNCHDETA" )
 436         {
 437           TList* list = static_cast<TList*>(Event()->FindListObject("NCH"));
 438           if (list)
 439           {
 440             Int_t i(-1);
 441             Bool_t parFound(kFALSE);
 442             while ( i < list->GetEntries() - 1 && !parFound )
 443             {
 444               i++;
 445               while ( list->At(i)->IsA() != TParameter<Double_t>::Class() && i < list->GetEntries() - 1 ) // In case there is a diferent object, just to skip it
 446               {
 447                 i++;
 448               }
 449
 450               TParameter<Double_t>* p = static_cast<TParameter<Double_t>*>(list->At(i));
 451
 452               if ( TString(p->GetName()).Contains("dNchdEta") )
 453               {
 454                 parFound = kTRUE;
 455                 ok = r->IsInRange(p->GetVal());
 456               }
 457             }
 458           }
 459 //          else AliFatal("No dNchdEta info on Event");
 460
 461         }
 462         else if ( r->Quantity() == "NTRCORR" || r->Quantity() == "RELNTRCORR" )
 463         {
 464           TList* list = static_cast<TList*>(Event()->FindListObject("NCH"));
 465           if (list)
 466           {
 467             Int_t i(-1);
 468             Bool_t parFound(kFALSE);
 469             while ( i < list->GetEntries() - 1 && !parFound )
 470             {
 471               i++;
 472               while ( list->At(i)->IsA() != TParameter<Double_t>::Class() && i < list->GetEntries() - 1 ) // In case there is a diferent object, just to skip it
 473               {
 474                 i++;
 475               }
 476
 477               TParameter<Double_t>* p = static_cast<TParameter<Double_t>*>(list->At(i));
 478
 479               if ( TString(p->GetName()).Contains("NtrCorr") )
 480               {
 481                 parFound = kTRUE;
 482                 if ( r->Quantity() == "NTRCORR" ) ok = r->IsInRange(p->GetVal());
 483                 else ok = r->IsInRange(p->GetVal()/5.97);
 484               }
 485             }
 486           }
 487 //          else AliFatal("No ntrcorr info on Event");
 488
 489         }
 490       }
 491
 492       if ( ok || okMC )
 493       {
 494         TString minvName = GetMinvHistoName(*r,kFALSE);
 495
 496         if (!IsHistogramDisabled(minvName.Data()))
 497         {
 498           TH1* h(0x0);
 499           if ( ok )
 500           {
 501             h = Histo(eventSelection,triggerClassName,centrality,pairCutName,minvName.Data());
 502
 503             if (!h)
 504             {
 505               AliError(Form("Could not get %s",minvName.Data()));
 506               //continue;
 507             }
 508             else h->Fill(pair4Momentum.M(),inputWeight);
 509           }
 510           if( okMC )
 511           {
 512             h = MCHisto(eventSelection,triggerClassName,centrality,pairCutName,minvName.Data());
 513
 514             if (!h)
 515             {
 516               AliError(Form("Could not get MC %s",minvName.Data()));
 517               //continue;
 518             }
 519             else h->Fill(pair4MomentumMC->M(),inputWeightMC);
 520           }
 521
 522           if ( fcomputeMeanPt )
 523           {
 524             TString hprofName("");
 525
 526             if ( ok )
 527             {
 528               hprofName= Form("MeanPtVs%s",minvName.Data());
 529
 530               TProfile* hprof = Prof(eventSelection,triggerClassName,centrality,pairCutName,hprofName.Data());
 531
 532               if ( !hprof )
 533               {
 534                 AliError(Form("Could not get %s",hprofName.Data()));
 535               }
 536
 537               else
 538               {
 539                 //              hprof->Approximate(); //I dont think its necessary here
 540                 hprof->Fill(pair4Momentum.M(),pair4Momentum.Pt(),inputWeight);
 541               }
 542             }
 543             if ( okMC )
 544             {
 545               hprofName= Form("MeanPtVs%s",minvName.Data());
 546
 547               TProfile* hprof = MCProf(eventSelection,triggerClassName,centrality,pairCutName,hprofName.Data());
 548
 549               if ( !hprof )
 550               {
 551                 AliError(Form("Could not get MC %s",hprofName.Data()));
 552               }
 553
 554               else
 555               {
 556                 //              hprof->Approximate(); //I dont think its necessary here
 557                 hprof->Fill(pair4MomentumMC->M(),pair4MomentumMC->Pt(),inputWeightMC);
 558               }
 559             }
 560
 561           }
 562
 563         }
 564
 565         if ( ShouldCorrectDimuonForAccEff() )
 566         {
 567
 568           Double_t AccxEff(0);
 569           Bool_t okAccEff(kFALSE);
 570           if ( ok )
 571           {
 572             AccxEff = GetAccxEff(pair4Momentum.Pt(),pair4Momentum.Rapidity());
 573             if ( AccxEff <= 0.0 )
 574             {
 575               AliError(Form("AccxEff < 0 for pt = %f & y = %f ",pair4Momentum.Pt(),pair4Momentum.Rapidity()));
 576               //            continue;
 577             }
 578             else okAccEff = kTRUE;
 579           }
 580
 581           Double_t AccxEffMC(0);
 582           Bool_t okAccEffMC(kFALSE);
 583           if ( okMC )
 584           {
 585              AccxEffMC= GetAccxEff(pair4MomentumMC->Pt(),pair4MomentumMC->Rapidity());
 586             if ( AccxEffMC <= 0.0 )
 587             {
 588               AliError(Form("AccxEff < 0 for MC pair with pt = %f & y = %f ",pair4MomentumMC->Pt(),pair4MomentumMC->Rapidity()));
 589               //            continue;
 590             }
 591             else okAccEffMC = kTRUE;
 592           }
 593
 594           minvName = GetMinvHistoName(*r,kTRUE);
 595
 596           if (!IsHistogramDisabled(minvName.Data()))
 597           {
 598             TH1* hCorr = Histo(eventSelection,triggerClassName,centrality,pairCutName,minvName.Data());
 599
 600             if (!hCorr)
 601             {
 602               AliError(Form("Could not get %sr",minvName.Data()));
 603             }
 604
 605             else  if ( okAccEff ) hCorr->Fill(pair4Momentum.M(),inputWeight/AccxEff);
 606
 607             if( okAccEffMC )
 608             {
 609               hCorr = MCHisto(eventSelection,triggerClassName,centrality,pairCutName,minvName.Data());
 610
 611               if (!hCorr)
 612               {
 613                 AliError(Form("Could not get MC %s",minvName.Data()));
 614                 //continue;
 615               }
 616               else hCorr->Fill(pair4MomentumMC->M(),inputWeightMC/AccxEffMC);
 617             }
 618
 619             if ( fcomputeMeanPt )
 620             {
 621               TString hprofCorrName("");
 622               if( ok )
 623               {
 624                 hprofCorrName = Form("MeanPtVs%s",minvName.Data());
 625
 626                 TProfile* hprofCorr = Prof(eventSelection,triggerClassName,centrality,pairCutName,hprofCorrName.Data());
 627
 628                 if ( !hprofCorr )
 629                 {
 630                   AliError(Form("Could not get %s",hprofCorrName.Data()));
 631                 }
 632                 else if ( okAccEff )
 633                 {
 634                   //                hprofCorr->Approximate(); //I dont know if its necessary here
 635                   hprofCorr->Fill(pair4Momentum.M(),pair4Momentum.Pt(),inputWeight/AccxEff);
 636                 }
 637               }
 638               if( okMC )
 639               {
 640                 hprofCorrName = Form("MeanPtVs%s",minvName.Data());
 641
 642                 TProfile* hprofCorr = MCProf(eventSelection,triggerClassName,centrality,pairCutName,hprofCorrName.Data());
 643
 644                 if ( !hprofCorr )
 645                 {
 646                   AliError(Form("Could not get MC %s",hprofCorrName.Data()));
 647                 }
 648                 else if ( okAccEffMC )
 649                 {
 650                   //                hprofCorr->Approximate(); //I dont know if its necessary here
 651                   hprofCorr->Fill(pair4MomentumMC->M(),pair4MomentumMC->Pt(),inputWeightMC/AccxEffMC);
 652                 }
 653               }
 654
 655             }
 656
 657           }
 658
 659         }
 660       }
 661     }
 662
 663     delete bins;
 664   }
 665
 666 }
 667
 668
 669 //_____________________________________________________________________________
 670 void AliAnalysisMuMuMinv::FillHistosForMCEvent(const char* eventSelection,const char* triggerClassName,const char* centrality)
 671 {
 672   // Fill the input Monte-Carlo histograms related to muons. Intended to be used on pure simulations so we wont use eventSelection, triggerClassName and centrality variables.
 673
 674   if ( !HasMC() ) return;
 675
 676   // Specific things for MC // These histos should go in the AliAnalysisMuMuGlobal task, but then we have to loop 2 times on the MCTracks...
 677 //  if (!Histo(MCInputPrefix(),"ALL","Pt"))
 678 //  {
 679 //    HistogramCollection()->Adopt(Form("/%s/ALL",MCInputPrefix()),new TH1F("Pt","Pt",200,0,25));
 680 //    HistogramCollection()->Adopt(Form("/%s/INYRANGE",MCInputPrefix()),new TH1F("Pt","Pt",200,0,25));
 681 //
 682 //    Double_t rapidityMin = -5;
 683 //    Double_t rapidityMax = -2;
 684 //    Int_t nbinsRapidity = GetNbins(rapidityMin,rapidityMax,0.05);
 685 //
 686 //    HistogramCollection()->Adopt(Form("/%s/ALL",MCInputPrefix()),new TH1F("Y","Y",nbinsRapidity,rapidityMin,rapidityMax));
 687 //    HistogramCollection()->Adopt(Form("/%s/INYRANGE",MCInputPrefix()),new TH1F("Y","Y",nbinsRapidity,rapidityMin,rapidityMax));
 688 //
 689 //    Double_t etaMin = -5;
 690 //    Double_t etaMax = -2;
 691 //    Int_t nbinsEta = GetNbins(etaMin,etaMax,0.05);
 692 //
 693 //    HistogramCollection()->Adopt(Form("/%s/ALL",MCInputPrefix()),new TH1F("Eta","Eta",nbinsEta,etaMin,etaMax));
 694 //    HistogramCollection()->Adopt(Form("/%s/INYRANGE",MCInputPrefix()),new TH1F("Eta","Eta",nbinsEta,etaMin,etaMax));
 695 //  }
 696
 697   Int_t nMCTracks = MCEvent()->GetNumberOfTracks();
 698
 699   TObjArray* bins = Binning()->CreateBinObjArray("psi","integrated,ptvsy,yvspt,pt,y,phi,nch,dnchdeta,ntrcorr,ntrcorrpt,ntrcorry,relntrcorr","");//We may include: ,v0a,v0acent.
 700   TIter nextBin(bins);
 701   AliAnalysisMuMuBinning::Range* r;
 702
 703   for ( Int_t i = 0; i < nMCTracks; ++i )
 704   {
 705     AliVParticle* part = MCEvent()->GetTrack(i);
 706
 707     if  (AliAnalysisMuonUtility::IsPrimary(part,MCEvent()) &&
 708          AliAnalysisMuonUtility::GetMotherIndex(part)==-1)
 709     {
 710       Double_t inputWeight = WeightDistribution(part->Pt(),part->Y());
 711
 712       MCHisto(eventSelection,triggerClassName,centrality,"Pt")->Fill(part->Pt(),inputWeight);
 713       MCHisto(eventSelection,triggerClassName,centrality,"Y")->Fill(part->Y(),inputWeight);
 714       MCHisto(eventSelection,triggerClassName,centrality,"Eta")->Fill(part->Eta());
 715
 716       if ( part->Y() < -2.5 && part->Y() > -4.0 )
 717       {
 718         MCHisto(eventSelection,triggerClassName,Form("%s/INYRANGE",centrality),"Pt")->Fill(part->Pt(),inputWeight);
 719         MCHisto(eventSelection,triggerClassName,Form("%s/INYRANGE",centrality),"Y")->Fill(part->Y(),inputWeight);
 720         MCHisto(eventSelection,triggerClassName,Form("%s/INYRANGE",centrality),"Eta")->Fill(part->Eta());
 721       }
 722
 723       nextBin.Reset();
 724
 725       while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) )
 726       {
 727         Bool_t ok(kFALSE);
 728
 729         if ( r->IsIntegrated() )
 730         {
 731           ok = kTRUE;
 732         }
 733         else if ( r->Is2D() )
 734         {
 735           if ( r->AsString().BeginsWith("PTVSY") )
 736           {
 737             ok = r->IsInRange(part->Y(),part->Pt());
 738           }
 739           else if ( r->AsString().BeginsWith("YVSPT") )
 740           {
 741             ok = r->IsInRange(part->Pt(),part->Y());
 742           }
 743           else
 744           {
 745             AliError(Form("Don't know how to deal with 2D bin %s",r->AsString().Data()));
 746           }
 747         }
 748         else
 749         {
 750           if ( r->Quantity() == "PT" )
 751           {
 752             ok = r->IsInRange(part->Pt());
 753           }
 754           else if ( r->Quantity() == "Y" )
 755           {
 756             ok = r->IsInRange(part->Y());
 757           }
 758           else if ( r->Quantity() == "PHI" )
 759           {
 760             ok = r->IsInRange(part->Phi());
 761           }
 762         }
 763
 764         if ( ok )
 765         {
 766           TString hname = GetMinvHistoName(*r,kFALSE);
 767
 768           if (!IsHistogramDisabled(hname.Data()))
 769           {
 770
 771             TH1* h = MCHisto(eventSelection,triggerClassName,centrality,hname.Data());
 772
 773             if (!h)
 774             {
 775               AliError(Form("Could not get /%s/%s/%s/%s/ %s",MCInputPrefix(),eventSelection,triggerClassName,centrality,hname.Data()));
 776               continue;
 777             }
 778
 779             h->Fill(part->M(),inputWeight);
 780
 781             if ( part->Y() < -2.5 && part->Y() > -4.0 )
 782             {
 783               h = MCHisto(eventSelection,triggerClassName,Form("%s/INYRANGE",centrality),hname.Data());
 784               if (!h)
 785               {
 786                 AliError(Form("Could not get /%s/%s/%s/%s/INYRANGE %s",MCInputPrefix(),eventSelection,triggerClassName,centrality,hname.Data()));
 787                 continue;
 788               }
 789               h->Fill(part->M(),inputWeight);
 790             }
 791
 792           }
 793
 794           if ( fcomputeMeanPt )
 795           {
 796             TString hprofName= Form("MeanPtVs%s",hname.Data());
 797
 798             TProfile* hprof = MCProf(eventSelection,triggerClassName,centrality,hprofName.Data());
 799
 800             if ( !hprof )
 801             {
 802               AliError(Form("Could not get %s",hprofName.Data()));
 803             }
 804             else
 805             {
 806               //              hprof->Approximate(); //I dont think its necessary here
 807               hprof->Fill(part->M(),part->Pt(),inputWeight);
 808             }
 809
 810             if ( part->Y() < -2.5 && part->Y() > -4.0 )
 811             {
 812               hprof = MCProf(eventSelection,triggerClassName,Form("%s/INYRANGE",centrality),hprofName.Data());
 813
 814               if ( !hprof )
 815               {
 816                 AliError(Form("Could not get %s",hprofName.Data()));
 817               }
 818               else
 819               {
 820                 //              hprof->Approximate(); //I dont think its necessary here
 821                 hprof->Fill(part->M(),part->Pt(),inputWeight);
 822               }
 823
 824             }
 825
 826
 827           }
 828
 829 //          if ( ShouldCorrectDimuonForAccEff() ) //What is the sense of this?? We should not correct the input
 830 //          {
 831 //            Double_t AccxEff = GetAccxEff(part->Pt(),part->Y());
 832 //            if ( AccxEff <= 0.0 )
 833 //            {
 834 //              AliError(Form("AccxEff < 0 for pt = %f & y = %f ",part->Pt(),part->Y()));
 835 //              continue;
 836 //            }
 837 //            hname = GetMinvHistoName(*r,kTRUE);
 838 //
 839 //            if (!IsHistogramDisabled(hname.Data()))
 840 //            {
 841 //
 842 //              TH1* h = MCHisto(eventSelection,triggerClassName,centrality,hname.Data());
 843 //
 844 //              if (!h)
 845 //              {
 846 //                AliError(Form("Could not get /%s/%s/%s/%s/ %s",MCInputPrefix(),eventSelection,triggerClassName,centrality,hname.Data()));
 847 //                continue;
 848 //              }
 849 //
 850 //              h->Fill(part->M(),inputWeight/AccxEff);
 851 //
 852 //              if ( part->Y() < -2.5 && part->Y() > -4.0 )
 853 //              {
 854 //                h = MCHisto(eventSelection,triggerClassName,Form("%s/INYRANGE",centrality),hname.Data());
 855 //                if (!h)
 856 //                {
 857 //                  AliError(Form("Could not get /%s/%s/%s/%s/INYRANGE %s",MCInputPrefix(),eventSelection,triggerClassName,Form("%s/INYRANGE",centrality),hname.Data()));
 858 //                  continue;
 859 //                }
 860 //                h->Fill(part->M(),inputWeight./AccxEff);
 861 //              }
 862 //
 863 //            }
 864 //
 865 //          }
 866
 867         }
 868       }
 869     }
 870   }
 871
 872   delete bins;
 873
 874
 875 //  for ( Int_t i = 0; i < nMCTracks; ++i )
 876 //  {
 877 //    AliVParticle* part = MCEvent()->GetTrack(i);
 878 //
 879 //    if  (AliAnalysisMuonUtility::IsPrimary(part,MCEvent()) &&
 880 //         AliAnalysisMuonUtility::GetMotherIndex(part)==-1)
 881 //    {
 882 //
 883 //      Histo(MCInputPrefix(),"ALL","Pt")->Fill(part->Pt());
 884 //      Histo(MCInputPrefix(),"ALL","Y")->Fill(part->Y());
 885 //      Histo(MCInputPrefix(),"ALL","Eta")->Fill(part->Eta());
 886 //
 887 //      if ( part->Y() < -2.5 && part->Y() > -4.0 )
 888 //      {
 889 //        Histo(MCInputPrefix(),"INYRANGE","Pt")->Fill(part->Pt());
 890 //        Histo(MCInputPrefix(),"INYRANGE","Y")->Fill(part->Y());
 891 //        Histo(MCInputPrefix(),"INYRANGE","Eta")->Fill(part->Eta());
 892 //      }
 893 //
 894 //      nextBin.Reset();
 895 //
 896 //      while ( ( r = static_cast<AliAnalysisMuMuBinning::Range*>(nextBin()) ) )
 897 //      {
 898 //        Bool_t ok(kFALSE);
 899 //
 900 //        if ( r->IsIntegrated() )
 901 //        {
 902 //          ok = kTRUE;
 903 //        }
 904 //        else if ( r->Is2D() )
 905 //        {
 906 //          if ( r->AsString().BeginsWith("PTVSY") )
 907 //          {
 908 //            ok = r->IsInRange(part->Y(),part->Pt());
 909 //          }
 910 //          else if ( r->AsString().BeginsWith("YVSPT") )
 911 //          {
 912 //            ok = r->IsInRange(part->Pt(),part->Y());
 913 //          }
 914 //          else
 915 //          {
 916 //            AliError(Form("Don't know how to deal with 2D bin %s",r->AsString().Data()));
 917 //          }
 918 //        }
 919 //        else
 920 //        {
 921 //          if ( r->Quantity() == "PT" )
 922 //          {
 923 //            ok = r->IsInRange(part->Pt());
 924 //          }
 925 //          else if ( r->Quantity() == "Y" )
 926 //          {
 927 //            ok = r->IsInRange(part->Y());
 928 //          }
 929 //          else if ( r->Quantity() == "PHI" )
 930 //          {
 931 //            ok = r->IsInRange(part->Phi());
 932 //          }
 933 //        }
 934 //
 935 //        if ( ok )
 936 //        {
 937 //          TString hname = GetMinvHistoName(*r,kFALSE);
 938 //
 939 //          if (!IsHistogramDisabled(hname.Data()))
 940 //          {
 941 //
 942 //            TH1* h = Histo(MCInputPrefix(),"ALL",hname.Data());
 943 //
 944 //            if (!h)
 945 //            {
 946 //              AliError(Form("Could not get ALL %s",hname.Data()));
 947 //              continue;
 948 //            }
 949 //
 950 //            h->Fill(part->M());
 951 //
 952 //            if ( part->Y() < -2.5 && part->Y() > -4.0 )
 953 //            {
 954 //              h = Histo(MCInputPrefix(),"INYRANGE",hname.Data());
 955 //              if (!h)
 956 //              {
 957 //                AliError(Form("Could not get INYRANGE %s",hname.Data()));
 958 //                continue;
 959 //              }
 960 //              h->Fill(part->M());
 961 //            }
 962 //
 963 //          }
 964 //
 965 //          if ( ShouldCorrectDimuonForAccEff() )
 966 //          {
 967 //            Double_t AccxEff = GetAccxEff(part->Pt(),part->Y());
 968 //            if ( AccxEff <= 0.0 )
 969 //            {
 970 //              AliError(Form("AccxEff < 0 for pt = %f & y = %f ",part->Pt(),part->Y()));
 971 //              continue;
 972 //            }
 973 //            hname = GetMinvHistoName(*r,kTRUE);
 974 //
 975 //            if (!IsHistogramDisabled(hname.Data()))
 976 //            {
 977 //
 978 //              TH1* h = Histo(MCInputPrefix(),"ALL",hname.Data());
 979 //
 980 //              if (!h)
 981 //              {
 982 //                AliError(Form("Could not get ALL %s",hname.Data()));
 983 //                continue;
 984 //              }
 985 //
 986 //              h->Fill(part->M(),1./AccxEff);
 987 //
 988 //              if ( part->Y() < -2.5 && part->Y() > -4.0 )
 989 //              {
 990 //                h = Histo(MCInputPrefix(),"INYRANGE",hname.Data());
 991 //                if (!h)
 992 //                {
 993 //                  AliError(Form("Could not get INYRANGE %s",hname.Data()));
 994 //                  continue;
 995 //                }
 996 //                h->Fill(part->M(),1./AccxEff);
 997 //              }
 998 //
 999 //            }
1000 //
1001 //          }
1002 //
1003 //        }
1004 //      }
1005 //    }
1006 //
1007 //    delete bins;
1008 //  }
1009 }
1010
1011 //_____________________________________________________________________________
1012 TString AliAnalysisMuMuMinv::GetMinvHistoName(const AliAnalysisMuMuBinning::Range& r, Bool_t accEffCorrected) const
1013 {
1014   return TString::Format("MinvUS%s%s",r.AsString().Data(),
1015                          accEffCorrected ? "_AccEffCorr" : "");
1016 }
1017
1018
1019 //_____________________________________________________________________________
1020 Double_t AliAnalysisMuMuMinv::GetAccxEff(Double_t pt,Double_t rapidity)
1021 {
1022   if (!fAccEffHisto)
1023   {
1024     AliError("ERROR: No AccxEff histo");
1025     return 0;
1026   }
1027   Int_t bin = fAccEffHisto->FindBin(pt,rapidity);
1028   Double_t accXeff = fAccEffHisto->GetBinContent(bin);
1029
1030   return accXeff;
1031 }
1032
1033 //_____________________________________________________________________________
1034 Double_t AliAnalysisMuMuMinv::WeightDistribution(Double_t pt,Double_t rapidity)
1035 {
1036   //Return a weight for a dimuon pt and y, which depend on the varied distributions.
1037   // FIXME: hard coded, find a clean way to fix the distribution parameters from outside
1038
1039   if (!HasMC() ) return 1.;
1040
1041   //================ p-Pb ==============//
1042   //value for input distribution: this is the nominal pt and y distribution
1043   Double_t paryPPB[2] = {1.0,0.174};
1044   Double_t parptPPB[3] = {1.0,0.0557,3.52};
1045
1046   Double_t paryHardPPB = 0.1344, parySoftPPB = 0.1971;
1047   Double_t par1ptHardPPB = 5.51e-2, par2ptHardPPB = 3.47,
1048   par1ptSoftPPB = 5.67e-2, par2ptSoftPPB = 3.68;
1049   //systematic 1: hardest in y  x softest in pt
1050   Double_t pary1[2] = {1.0,paryHardPPB};
1051   Double_t parpt1[3] = {1.0,par1ptSoftPPB,par2ptSoftPPB};
1052   //systematic 2: hardest in y x hardest in pt
1053   Double_t pary2[2] = {1.0,paryHardPPB};
1054   Double_t parpt2[3] = {1.0,par1ptHardPPB,par2ptHardPPB};
1055   //systematic 3: softest in y  x softest in pt
1056   Double_t pary3[2] = {1.0, parySoftPPB};
1057   Double_t parpt3[3] = {1.0,par1ptSoftPPB,par2ptSoftPPB};
1058   //systematic 4: softest in y  x hardest in pt
1059   Double_t pary4[2] = {1.0,parySoftPPB};
1060   Double_t parpt4[3] = {1.0,par1ptHardPPB,par2ptHardPPB};
1061
1062   Double_t funcptvalPPB = powerLaw3Par(&pt,parptPPB);
1063   Double_t funcyvalPPB = normPol12Par(&rapidity,paryPPB);
1064
1065   //================ Pb-p ==============//
1066   //value for input distribution
1067   Double_t paryPBP[2] = {1.0,0.189};
1068   Double_t parptPBP[3] = {1.0,0.0592,3.92};
1069
1070   Double_t paryHardPBP = 0.1517, parySoftPBP = 0.2191;
1071   Double_t par1ptHardPBP = 5.58e-2, par2ptHardPBP = 3.83,
1072   par1ptSoftPBP = 5.59e-2, par2ptSoftPBP = 4.31;
1073   //systematic 5: hardest in y  x softest in pt
1074   Double_t pary5[2] = {1.0,paryHardPBP};
1075   Double_t parpt5[3] = {1.0,par1ptSoftPBP,par2ptSoftPBP};
1076   //systematic 6: hardest in y x hardest in pt
1077   Double_t pary6[2] = {1.0,paryHardPBP};
1078   Double_t parpt6[3] = {1.0,par1ptHardPBP,par2ptHardPBP};
1079   //systematic 7: softest in y  x softest in pt
1080   Double_t pary7[2] = {1.0, parySoftPBP};
1081   Double_t parpt7[3] = {1.0,par1ptSoftPBP,par2ptSoftPBP};
1082   //systematic 8: softest in y  x hardest in pt
1083   Double_t pary8[2] = {1.0,parySoftPBP};
1084   Double_t parpt8[3] = {1.0,par1ptHardPBP,par2ptHardPBP};
1085
1086   Double_t funcptvalPBP = powerLaw3Par(&pt,parptPBP);
1087   Double_t funcyvalPBP = normPol12Par(&rapidity,paryPBP);
1088
1089   //================ pp ==============//
1090   //value for input distribution
1091   Double_t paryPP[2] = {3.0,0.514/3.};
1092   Double_t parptPP[3] = {1.0,0.0546,3.90};
1093
1094   Double_t paryHardPP = 0.4125/3., parySoftPP = 0.5958/3.;
1095   Double_t par1ptHardPP = 4.78e-2, par2ptHardPP = 3.65,//4.84e-2/3.45
1096   par1ptSoftPP = 6.12e-2, par2ptSoftPP = 4.31;//5.47e-2//4.29
1097   //systematic 9: hardest in y  x softest in pt
1098   Double_t pary9[2] = {3.0,paryHardPP};
1099   Double_t parpt9[3] = {1.0,par1ptSoftPP,par2ptSoftPP};
1100   //systematic 10: hardest in y x hardest in pt
1101   Double_t pary10[2] = {3.0,paryHardPP};
1102   Double_t parpt10[3] = {1.0,par1ptHardPP,par2ptHardPP};
1103   //systematic 11: softest in y  x softest in pt
1104   Double_t pary11[2] = {3.0, parySoftPP};
1105   Double_t parpt11[3] = {1.0,par1ptSoftPP,par2ptSoftPP};
1106   //systematic 12: softest in y  x hardest in pt
1107   Double_t pary12[2] = {3.0,parySoftPP};
1108   Double_t parpt12[3] = {1.0,par1ptHardPP,par2ptHardPP};
1109
1110   Double_t funcptvalPP = powerLaw3Par(&pt,parptPP);
1111   Double_t funcyvalPP = normPol12Par(&rapidity,paryPP);
1112
1113   //______
1114   Double_t weight(1.),funcptsyst(0.),funcysyst(0.);
1115   switch ( fsystLevel )
1116   {
1117     case 0:
1118       weight = 1;
1119       break;
1120     case 1:
1121       funcptsyst = powerLaw3Par(&pt,parpt1);
1122       if ( funcptvalPPB > 0 && funcptsyst > 0 ) weight = funcptsyst/funcptvalPPB;
1123       else  weight = 1;
1124       funcysyst = normPol12Par(&rapidity,pary1);
1125       if ( funcyvalPPB > 0 && funcysyst > 0 ) weight *= funcysyst/funcyvalPPB;
1126       break;
1127     case 2:
1128       funcptsyst = powerLaw3Par(&pt,parpt2);
1129       if ( funcptvalPPB > 0 && funcptsyst > 0 ) weight = funcptsyst/funcptvalPPB;
1130       else  weight = 1;
1131       funcysyst = normPol12Par(&rapidity,pary2);
1132       if ( funcyvalPPB > 0 && funcysyst > 0 ) weight *= funcysyst/funcyvalPPB;
1133       break;
1134     case 3:
1135       funcptsyst = powerLaw3Par(&pt,parpt3);
1136       if ( funcptvalPPB > 0 && funcptsyst > 0 ) weight = funcptsyst/funcptvalPPB;
1137       else  weight = 1;
1138       funcysyst = normPol12Par(&rapidity,pary3);
1139       if ( funcyvalPPB > 0 && funcysyst > 0 ) weight *= funcysyst/funcyvalPPB;
1140       break;
1141     case 4:
1142       funcptsyst = powerLaw3Par(&pt,parpt4);
1143       if ( funcptvalPPB > 0 && funcptsyst > 0 ) weight = funcptsyst/funcptvalPPB;
1144       else  weight = 1;
1145       funcysyst = normPol12Par(&rapidity,pary4);
1146       if ( funcyvalPPB > 0 && funcysyst > 0 ) weight *= funcysyst/funcyvalPPB;
1147       break;
1148     case 5:
1149       funcptsyst = powerLaw3Par(&pt,parpt5);
1150       if ( funcptvalPBP > 0 && funcptsyst > 0 ) weight = funcptsyst/funcptvalPBP;
1151       else  weight = 1;
1152       funcysyst = normPol12Par(&rapidity,pary5);
1153       if ( funcyvalPBP > 0 && funcysyst > 0 ) weight *= funcysyst/funcyvalPBP;
1154       break;
1155     case 6:
1156       funcptsyst = powerLaw3Par(&pt,parpt6);
1157       if ( funcptvalPBP > 0 && funcptsyst > 0 ) weight = funcptsyst/funcptvalPBP;
1158       else  weight = 1;
1159       funcysyst = normPol12Par(&rapidity,pary6);
1160       if ( funcyvalPBP > 0 && funcysyst > 0 ) weight *= funcysyst/funcyvalPBP;
1161       break;
1162     case 7:
1163       funcptsyst = powerLaw3Par(&pt,parpt7);
1164       if ( funcptvalPBP > 0 && funcptsyst > 0 ) weight = funcptsyst/funcptvalPBP;
1165       else  weight = 1;
1166       funcysyst = normPol12Par(&rapidity,pary7);
1167       if ( funcyvalPBP > 0 && funcysyst > 0 ) weight *= funcysyst/funcyvalPBP;
1168       break;
1169     case 8:
1170       funcptsyst = powerLaw3Par(&pt,parpt8);
1171       if ( funcptvalPBP > 0 && funcptsyst > 0 ) weight = funcptsyst/funcptvalPBP;
1172       else  weight = 1;
1173       funcysyst = normPol12Par(&rapidity,pary8);
1174       if ( funcyvalPBP > 0 && funcysyst > 0 ) weight *= funcysyst/funcyvalPBP;
1175       break;
1176     case 9:
1177       funcptsyst = powerLaw3Par(&pt,parpt9);
1178       if ( funcptvalPP > 0 && funcptsyst > 0 ) weight = funcptsyst/funcptvalPP;
1179       else  weight = 1;
1180       funcysyst = normPol12Par(&rapidity,pary9);
1181       if ( funcyvalPP > 0 && funcysyst > 0 ) weight *= funcysyst/funcyvalPP;
1182       break;
1183     case 10:
1184       funcptsyst = powerLaw3Par(&pt,parpt10);
1185       if ( funcptvalPP > 0 && funcptsyst > 0 ) weight = funcptsyst/funcptvalPP;
1186       else  weight = 1;
1187       funcysyst = normPol12Par(&rapidity,pary10);
1188       if ( funcyvalPP > 0 && funcysyst > 0 ) weight *= funcysyst/funcyvalPP;
1189       break;
1190     case 11:
1191       funcptsyst = powerLaw3Par(&pt,parpt11);
1192       if ( funcptvalPP > 0 && funcptsyst > 0 ) weight = funcptsyst/funcptvalPP;
1193       else  weight = 1;
1194       funcysyst = normPol12Par(&rapidity,pary11);
1195       if ( funcyvalPP > 0 && funcysyst > 0 ) weight *= funcysyst/funcyvalPP;
1196       break;
1197     case 12:
1198       funcptsyst = powerLaw3Par(&pt,parpt12);
1199       if ( funcptvalPP > 0 && funcptsyst > 0 ) weight = funcptsyst/funcptvalPP;
1200       else  weight = 1;
1201       funcysyst = normPol12Par(&rapidity,pary12);
1202       if ( funcyvalPP > 0 && funcysyst > 0 ) weight *= funcysyst/funcyvalPP;
1203       break;
1204   }
1205
1206   return weight;
1207 }
1208
1209 //________________________________________________________________________
1210 Double_t AliAnalysisMuMuMinv::powerLaw3Par(Double_t *x, Double_t *par)
1211 {
1212   //3 parameters
1213   Double_t arg = 0;
1214
1215   arg = par[0]*x[0] / TMath::Power( 1 + par[1]*x[0]*x[0], par[2]);
1216
1217   return arg;
1218 }
1219
1220
1221 //________________________________________________________________________
1222 Double_t AliAnalysisMuMuMinv::normPol12Par(Double_t *x, Double_t *par)
1223 {
1224   //2 parameters
1225   Double_t arg1 = 0;
1226
1227   arg1 = par[0] * ( 1 + par[1]*x[0] );
1228
1229
1230   return arg1;
1231 }
1232
1233 //_____________________________________________________________________________
1234 Bool_t AliAnalysisMuMuMinv::IsPtInRange(const AliVParticle& t1, const AliVParticle& t2, Double_t& ptmin, Double_t& ptmax) const
1235 {
1236   /// Whether the pair passes the rapidity cut
1237
1238   TLorentzVector p1(t1.Px(),t1.Py(),t1.Pz(),TMath::Sqrt(AliAnalysisMuonUtility::MuonMass2()+t1.P()*t1.P()));
1239   TLorentzVector p2(t2.Px(),t2.Py(),t2.Pz(),TMath::Sqrt(AliAnalysisMuonUtility::MuonMass2()+t2.P()*t2.P()));
1240
1241   TLorentzVector total(p1+p2);
1242
1243   Double_t pt = total.Pt();
1244
1245   return  ( pt < ptmax && pt > ptmin );
1246 }
1247
1248 //_____________________________________________________________________________
1249 void AliAnalysisMuMuMinv::NameOfIsPtInRange(TString& name, Double_t& ptmin, Double_t& ptmax) const
1250 {
1251   name.Form("PAIRPTIN%2.1f-%2.1f",ptmin,ptmax);
1252 }
1253
1254 //_____________________________________________________________________________
1255 Bool_t AliAnalysisMuMuMinv::IsRapidityInRange(const AliVParticle& t1, const AliVParticle& t2) const
1256 {
1257   TLorentzVector p1(t1.Px(),t1.Py(),t1.Pz(),TMath::Sqrt(AliAnalysisMuonUtility::MuonMass2()+t1.P()*t1.P()));
1258   TLorentzVector p2(t2.Px(),t2.Py(),t2.Pz(),TMath::Sqrt(AliAnalysisMuonUtility::MuonMass2()+t2.P()*t2.P()));
1259
1260   TLorentzVector total(p1+p2);
1261
1262   Double_t y = total.Rapidity();
1263
1264   return  ( y < -2.5 && y > -4.0 );
1265 }