PWGPP/MUON/lite/AliAnalysisTaskMuonCuts.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id: AliAnalysisTaskMuonCuts.cxx 47782 2011-02-24 18:37:31Z martinez $ */
  17
  18 /// Task for tuning cuts for single muons in the spectrometer.
  19 /// The output is a list of histograms.
  20 /// Task need to be run in steps:
  21 /// Step 1: find <DCAx> and <DCAy> in MC
  22 /// Step 2: find signal sigma_pDCA in MC
  23 /// Step 3: find <DCAx> and <DCAy> in data
  24 /// Step 4: apply sigma_pDCA cut on simulations and check the
  25 ///         resulting kinematic variables
  26 ///
  27 /// \author Diego Stocco
  28
  29 #define AliAnalysisTaskMuonCuts_cxx
  30
  31 #include "AliAnalysisTaskMuonCuts.h"
  32
  33 // ROOT includes
  34 #include "TROOT.h"
  35 #include "TH1.h"
  36 #include "TH2.h"
  37 #include "TH3.h"
  38 #include "TAxis.h"
  39 #include "TCanvas.h"
  40 #include "TLegend.h"
  41 #include "TMath.h"
  42 #include "TObjString.h"
  43 #include "TObjArray.h"
  44 #include "TF1.h"
  45 #include "TStyle.h"
  46 //#include "TMCProcess.h"
  47
  48 // STEER includes
  49 #include "AliAODEvent.h"
  50 #include "AliAODTrack.h"
  51 #include "AliAODMCParticle.h"
  52 #include "AliMCEvent.h"
  53 #include "AliMCParticle.h"
  54 //#include "AliStack.h"
  55 #include "AliESDEvent.h"
  56 #include "AliESDMuonTrack.h"
  57 #include "AliInputEventHandler.h"
  58
  59 #include "AliAnalysisManager.h"
  60
  61 // PWG includes
  62 #include "AliMergeableCollection.h"
  63 #include "AliMuonTrackCuts.h"
  64
  65
  66 /// \cond CLASSIMP
  67 ClassImp(AliAnalysisTaskMuonCuts) // Class implementation in ROOT context
  68 /// \endcond
  69
  70
  71 //________________________________________________________________________
  72 AliAnalysisTaskMuonCuts::AliAnalysisTaskMuonCuts() :
  73   AliVAnalysisMuon(),
  74   fHistoTypeKeys(0x0),
  75   fThetaAbsKeys(0x0),
  76   fSigmaCuts(TArrayD())
  77 {
  78   /// Default ctor.
  79 }
  80
  81 //________________________________________________________________________
  82 AliAnalysisTaskMuonCuts::AliAnalysisTaskMuonCuts(const char *name, const AliMuonTrackCuts& cuts) :
  83   AliVAnalysisMuon(name, cuts),
  84   fHistoTypeKeys(0x0),
  85   fThetaAbsKeys(0x0),
  86   fSigmaCuts(TArrayD())
  87 {
  88   //
  89   /// Constructor.
  90   //
  91   TString histoTypeKeys = "hDCAxVsP hDCAyVsP hPdcaVsP hPDCAVsPCheck hDCAVsPCheck hChiProbVsP hSigmaVsPt hSigmaVsEta";
  92   fHistoTypeKeys = histoTypeKeys.Tokenize(" ");
  93
  94   TString thetaAbsKeys = "ThetaAbs23 ThetaAbs310";
  95   fThetaAbsKeys = thetaAbsKeys.Tokenize(" ");
  96
  97   SetSigmaCuts();
  98 }
  99
 100
 101 //________________________________________________________________________
 102 AliAnalysisTaskMuonCuts::~AliAnalysisTaskMuonCuts()
 103 {
 104   //
 105   /// Destructor
 106   //
 107
 108   delete fHistoTypeKeys;
 109   delete fThetaAbsKeys;
 110 }
 111
 112 //___________________________________________________________________________
 113 void AliAnalysisTaskMuonCuts::MyUserCreateOutputObjects()
 114 {
 115   TH1* histo = 0x0;
 116   TString histoName = "";
 117   for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
 118     for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
 119       histoName = GetHistoName(kDCAxVsP, itheta, isrc);
 120       histo = new TH2F(histoName.Data(), histoName.Data(), 50, 0., 200., 100, -50., 50.);
 121       histo->SetXTitle("p (GeV/c)");
 122       histo->SetYTitle("DCA_{x} (cm)");
 123       AddObjectToCollection(histo);
 124
 125       histoName = GetHistoName(kDCAyVsP, itheta, isrc);
 126       histo = new TH2F(histoName.Data(), histoName.Data(), 50, 0., 200., 100, -50., 50.);
 127       histo->SetXTitle("p (GeV/c)");
 128       histo->SetYTitle("DCA_{y} (cm)");
 129       AddObjectToCollection(histo);
 130
 131       histoName = GetHistoName(kPdcaVsP, itheta, isrc);
 132       histo = new TH2F(histoName.Data(), histoName.Data(), 100, 0., 400., 50, 0., 500.);
 133       histo->SetXTitle("p (GeV/c)");
 134       histo->SetYTitle("p#timesDCA (cm #times GeV/c)");
 135       AddObjectToCollection(histo);
 136
 137       histoName = GetHistoName(kPDCAVsPCheck, itheta, isrc);
 138       histo = new TH2F(histoName.Data(), histoName.Data(), 100, 0., 800., 100, 0., 5000.);
 139       histo->SetXTitle("p (GeV/c)");
 140       histo->SetYTitle("p#timesDCA (cm #times GeV/c)");
 141       AddObjectToCollection(histo);
 142
 143       histoName = GetHistoName(kDCAVsPCheck, itheta, isrc);
 144       histo = new TH2F(histoName.Data(), histoName.Data(), 100, 0., 800., 100, 0., 200.);
 145       histo->SetXTitle("p (GeV/c)");
 146       histo->SetYTitle("DCA (cm)");
 147       AddObjectToCollection(histo);
 148
 149       histoName = GetHistoName(kChiProbVsP, itheta, isrc);
 150       histo = new TH2F(histoName.Data(), histoName.Data(), 50, 0., 200., 50, 0., 1.);
 151       histo->SetXTitle("p (GeV/c)");
 152       histo->SetYTitle("Chisquare prob.");
 153       AddObjectToCollection(histo);
 154
 155       histoName = GetHistoName(kSigmaVsPt, itheta, isrc);
 156       histo = new TH2F(histoName.Data(), histoName.Data(), 100, 0., 100., fSigmaCuts.GetSize(), 0.5, 0.5+(Double_t)fSigmaCuts.GetSize());
 157       histo->SetXTitle("p_{t} (GeV/c)");
 158       histo->SetYTitle("#sigma_{p#timesDCA}");
 159       for ( Int_t ibin=0; ibin<fSigmaCuts.GetSize(); ++ibin ) {
 160         histo->GetYaxis()->SetBinLabel(ibin+1,Form("%g", fSigmaCuts[ibin]));
 161       }
 162       AddObjectToCollection(histo);
 163
 164       histoName = GetHistoName(kSigmaVsEta, itheta, isrc);
 165       histo = new TH2F(histoName.Data(), histoName.Data(), 25, -4.5, -2., fSigmaCuts.GetSize(), 0.5, 0.5+(Double_t)fSigmaCuts.GetSize());
 166       histo->SetXTitle("#eta");
 167       histo->SetYTitle("#sigma_{p#timesDCA}");
 168       for ( Int_t ibin=0; ibin<fSigmaCuts.GetSize(); ++ibin ) {
 169         histo->GetYaxis()->SetBinLabel(ibin+1,Form("%g", fSigmaCuts[ibin]));
 170       }
 171       AddObjectToCollection(histo);
 172     } // loop on track sources
 173   } // loop on theta abs
 174
 175   fMuonTrackCuts->Print();
 176
 177 }
 178
 179 //________________________________________________________________________
 180 TString AliAnalysisTaskMuonCuts::GetHistoName(Int_t histoTypeIndex, Int_t thetaAbsIndex, Int_t srcIndex)
 181 {
 182   /// Get local histogram name
 183   TString histoName = Form("%s%s%s", fHistoTypeKeys->At(histoTypeIndex)->GetName(), fThetaAbsKeys->At(thetaAbsIndex)->GetName(), fSrcKeys->At(srcIndex)->GetName());
 184
 185   return histoName;
 186 }
 187
 188 //________________________________________________________________________
 189 void AliAnalysisTaskMuonCuts::ProcessEvent(TString physSel, const TObjArray& selectTrigClasses, TString centrality)
 190 {
 191   //
 192   /// Fill histogram
 193   //
 194
 195   if ( GetVertexSPD()->GetNContributors() < fMinNvtxContirbutors ) return;
 196
 197   TString histoName = "";
 198   AliVParticle* track = 0x0;
 199   Int_t nTracks = GetNTracks();
 200   for (Int_t itrack = 0; itrack < nTracks; itrack++) {
 201     track = GetTrack(itrack);
 202     fMuonTrackCuts->SetNSigmaPdca(1.e10);
 203     if ( ! fMuonTrackCuts->IsSelected(track) ) continue;
 204
 205     Double_t rAbsEnd =  ( fAODEvent ) ? ((AliAODTrack*)track)->GetRAtAbsorberEnd(): ((AliESDMuonTrack*)track)->GetRAtAbsorberEnd();
 206     Double_t thetaAbsEndDeg = TMath::ATan( rAbsEnd / 505. ) * TMath::RadToDeg();
 207     Int_t thetaAbsBin = ( thetaAbsEndDeg < 3. ) ? kThetaAbs23 : kThetaAbs310;
 208
 209     Int_t trackSrc = GetParticleType(track);
 210
 211     TVector3 dcaAtVz = fMuonTrackCuts->GetCorrectedDCA(track);
 212     Double_t pTotMean = fMuonTrackCuts->GetAverageMomentum(track);
 213     Double_t dca = dcaAtVz.Mag();
 214     Double_t pDca = pTotMean * dca;
 215
 216     Double_t chi2 = pDca / fMuonTrackCuts->GetSigmaPdca(rAbsEnd) ;
 217     chi2 *= chi2;
 218     Double_t chiProb = TMath::Prob(chi2, 2);
 219
 220     Double_t pTot = track->P();
 221     Double_t pt = track->Pt();
 222     Double_t eta = track->Eta();
 223
 224     for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
 225       TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
 226       histoName = GetHistoName(kDCAxVsP, thetaAbsBin, trackSrc);
 227       ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, dcaAtVz.X());
 228
 229       histoName = GetHistoName(kDCAyVsP, thetaAbsBin, trackSrc);
 230       ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, dcaAtVz.Y());
 231
 232       histoName = GetHistoName(kPdcaVsP, thetaAbsBin, trackSrc);
 233       ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, pDca);
 234
 235       histoName = GetHistoName(kPDCAVsPCheck, thetaAbsBin, trackSrc);
 236       ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, pDca);
 237
 238       histoName = GetHistoName(kDCAVsPCheck, thetaAbsBin, trackSrc);
 239       ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, dca);
 240
 241       histoName = GetHistoName(kChiProbVsP, thetaAbsBin, trackSrc);
 242       ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, chiProb);
 243     } // loop on selected trigger classes
 244
 245     for ( Int_t isigma=0; isigma<fSigmaCuts.GetSize(); ++isigma) {
 246       fMuonTrackCuts->SetNSigmaPdca(fSigmaCuts[isigma]);
 247       if ( ! fMuonTrackCuts->IsSelected(track) ) continue;
 248       for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
 249         TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
 250         histoName = GetHistoName(kSigmaVsPt, thetaAbsBin, trackSrc);
 251         ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pt, isigma+1);
 252         histoName = GetHistoName(kSigmaVsEta, thetaAbsBin, trackSrc);
 253         ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(eta, isigma+1);
 254       } // loop on selected trigger classes
 255     } // loop on sigmas
 256   }
 257 }
 258
 259 //________________________________________________________________________
 260 void AliAnalysisTaskMuonCuts::SetSigmaCuts(Int_t nSigmaCuts, Double_t* sigmaCuts)
 261 {
 262   /// Set number of sigmas
 263   if ( ! sigmaCuts || nSigmaCuts < 0 ) {
 264 //     if ( defaultChiSquare ) {
 265 //       Double_t cuts[] = {0.3, 0.1, 0.05, 0.03, 0.01, 0.005, 0.003, 0.001, 0.0005, 0.0003, 0.0001, 0.};
 266 //       Int_t nCuts = sizeof(cuts)/sizeof(cuts[0]);
 267 //       fSigmaCuts.Set(nCuts, cuts);
 268 //     }
 269 //     else {
 270       Double_t cuts[] = {2., 3., 4., 5., 6., 7., 10., 15., 20., 25., 30., 1.e10};
 271       Int_t nCuts = sizeof(cuts)/sizeof(cuts[0]);
 272       fSigmaCuts.Set(nCuts, cuts);
 273       //    }
 274   }
 275   else {
 276     fSigmaCuts.Set(nSigmaCuts, sigmaCuts);
 277   }
 278 }
 279
 280 //________________________________________________________________________
 281 void AliAnalysisTaskMuonCuts::Terminate(Option_t *) {
 282   //
 283   /// Draw some histogram at the end.
 284   //
 285
 286   AliVAnalysisMuon::Terminate("");
 287
 288   if ( ! fMergeableCollection ) return;
 289
 290   TString physSel = fTerminateOptions->At(0)->GetName();
 291   TString trigClassName = fTerminateOptions->At(1)->GetName();
 292   TString centralityRange = fTerminateOptions->At(2)->GetName();
 293   TString furtherOpt = fTerminateOptions->At(3)->GetName();
 294   furtherOpt.ToUpper();
 295
 296   furtherOpt.ReplaceAll("  ", " ");
 297   furtherOpt.ReplaceAll(" =", "=");
 298   furtherOpt.ReplaceAll("= ", "=");
 299   TObjArray* optArray = furtherOpt.Tokenize(" ");
 300   Double_t refSigmaCut = 6.;
 301   for ( Int_t iopt=0; iopt<optArray->GetEntries(); ++iopt ) {
 302     TString currOpt = optArray->At(iopt)->GetName();
 303     if ( currOpt.Contains("REFSIGMA") ) {
 304       currOpt.Remove(0,currOpt.Index("=")+1);
 305       refSigmaCut = currOpt.Atof();
 306     }
 307   }
 308   delete optArray;
 309
 310   Int_t srcColors[kNtrackSources] = {kBlack, kRed, kGreen, kBlue, kViolet, 7, kOrange};
 311
 312   TCanvas* can = 0x0;
 313   Int_t xshift = 100;
 314   Int_t yshift = 100;
 315   Int_t igroup1 = -1;
 316   Int_t igroup2 = 0;
 317
 318   //////////////
 319   // Reco DCA //
 320   //////////////
 321   igroup1++;
 322   igroup2 = 0;
 323   TString histoName = "", currName = "", histoPattern = "", drawOpt = "";
 324   currName = Form("%s_recoDCA", GetName());
 325   can = new TCanvas(currName.Data(),"Reco DCA",igroup1*xshift,igroup2*yshift,600,600);
 326   can->Divide(2,2);
 327   igroup2++;
 328   Int_t recoDcaHisto[2] = {kDCAxVsP, kDCAyVsP};
 329   TString dcaName[2] = {"DCAx", "DCAy"};
 330   Double_t averageDca[4] = {0., 0.};
 331   printf("\nAverage reconstructed DCA:\n");
 332   TF1* fitFuncMeanDca = new TF1("fitFuncMeanDca","gausn",-20.,20.);
 333   fitFuncMeanDca->SetParNames("Norm", "Mean", "Sigma");
 334   for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
 335     for ( Int_t ihisto=0; ihisto<2; ++ihisto ) {
 336       TH2* histo = 0x0;
 337       histoPattern = "";
 338       histoPattern = Form("%s & %s", fHistoTypeKeys->At(recoDcaHisto[ihisto])->GetName(), fThetaAbsKeys->At(itheta)->GetName());
 339       histo = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
 340       if ( ! histo ) continue;
 341
 342       TH1* meanDcaVsP = histo->ProjectionX(Form("mean%sVsP_%s", dcaName[ihisto].Data(), fThetaAbsKeys->At(itheta)->GetName()));
 343       meanDcaVsP->Reset();
 344       meanDcaVsP->SetTitle(Form("Mean %s vs. p %s", dcaName[ihisto].Data(), fThetaAbsKeys->At(itheta)->GetName()));
 345       meanDcaVsP->SetYTitle(Form("<%s> (cm)", dcaName[ihisto].Data()));
 346       meanDcaVsP->SetStats(kFALSE);
 347
 348       Int_t nPbins = histo->GetXaxis()->GetNbins();
 349       //Int_t nPadX = (Int_t)TMath::Sqrt(nPbins);
 350       //Int_t nPadY = nPadX;
 351       //if ( nPadX * nPadY < nPbins ) nPadX++;
 352       TCanvas* meanDcaFitCan = 0x0;
 353       Int_t nPadX = 5;
 354       Int_t nPadY = 5;
 355       Int_t ipad = 0;
 356       Int_t ican = 0;
 357
 358       for ( Int_t ibin=2; ibin<=nPbins; ++ibin ) {
 359         currName = Form("hMean%s_%s_%s", dcaName[ihisto].Data(), physSel.Data(), trigClassName.Data());
 360         Int_t minBin = ( ibin == 0 ) ? 1 : ibin;
 361         Int_t maxBin = ( ibin == 0 ) ? nPbins : ibin;
 362         if ( ibin > 0 ) currName += Form("_pBin%i", ibin);
 363         TH1* projHisto = histo->ProjectionY(currName.Data(), minBin, maxBin, "e");
 364         projHisto->SetTitle(Form("%s %s %g < p < %g (GeV/c)", dcaName[ihisto].Data(), fThetaAbsKeys->At(itheta)->GetName(), meanDcaVsP->GetXaxis()->GetBinLowEdge(minBin), meanDcaVsP->GetXaxis()->GetBinUpEdge(maxBin)));
 365
 366         if ( projHisto->GetEntries() == 0 ) continue;
 367         if ( ipad % (nPadX*nPadY) == 0 ) {
 368           currName = histo->GetName();
 369           currName += Form("Fit_can_%i", ican++);
 370           meanDcaFitCan = new TCanvas(currName.Data(), currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
 371           meanDcaFitCan->Divide(nPadX,nPadY);
 372           ipad = 0;
 373         }
 374         meanDcaFitCan->cd(++ipad);
 375         gPad->SetLogy();
 376         if ( projHisto->Integral() > 50 ) {
 377           fitFuncMeanDca->SetParameter(0, projHisto->Integral());
 378           fitFuncMeanDca->SetParameter(1, projHisto->GetMean());
 379           fitFuncMeanDca->SetParameter(2, projHisto->GetRMS());
 380           Double_t fitDcaLim = ( ibin == 0 ) ?  40. :  40./((Double_t)ibin);
 381           fitDcaLim = TMath::Max(5., fitDcaLim);
 382           projHisto->Fit(fitFuncMeanDca, "RQ", "e", -fitDcaLim, fitDcaLim);
 383           Double_t chi2 = fitFuncMeanDca->GetChisquare();
 384           Double_t ndf = fitFuncMeanDca->GetNDF();
 385           if ( ndf <= 0.) continue;
 386           if ( chi2 / ndf > 100. ) continue;
 387           Double_t meanDca = fitFuncMeanDca->GetParameter(1);
 388           Double_t meanDcaErr = fitFuncMeanDca->GetParError(1);
 389           meanDcaVsP->SetBinContent(ibin, meanDca);
 390           meanDcaVsP->SetBinError(ibin, meanDcaErr);
 391           //if ( ibin == 0 ) printf("%s  meanDca = %g +- %g\n", fThetaAbsKeys->At(itheta)->GetName(),  meanDca, meanDcaErr);
 392         }
 393         else projHisto->Draw("e");
 394       } // loop on momentum bins
 395       can->cd(2*itheta+ihisto+1);
 396       //meanDcaVsP->SetLineColor(srcColors[isrc]);
 397       //meanDcaVsP->SetMarkerColor(srcColors[isrc]);
 398       //meanDcaVsP->SetMarkerStyle(20+isrc);
 399       meanDcaVsP->Fit("pol0","Q");
 400       TF1* trendFit = (TF1*)meanDcaVsP->GetListOfFunctions()->FindObject("pol0");
 401       if ( trendFit ) {
 402         averageDca[2*itheta+ihisto] = trendFit->GetParameter(0);
 403         printf("  %s: mean %s = %g +- %g  (chi2/%i = %g)\n",  fThetaAbsKeys->At(itheta)->GetName(), dcaName[ihisto].Data(), trendFit->GetParameter(0), trendFit->GetParError(0), trendFit->GetNDF(), trendFit->GetChisquare()/((Double_t)trendFit->GetNDF()));
 404       }
 405         //drawOpt = "esame";
 406       //leg->AddEntry(meanDcaVsP, fSrcKeys->At(isrc)->GetName(), "lp");
 407       //} // loop on src
 408     //can->cd(itheta+1);
 409     //leg->Draw("same");
 410
 411       //can->cd(ipad++);
 412       //histo->Draw();
 413       //meanDca[ihisto] = histo->GetMean();
 414     } // loop on histo type
 415   } // loop on theta abs
 416
 417   for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
 418     printf("muonCuts->SetMeanDCA(%g, %g, 0.); // %s\n", averageDca[2*itheta], averageDca[2*itheta+1], fThetaAbsKeys->At(itheta)->GetName());
 419   }
 420
 421   //////////////
 422   // Fit pDCA //
 423   //////////////
 424   Double_t nSigmaCut = fMuonTrackCuts->GetNSigmaPdca(); //6.;
 425   Double_t sigmaMeasCut[2] = { fMuonTrackCuts->GetSigmaPdca(AliMuonTrackCuts::kThetaAbs23), fMuonTrackCuts->GetSigmaPdca(AliMuonTrackCuts::kThetaAbs23)}; //{99., 54.}; //{120., 63.};
 426   Double_t relPResolution = fMuonTrackCuts->GetRelPResolution(); //4.5e-4; //6.e-3;//8.e-4;
 427   Double_t angleResolution = 535.*fMuonTrackCuts->GetSlopeResolution(); //6.e-4;
 428   Double_t pMinCut = 0.1;
 429   Double_t pMaxCut =  800.;
 430   const Int_t kNCutFuncs = 2;
 431   Int_t nShowFuncs = 1;
 432   TString cutFormula = "[1]*TMath::Sqrt( ( [0] / ( 1. - [1]*[2]*x / ( 1.+[1]*[2]*x ) ) ) * ( [0] / ( 1. - [1]*[2]*x / ( 1.+[1]*[2]*x ) ) ) + [3]*[3]*x*x)";
 433   Double_t cutParam[kNCutFuncs][4] = {{sigmaMeasCut[0], nSigmaCut, relPResolution, angleResolution}, {sigmaMeasCut[0], nSigmaCut, 0., 0.32}};
 434   Int_t cutColor[kNCutFuncs] = {kBlack, kRed};
 435   igroup1++;
 436   igroup2 = 0;
 437   can = new TCanvas("pdcaSigmaFit","Sigma vs P fit",igroup1*xshift,igroup2*yshift,600,600);
 438   can->Divide(2,1);
 439   igroup2++;
 440   TF1* fitFunc = new TF1("fitFunc", "x*gausn", 0., 400.);
 441   fitFunc->SetParNames("Norm", "Mean", "Sigma");
 442   gStyle->SetOptFit(1111);
 443   Double_t xMinFit[2] = {0., 0.};
 444   Double_t xMaxFit[2] = {320., 150.}; // {360., 180.};
 445   printf("\nSigma p x DCA:\n");
 446   Double_t averageSigmaPdca[kNtrackSources*kNthetaAbs] = {0.};
 447   for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
 448     TLegend* leg = new TLegend(0.7, 0.7, 0.9, 0.9);
 449     leg->SetBorderSize(1);
 450     for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
 451       histoPattern = GetHistoName(kPdcaVsP, itheta, isrc);
 452       TH2* histo = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
 453       if ( ! histo ) continue;
 454       if ( histo->Integral() < 200 ) {
 455         delete histo;
 456         continue;
 457       }
 458
 459       TH1* sigmaVsP = histo->ProjectionX(Form("sigma%s_%s_%s", fHistoTypeKeys->At(kPdcaVsP)->GetName(), fThetaAbsKeys->At(itheta)->GetName(), fSrcKeys->At(isrc)->GetName()));
 460       sigmaVsP->Reset();
 461       sigmaVsP->SetTitle(Form("#sigma_{p#timesDCA} vs. p %s %s", fThetaAbsKeys->At(itheta)->GetName(), fSrcKeys->At(isrc)->GetName()));
 462       sigmaVsP->SetYTitle("#sigma_{p#timesDCA} (cm #times GeV/c)");
 463       sigmaVsP->SetStats(kFALSE);
 464
 465       Int_t nPbins = histo->GetXaxis()->GetNbins();
 466       //Int_t nPadX = (Int_t)TMath::Sqrt(nPbins);
 467       //Int_t nPadY = nPadX;
 468       //if ( nPadX * nPadY < nPbins ) nPadX++;
 469       TCanvas* pdcaFitCan = 0x0;
 470       Int_t nPadX = 5;
 471       Int_t nPadY = 5;
 472       Int_t ipad = 0;
 473       Int_t ican = 0;
 474
 475       for ( Int_t ibin=2; ibin<=nPbins; ++ibin ) {
 476         currName = Form("hPDCA_%s_%s_%s", fSrcKeys->At(isrc)->GetName(), physSel.Data(), trigClassName.Data());
 477         Int_t minBin = ( ibin == 0 ) ? 1 : ibin;
 478         Int_t maxBin = ( ibin == 0 ) ? nPbins : ibin;
 479         if ( ibin > 0 ) currName += Form("_pBin%i", ibin);
 480         TH1* projHisto = histo->ProjectionY(currName.Data(), minBin, maxBin, "e");
 481         if ( projHisto->GetEntries() == 0 ) continue;
 482         projHisto->SetTitle(Form("P DCA %s %s %g < p < %g (GeV/c)", fThetaAbsKeys->At(itheta)->GetName(), fSrcKeys->At(isrc)->GetName(), sigmaVsP->GetXaxis()->GetBinLowEdge(minBin), sigmaVsP->GetXaxis()->GetBinUpEdge(maxBin)));
 483         if ( ipad % (nPadX*nPadY) == 0 ) {
 484           currName = histo->GetName();
 485           currName += Form("Fit_can_%i", ican++);
 486           pdcaFitCan = new TCanvas(currName.Data(), currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
 487           pdcaFitCan->Divide(nPadX,nPadY);
 488           ipad = 0;
 489         }
 490         pdcaFitCan->cd(++ipad);
 491         if ( projHisto->Integral() > 0.) projHisto->Scale(1./projHisto->Integral());
 492         gPad->SetLogy();
 493         if ( projHisto->GetEntries() > 50 ) {
 494           fitFunc->SetParameter(0, projHisto->Integral());
 495           fitFunc->SetParameter(1, projHisto->GetMean());
 496           fitFunc->SetParameter(2, projHisto->GetRMS());
 497           projHisto->Fit(fitFunc, "RQ", "e", xMinFit[itheta], xMaxFit[itheta]);
 498           Double_t chi2 = fitFunc->GetChisquare();
 499           Double_t ndf = fitFunc->GetNDF();
 500           if ( ndf <= 0.) continue;
 501           if ( chi2 / ndf > 100. ) continue;
 502           Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
 503           Double_t sigmaErr = fitFunc->GetParError(2);
 504           sigmaVsP->SetBinContent(ibin, sigma);
 505           sigmaVsP->SetBinError(ibin, sigmaErr);
 506           //if ( ibin == 0 ) printf("%s %s  sigma = %g +- %g\n", fThetaAbsKeys->At(itheta)->GetName(), fSrcKeys->At(isrc)->GetName(), sigma, sigmaErr);
 507         }
 508         else projHisto->Draw("e");
 509       } // loop on momentum bins
 510       can->cd(itheta+1);
 511       sigmaVsP->SetLineColor(srcColors[isrc]);
 512       sigmaVsP->SetMarkerColor(srcColors[isrc]);
 513       sigmaVsP->SetMarkerStyle(20+isrc);
 514       drawOpt = "e";
 515       if ( gPad->GetListOfPrimitives()->GetEntries() != 0 ) drawOpt += "same";
 516       sigmaVsP->Draw(drawOpt.Data());
 517       sigmaVsP->Fit("pol0","Q",drawOpt.Data());
 518       TF1* trendFit = (TF1*)sigmaVsP->GetListOfFunctions()->FindObject("pol0");
 519       if ( trendFit ) {
 520         averageSigmaPdca[kNtrackSources*itheta+isrc] = trendFit->GetParameter(0);
 521         printf("  %s %s: Sigma = %g +- %g  (chi2/%i = %g)\n",  fThetaAbsKeys->At(itheta)->GetName(), fSrcKeys->At(isrc)->GetName(), trendFit->GetParameter(0), trendFit->GetParError(0), trendFit->GetNDF(), trendFit->GetChisquare()/((Double_t)trendFit->GetNDF()));
 522       }
 523       leg->AddEntry(sigmaVsP, fSrcKeys->At(isrc)->GetName(), "lp");
 524
 525       // Plot 2D function for check!
 526       histoPattern = GetHistoName(kPDCAVsPCheck, itheta, isrc);
 527       TH2* histoCheck = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
 528       if ( ! histoCheck ) histoCheck = histo; // needed for old data
 529       currName = histoCheck->GetName();
 530       currName.Append("_plotCut");
 531       TCanvas* pdcaVsPcan = new TCanvas(currName.Data(), currName.Data(), igroup1*xshift,(igroup2+1)*yshift,600,600);
 532       pdcaVsPcan->SetLogz();
 533       pdcaVsPcan->SetRightMargin(0.12);
 534       histoCheck->Draw("COLZ");
 535
 536       for ( Int_t icut=0; icut<nShowFuncs; ++icut ) {
 537         currName = Form("%s_cutFunc%i", histo->GetName(), icut);
 538         TF1* cutFunction = new TF1(currName.Data(),cutFormula.Data(), pMinCut, pMaxCut);
 539         cutParam[icut][0] = sigmaMeasCut[itheta];
 540         cutParam[icut][1] = nSigmaCut;
 541         cutFunction->SetParameters(cutParam[icut]);
 542         cutFunction->SetLineWidth(2);
 543         cutFunction->SetLineColor(cutColor[icut]);
 544         cutFunction->Draw("same");
 545       } // loop on cut func
 546     } // loop on src
 547     can->cd(itheta+1);
 548     leg->Draw("same");
 549
 550     for ( Int_t icut=0; icut<nShowFuncs; ++icut ) {
 551       currName = Form("sigmaCut_%s_%i", fThetaAbsKeys->At(itheta)->GetName(), icut);
 552       TF1* cutFunction = new TF1(currName.Data(), cutFormula.Data(), pMinCut, pMaxCut);
 553       cutParam[icut][0] = sigmaMeasCut[itheta];
 554       cutParam[icut][1] = 1.;
 555       cutFunction->SetParameters(cutParam[icut]);
 556       cutFunction->SetLineColor(cutColor[icut]);
 557       cutFunction->SetLineWidth(2);
 558       cutFunction->Draw("same");
 559     }
 560   } // loop on theta abs
 561
 562   for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
 563     printf("muonCuts->SetSigmaPdca(%g, %g); // %s\n", averageSigmaPdca[isrc], averageSigmaPdca[kNtrackSources+isrc], fSrcKeys->At(isrc)->GetName());
 564   }
 565   printf("\n");
 566
 567   igroup2++;
 568   for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
 569     for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
 570       histoPattern = GetHistoName(kDCAVsPCheck, itheta, isrc);
 571       TH2* histoCheck = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
 572       if ( ! histoCheck ) continue;
 573       currName = histoCheck->GetName();
 574       currName.Append("_plotCut");
 575       TCanvas* pdcaVsPcan = new TCanvas(currName.Data(), currName.Data(), igroup1*xshift,(igroup2+1)*yshift,600,600);
 576       pdcaVsPcan->SetRightMargin(0.12);
 577       pdcaVsPcan->SetLogz();
 578       histoCheck->Draw("COLZ");
 579
 580       for ( Int_t icut=0; icut<nShowFuncs; ++icut ) {
 581         currName = histoCheck->GetName();
 582         currName.Append(Form("_cutFunc%i",icut));
 583         TString currFormula = cutFormula;
 584         currFormula.Append("/x");
 585         TF1* cutFunction = new TF1(currName.Data(),currFormula.Data(), pMinCut, pMaxCut);
 586         cutParam[icut][0] = sigmaMeasCut[itheta];
 587         cutParam[icut][1] = nSigmaCut;
 588         cutFunction->SetParameters(cutParam[icut]);
 589         cutFunction->SetLineWidth(2);
 590         cutFunction->SetLineColor(cutColor[icut]);
 591         cutFunction->Draw("same");
 592       } // loop on cut functions
 593     } // loop on src
 594   } //loop on theta
 595
 596
 597   ///////////////////////////
 598   // Check Chi square prob //
 599   ///////////////////////////
 600   igroup1++;
 601   igroup2 = 0;
 602   for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
 603     for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
 604       histoPattern = GetHistoName(kChiProbVsP, itheta, isrc);
 605       TH2* histo = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
 606       if ( ! histo ) continue;
 607
 608       Int_t nPbins = histo->GetXaxis()->GetNbins();
 609       //Int_t nPadX = (Int_t)TMath::Sqrt(nPbins);
 610       //Int_t nPadY = nPadX;
 611       //if ( nPadX * nPadY < nPbins ) nPadX++;
 612       TCanvas* chi2probCan = 0x0;
 613       Int_t nPadX = 5;
 614       Int_t nPadY = 5;
 615       Int_t ipad = 0;
 616       Int_t ican = 0;
 617
 618       for ( Int_t ibin=0; ibin<=nPbins; ++ibin ) {
 619         currName = Form("hChiProb_%s_%s_%s", fSrcKeys->At(isrc)->GetName(), physSel.Data(), trigClassName.Data());
 620         Int_t minBin = ( ibin == 0 ) ? 1 : ibin;
 621         Int_t maxBin = ( ibin == 0 ) ? nPbins : ibin;
 622         if ( ibin > 0 ) currName += Form("_pBin%i", ibin);
 623         TH1* projHisto = histo->ProjectionY(currName.Data(), minBin, maxBin);
 624         projHisto->SetTitle(Form("Chisquare prob %s %s %g < p < %g (GeV/c)", fThetaAbsKeys->At(itheta)->GetName(), fSrcKeys->At(isrc)->GetName(), histo->GetXaxis()->GetBinLowEdge(minBin), histo->GetXaxis()->GetBinUpEdge(maxBin)));
 625         if ( projHisto->GetEntries() == 0 ) continue;
 626         if ( ipad % (nPadX*nPadY) == 0 ) {
 627           currName = histo->GetName();
 628           currName += Form("Fit_can_%i", ican++);
 629           chi2probCan = new TCanvas(currName.Data(), currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
 630           chi2probCan->Divide(nPadX,nPadY);
 631           ipad = 0;
 632         }
 633         chi2probCan->cd(++ipad);
 634         gPad->SetLogy();
 635         projHisto->Draw("e");
 636       } // loop on momentum bins
 637     } // loop on src
 638   } // loop on theta abs
 639
 640
 641   //////////////////////
 642   // Check sigma cuts //
 643   //////////////////////
 644   printf("\nReference sigma cut %g\n", refSigmaCut);
 645
 646   Float_t fracOfHeight = 0.35;
 647   Float_t rightMargin = 0.03;
 648   Int_t cutColors[14] = {kBlack, kRed, kBlue, kGreen, kCyan, kMagenta, kOrange, kViolet, kSpring, kGray, kSpring, kAzure, kPink, kYellow};
 649   Int_t* orderCuts = 0x0;
 650   Int_t nSigmaCuts = 0;
 651   Int_t checkHistos[2] = {kSigmaVsPt, kSigmaVsEta};
 652   Bool_t useCustomSigma = furtherOpt.Contains("CUSTOMSIGMA");
 653
 654   igroup1++;
 655   igroup2 = 0;
 656
 657   for ( Int_t ihisto=0; ihisto<2; ++ihisto ) {
 658     for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
 659       can = 0x0;
 660       for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
 661         histoPattern = GetHistoName(checkHistos[ihisto], itheta, isrc);
 662         TH2* histo = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
 663         if ( ! histo ) continue;
 664         if ( histo->Integral() == 0. ) {
 665           delete histo;
 666           continue;
 667         }
 668         if ( ! orderCuts ) {
 669           // Re-order axis
 670           TAxis* axis = histo->GetYaxis();
 671           nSigmaCuts = ( useCustomSigma ) ? fSigmaCuts.GetSize() : axis->GetNbins();
 672           orderCuts = new Int_t[nSigmaCuts];
 673           Int_t countBin = 0;
 674           for ( Int_t isigma=0; isigma<axis->GetNbins(); ++isigma ) {
 675             TString currLabel = axis->GetBinLabel(isigma+1);
 676             Double_t currSigma = currLabel.Atof();
 677             if ( useCustomSigma ) {
 678               Bool_t foundMatch = kFALSE;
 679               for ( Int_t jsigma=0; jsigma<fSigmaCuts.GetSize(); ++jsigma ) {
 680                 if ( TMath::Abs(currSigma - fSigmaCuts[jsigma]) < 1e-4 ) {
 681                   foundMatch = kTRUE;
 682                   break;
 683                 }
 684               }
 685               if ( ! foundMatch ) continue;
 686             }
 687             Int_t currBin = ( TMath::Abs(currSigma - refSigmaCut) < 1e-4) ? 0 : ++countBin;
 688             orderCuts[currBin] = isigma+1;
 689           }
 690         }
 691         currName = histo->GetName();
 692         currName.Append("_can");
 693         can = new TCanvas(currName.Data(), currName.Data(), igroup1*xshift,igroup2*yshift, 600, 600);
 694         TLegend* leg = new TLegend(0.6,0.6,0.9,0.9);
 695         leg->SetBorderSize(1);
 696         can->Divide(1,2,0,0);
 697         can->cd(1);
 698         gPad->SetPad(0., fracOfHeight, 0.99, 0.99);
 699         //gPad->SetTopMargin(0.08);
 700         gPad->SetTopMargin(0.03);
 701         gPad->SetRightMargin(rightMargin);
 702         gPad->SetLogy();
 703
 704         can->cd(2);
 705         gPad->SetPad(0., 0., 0.99, fracOfHeight);
 706         gPad->SetRightMargin(rightMargin);
 707         gPad->SetBottomMargin(0.08/fracOfHeight);
 708
 709         TH1* refCutHisto = 0x0;
 710         TString legTitle = "";
 711         for ( Int_t isigma=0; isigma<nSigmaCuts; ++isigma ) {
 712           currName = histo->GetName();
 713           currName.Append(Form("_sigma%i", isigma));
 714           Int_t currBin = orderCuts[isigma];
 715           TH1* projectHisto = histo->ProjectionX(currName.Data(), currBin, currBin);
 716           projectHisto->SetLineColor(cutColors[isigma]);
 717           projectHisto->SetMarkerColor(cutColors[isigma]);
 718           projectHisto->SetMarkerStyle(20+isigma);
 719           projectHisto->SetStats(kFALSE);
 720           projectHisto->SetTitle("");
 721           can->cd(1);
 722           drawOpt = "e";
 723           if ( gPad->GetListOfPrimitives()->GetEntries() != 0 ) drawOpt += "same";
 724           projectHisto->Draw(drawOpt.Data());
 725           TString currLabel = histo->GetYaxis()->GetBinLabel(currBin);
 726           TString cutName = ( TMath::Abs(currLabel.Atof() ) > 1000. ) ? "Total" :  Form("%s #sigma cut", currLabel.Data());
 727           leg->AddEntry(projectHisto, cutName.Data(), "lp");
 728           if ( ! refCutHisto ) {
 729             refCutHisto = projectHisto;
 730             legTitle = Form("(pass n #sigma) / (pass %s #sigma)", currLabel.Data());
 731             continue;
 732           }
 733           currName.Append("_ratio");
 734           TH1* histoRatio = (TH1*)projectHisto->Clone(currName.Data());
 735           histoRatio->Sumw2();
 736           histoRatio->Divide(refCutHisto);
 737           histoRatio->SetTitle("");
 738           histoRatio->GetYaxis()->SetTitle(legTitle.Data());
 739           histoRatio->GetXaxis()->SetLabelSize(0.04/fracOfHeight);
 740           histoRatio->GetXaxis()->SetTitleSize(0.035/fracOfHeight);
 741           histoRatio->GetYaxis()->SetLabelSize(0.03/fracOfHeight);
 742           histoRatio->GetYaxis()->SetTitleSize(0.03/fracOfHeight);
 743           histoRatio->GetXaxis()->SetTitleOffset(1.);
 744           histoRatio->GetYaxis()->SetTitleOffset(0.6);
 745           histoRatio->GetYaxis()->SetRangeUser(0.5,1.5);
 746
 747           can->cd(2);
 748           drawOpt = "e";
 749           if ( gPad->GetListOfPrimitives()->GetEntries() != 0 ) drawOpt += "same";
 750           histoRatio->Draw(drawOpt.Data());
 751         }// loop on sigma cuts
 752         can->cd(1);
 753         leg->Draw("same");
 754       } // loop on theta abs
 755     } // loop on src
 756   } // loop on histo type
 757
 758   igroup1++;
 759   igroup2=0;
 760   Double_t ptMin[] = {0., 2., 4., 15.};
 761   Int_t nPtMins = sizeof(ptMin)/sizeof(ptMin[0]);
 762   for ( Int_t iptmin=0; iptmin<nPtMins; ++iptmin) {
 763     for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
 764       can = 0x0;
 765       for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
 766         TLegend* leg = 0x0;
 767         histoName = GetHistoName(kSigmaVsPt, itheta, isrc);
 768         for ( Int_t icent=1; icent<=fCentralityClasses->GetNbins(); ++icent ) {
 769           TH2* histo = (TH2*)GetSum(physSel, trigClassName, fCentralityClasses->GetBinLabel(icent), histoName);
 770           if ( ! histo ) continue;
 771           Int_t ptMinBin = histo->GetXaxis()->FindBin(ptMin[iptmin]);
 772           Int_t ptMaxBin = histo->GetXaxis()->GetNbins()+1;
 773           currName = histo->GetName();
 774           currName += Form("_%s_ptMin%i", fCentralityClasses->GetBinLabel(icent), TMath::Nint(ptMin[iptmin]));
 775           TH1* projectHisto = histo->ProjectionY(currName.Data(), ptMinBin, ptMaxBin);
 776           if ( projectHisto->GetMaximum() < 50. ) {
 777             delete projectHisto;
 778             continue;
 779           }
 780           if ( ! can ) {
 781             currName = Form("CutEffect_%s_ptMin%i", fSrcKeys->At(isrc)->GetName(), TMath::Nint(ptMin[iptmin]));
 782             can = new TCanvas(currName.Data(), currName.Data(), igroup1*xshift,igroup2*yshift, 600, 600);
 783             can->Divide(2,1);
 784           }
 785           if ( ! leg ) {
 786             leg = new TLegend(0.6,0.6,0.9,0.9);
 787             leg->SetBorderSize(1);
 788           }
 789           can->cd(itheta+1);
 790           projectHisto->SetTitle(Form("Cut effect %s %s %g < p_{t} < %g (GeV/c)", fSrcKeys->At(isrc)->GetName(), fThetaAbsKeys->At(itheta)->GetName(), histo->GetXaxis()->GetBinLowEdge(ptMinBin), histo->GetXaxis()->GetBinUpEdge(ptMaxBin)));
 791           projectHisto->SetLineColor(cutColors[icent-1]);
 792           projectHisto->SetMarkerColor(cutColors[icent-1]);
 793           projectHisto->SetMarkerStyle(19+icent);
 794           projectHisto->SetStats(0);
 795           drawOpt = "e";
 796           if ( gPad->GetListOfPrimitives()->GetEntries() != 0 ) drawOpt += "same";
 797           projectHisto->Draw(drawOpt.Data());
 798           leg->AddEntry(projectHisto, fCentralityClasses->GetBinLabel(icent), "lp");
 799           Double_t keptEvents = projectHisto->GetBinContent(orderCuts[0]);
 800           Double_t totalEvents = projectHisto->GetBinContent(orderCuts[nSigmaCuts-1]);
 801           Double_t accepted = ( totalEvents > 0. ) ? keptEvents / totalEvents : 1.;
 802           Double_t acceptedErr = ( totalEvents > 0. ) ? TMath::Sqrt(accepted*(1.-accepted)/totalEvents) : 1.;
 803           printf("%12s %11s %6s (pt>%g) rejected evts : %6.2f +- %6.3f %%\n", fSrcKeys->At(isrc)->GetName(), fThetaAbsKeys->At(itheta)->GetName(), fCentralityClasses->GetBinLabel(icent), ptMin[iptmin], (1.-accepted)*100., acceptedErr*100.);
 804           //printf("  rejected %g  total %g   (%s vs %s)\n",totalEvents-keptEvents,totalEvents,projectHisto->GetXaxis()->GetBinLabel(orderCuts[0]),projectHisto->GetXaxis()->GetBinLabel(orderCuts[nSigmaCuts-1]));
 805         } // loop on centrality
 806         if ( leg ) leg->Draw("same");
 807       } // loop on theta abs
 808     } // loop on sources
 809   } // loop on pt min
 810   delete [] orderCuts;
 811 }