PWGLF/STRANGENESS/LambdaK0/AliV0Module.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /***********************************************
  17
  18   Lambda Analysis Module
  19   ----------------------
  20
  21 This Class enables the analysis of output files
  22 created with AliAnalysisTaskExtractV0 and
  23 AliAnalysisTaskExtractPerformanceV0 grid tasks.
  24
  25 It constructs corrected Lambda, AntiLambda and
  26 K0Short spectra.
  27
  28 This version: 27th April 2012
  29
  30 To be compiled as a library
  31 (.L AliV0Module.cxx++ or something to that effect)
  32
  33 --- David Dobrigkeit Chinellato
  34     daviddc@ifi.unicamp.br
  35
  36 ***********************************************/
  37
  38 //--- For C++ ----
  39 #include <TApplication.h>
  40 #include <TMatrix.h>
  41 #include <TMatrixD.h>
  42 #include <TROOT.h>
  43 #include <TMath.h>
  44 #include <TF1.h>
  45 #include <TH1F.h>
  46 #include <TH2F.h>
  47 #include <TH3F.h>
  48 #include <TList.h>
  49 #include <TLine.h>
  50 #include <TFile.h>
  51 #include <TTree.h>
  52 #include <TCanvas.h>
  53 #include <TString.h>
  54 #include <THashList.h>
  55 #include <TDirectoryFile.h>
  56 #include <cstdlib>
  57 using namespace std;
  58
  59 //--- For ROOT ---
  60 #include "AliV0Module.h"
  61
  62 AliV0Module::AliV0Module()
  63 {
  64   //Dummy Constructor that sets some default values
  65   //so that nothing too horrible will happen (hopefully)
  66   fWhichParticle           = "Lambda"; //Default
  67   fCINT1BoverINELratio     = 1;
  68   fRapidityBoundary        = 0.5;
  69   fRealDataFile            = "";
  70   fMCDataFile              = "";
  71   fFeedDownDataFile        = "";
  72   fOutputDataFile          = "";
  73
  74   //Selections
  75   fCutV0Radius                     = -1;
  76   fCutDCANegToPV                   = -1;
  77   fCutDCAPosToPV                   = -1;
  78   fCutDCAV0Daughters               = 1000;
  79   fCutV0CosPA                      = -2;
  80   fCutProperLifetime               = 1e+6;
  81   fCutTPCPIDNSigmas                = 1e+6;
  82   fCutNSigmasForSignalExtraction   = 5;
  83   fCutLeastNumberOfCrossedRows     = 70;
  84   fCutDaughterEta                  = 0.8;
  85   fCutCompetingV0Rejection         = -1;
  86
  87   //Set Feeddown Treatment
  88   fFDSwitch = "UseMCRatio";
  89
  90   //Default: Use bin counting
  91   fFitBackgroundSwitch = kFALSE;
  92
  93   //Default: Min-Bias
  94   fPerformMultiplicityStudy = kFALSE;
  95   fLoMultBound       = -1;
  96   fHiMultBound       = 10000;
  97
  98   //Default: no cut in Armenteros
  99   fSpecialArmenterosCutK0s = kFALSE;
 100
 101   //Pt Bins: undefined
 102   fptbinnumb = -1;
 103 }
 104
 105 AliV0Module::AliV0Module(TString fParticleType)
 106 {
 107   // Allows definition of Particle Type in analysis.
 108   // Possible Options are "Lambda", "AntiLambda" and "K0Short".
 109   // If some other string is given, this constructor will
 110   // default to "Lambda".
 111   fWhichParticle = fParticleType;
 112   if(fWhichParticle!="Lambda"&&fWhichParticle!="AntiLambda"&&fWhichParticle!="K0Short"){
 113     cout<<"Particle Type "<<fParticleType<<" unknown. Set to lambda."<<endl;
 114     fWhichParticle = "Lambda";
 115   }
 116   fCINT1BoverINELratio     = 1;
 117   fRapidityBoundary        = 0.5;
 118   fRealDataFile            = "";
 119   fMCDataFile              = "";
 120   fFeedDownDataFile        = "";
 121   fOutputDataFile          = "";
 122
 123   //Selections
 124   fCutV0Radius                     = -1;
 125   fCutDCANegToPV                   = -1;
 126   fCutDCAPosToPV                   = -1;
 127   fCutDCAV0Daughters               = 1000;
 128   fCutV0CosPA                      = -2;
 129   fCutProperLifetime               = 1e+6;
 130   fCutTPCPIDNSigmas                = 1e+6;
 131   fCutNSigmasForSignalExtraction   = 5;
 132   fCutLeastNumberOfCrossedRows     = 70;
 133   fCutDaughterEta                  = 0.8;
 134   fCutCompetingV0Rejection         = -1;
 135
 136   //Set Feeddown Treatment
 137   fFDSwitch = "UseMCRatio";
 138
 139   //Default: Use bin counting
 140   fFitBackgroundSwitch = kFALSE;
 141
 142   //Default: Min-Bias
 143   fPerformMultiplicityStudy = kFALSE;
 144   fLoMultBound       = -1;
 145   fHiMultBound       = 10000;
 146
 147   //Default: no cut in Armenteros
 148   fSpecialArmenterosCutK0s = kFALSE;
 149
 150   //Pt Bins: undefined
 151   fptbinnumb = -1;
 152 }
 153
 154 /***********************************************
 155  --- Setters For Configuration ---
 156 ***********************************************/
 157
 158 // Filename Setters
 159 void AliV0Module::SetRealDataFile    ( TString RealDataFilename     ){
 160   //Set root file containing real data candidates.
 161   fRealDataFile = RealDataFilename; }
 162 void AliV0Module::SetMCDataFile      ( TString MCDataFilename       ){
 163   //Set root file containing Monte Carlo data (for efficiency computation).
 164   fMCDataFile   = MCDataFilename;   }
 165 void AliV0Module::SetFeedDownDataFile( TString FeedDownDataFilename ){
 166   //Set root file containing Monte Carlo data (for feeddown computation).
 167   fFeedDownDataFile = FeedDownDataFilename; }
 168 void AliV0Module::SetOutputFile      ( TString OutputFilename       ){
 169   //Set root filename for the analysis output.
 170   fOutputDataFile = OutputFilename;
 171   cout<<"[AliV0Module] Set output file to \'"<<OutputFilename<<"\'."<<endl;
 172 }
 173
 174 // Bin Limit Setter
 175 void AliV0Module::SetPtBinLimits(Long_t got_fptbinnumb,const Double_t *got_fptbinlimits){
 176   //Function to set pt binning. First argument is the number of pt bins, second is
 177   //an array with bin limits.
 178   fptbinnumb = got_fptbinnumb;
 179         for(int ix = 0;ix<fptbinnumb+1;ix++){
 180                 fptbinlimits[ix] = got_fptbinlimits[ix];
 181         }
 182         for(int ix = 0;ix<fptbinnumb;ix++){
 183                 fptX[ix] = (fptbinlimits[ix+1] + fptbinlimits[ix])/2.;
 184         }
 185   cout<<"[AliV0Module] Received "<<fptbinnumb<<" pt bins, set accordingly."<<endl;
 186 }
 187
 188 // Rapidity Window Setter
 189 void AliV0Module::SetRapidityWindow(Double_t got_fRapidityBoundary){
 190   //Set Rapidity Boundary used in analysis.
 191   //Value provided will be used as upper limit for the modulus of y (|y|< given value)
 192   fRapidityBoundary = got_fRapidityBoundary;
 193   cout<<"[AliV0Module] Received "<<got_fRapidityBoundary<<" as rapidity limits, set accordingly."<<endl;
 194 }
 195
 196 // CINT1B/INEL Setter for normalization to yields
 197 void AliV0Module::SetCINT1BoverINEL(Double_t got_fCINT1BoverINEL){
 198   //Set CINT1B/INEL ratio (determined by collaboration).
 199   fCINT1BoverINELratio = got_fCINT1BoverINEL;
 200   cout<<"[AliV0Module] Received CINT1B/INEL = "<<got_fCINT1BoverINEL<<" to normalize with, set accordingly."<<endl;
 201 }
 202
 203 // Topological Selection Setters
 204 void AliV0Module::SetCutV0Radius(Double_t cut){
 205   //Set minimum decay radius for the V0 in centimeters.
 206   //Note: this is (R_{2D}, cylindrical, centered around ALICE detector)
 207   fCutV0Radius = cut;
 208   cout<<"[AliV0Module] Received V0 Radius (min value) = "<<cut<<endl;
 209 }
 210 void AliV0Module::SetCutDCANegToPV(Double_t cut){
 211   //Set minimum distance of closest approach between V0 negative daughter
 212   //track and primary vertex (in centimeters).
 213   fCutDCANegToPV = cut;
 214   cout<<"[AliV0Module] Received DCA Negative track to PV (min value) = "<<cut<<endl;
 215 }
 216 void AliV0Module::SetCutDCAPosToPV(Double_t cut){
 217   //Set minimum distance of closest approach between V0 positive daughter
 218   //track and primary vertex (in centimeters).
 219   fCutDCAPosToPV = cut;
 220   cout<<"[AliV0Module] Received DCA Positive track to PV (min value) = "<<cut<<endl;
 221 }
 222 void AliV0Module::SetCutDCAV0Daughters(Double_t cut){
 223   //Set minimum distance of closest approach between V0 daughter
 224   //tracks, in sigmas. This is not in centimeters because if the
 225   //tracks have been determined with ITS refit the resolution is
 226   //greatly improved; thus, the cut may be tighter in that case.
 227   //Using sigmas will therefore take the tracking method in
 228   //consideration.
 229   fCutDCAV0Daughters = cut;
 230   cout<<"[AliV0Module] Received DCA V0 Daughters (max value) = "<<cut<<endl;
 231 }
 232 void AliV0Module::SetCutV0CosPA(Double_t cut){
 233   //Set minimum value for the cosine of pointing angle of the V0.
 234   fCutV0CosPA = cut;
 235   cout<<"[AliV0Module] Received V0 Cosine of Pointing Angle (min value) = "<<cut<<endl;
 236 }
 237
 238 // Other Selection Setters
 239 void AliV0Module::SetCutProperLifetime(Double_t cut){
 240   //Set maximum value for m*L/p variable for the V0.
 241   //This is the "proper lifetime" selection and is usually called a
 242   //"c*tau cut". Should be set to a value larger than the c*tau for
 243   //the V0 considered.
 244   fCutProperLifetime = cut;
 245   cout<<"[AliV0Module] Received proper lifetime cut (max value) = "<<cut<<endl;
 246 }
 247 void AliV0Module::SetCutTPCPIDNSigmas(Double_t cut){
 248   //Set maximum deviation from the expected energy loss in
 249   //the TPC, in multiples of sigmas as computed from the AliPIDResponse
 250   //object. Selection is only used in real data and should thus
 251   //be very loose to ensure negligible signal loss.
 252   fCutTPCPIDNSigmas = cut;
 253   cout<<"[AliV0Module] Received TPC N-sigmas selection (max dist from BB curve) = "<<cut<<endl;
 254 }
 255 void AliV0Module::SetCutSigmaForSignalExtraction(Double_t cut){
 256   //Set number of sigmas for the signal extraction method. The value
 257   //provided is the number of sigmas from the peak position used for
 258   //the peak sampling, i.e. the peak will be from [<m>-cut*sigma,<m>+cut*sigma]
 259   //while the background sampling regions will be from
 260   //[<m>-2*cut*sigma,<m>+2*cut*sigma.
 261   fCutNSigmasForSignalExtraction = cut;
 262   cout<<"[AliV0Module] Received N-sigmas for sig. ext.: peak is (-"<<cut<<",+"<<cut<<") in sigmas"<<endl;
 263 }
 264 void AliV0Module::SetCutLeastNumberOfCrossedRows(Double_t cut){
 265   //Set smallest allowed number of TPC clusters for the V0 daughter tracks.
 266   fCutLeastNumberOfCrossedRows = cut;
 267   cout<<"[AliV0Module] Received Least Nbr of crossed rows (min value) = "<<cut<<endl;
 268 }
 269 void AliV0Module::SetCutLeastNumberOfCrossedRowsOverFindable(Double_t cut){
 270   //Set smallest allowed number of TPC clusters for the V0 daughter tracks.
 271   fCutLeastNumberOfCrossedRowsOverFindable = cut;
 272   cout<<"[AliV0Module] Received Least Nbr of crossed rows over findable clusters (min value) = "<<cut<<endl;
 273 }
 274
 275 void AliV0Module::SetCutDaughterEta(Double_t cut){
 276   //Set maximum eta value allowed for the V0 daughter tracks (|eta|<cut).
 277   fCutDaughterEta = cut;
 278   cout<<"[AliV0Module] Received Daughter |eta| cut (max value) = "<<cut<<endl;
 279 }
 280
 281 void AliV0Module::SetCutCompetingV0Rejection(Double_t cut){
 282   //Set rejection window around invariant mass of competing V0 species.
 283   //If negative, no rejection will occur. If analysis is for Lambdas, K0s
 284   //will be rejected and vice-versa. Typical values revolve around
 285   //0.003 - 0.010 GeV/c^2.
 286   fCutCompetingV0Rejection = cut;
 287   cout<<"[AliV0Module] Received Competing V0 Rejection window of +/- (in GeV/c^2) = "<<cut<<endl;
 288 }
 289
 290 void AliV0Module::SetFeeddownTreatment ( TString fFDMethod ){
 291   //Set method used to compute charged and neutral Xi baryon feeddown to Lambda.
 292   //Methods allowed are:
 293   //"NoFD" - No Feeddown subtraction performed.
 294   //"DoubleChargedXi" - Feeddown is computed for the charged Xi, and is then
 295   //subtracted twice from real data. Assumes charged and neutral Xi enter
 296   //the analysis in the same way and are produced at a ratio 1:1.
 297   //"UseMCRatio" - The Feeddown matrix F_{ij} is filled not only with Lambdas
 298   //coming from charged Xis but also from neutral ones. The scaling is performed
 299   //according to the measured charged Xis, and since the matrix is more populated
 300   //a larger subtraction will occurr. This method assumes that MC correctly
 301   //reproduces the ratio between charged and neutral Xi baryons. Warning: this
 302   //should not be used with charged Xi triggered Monte Carlo data.
 303   fFDSwitch = fFDMethod;
 304   cout<<"[AliV0Module] Received Feeddown treatment method: "<<fFDMethod<<endl;
 305   if( fFDMethod == "NoFD" )
 306     cout<<"[AliV0Module] ---> No Feeddown correction will be performed."<<endl;
 307   if( fFDMethod == "DoubleChargedXi" )
 308     cout<<"[AliV0Module] ---> Feeddown performed by doubling charged Xi correction."<<endl;
 309   if( fFDMethod == "UseMCRatio" )
 310     cout<<"[AliV0Module] ---> Feeddown performed by using MC neutral/charged Xi."<<endl;
 311 }
 312
 313 void AliV0Module::SetFitBackground ( Bool_t fitBgSwitch ){
 314   //Turns on background fitting for signal extraction instead of pure
 315   //bin counting. Useful, among other things, for systematics.
 316   fFitBackgroundSwitch = fitBgSwitch;
 317 }
 318
 319 void AliV0Module::SetPerformMultiplicityStudy ( Bool_t lPerformMultStudy ){
 320   //Turns on selection according to multiplicity of the event.
 321   //Boundaries are set in charged track multiplicity (pp) or in
 322   //centrality percentiles (PbPb). This is a requirement for studying
 323   //PbPb.
 324   fPerformMultiplicityStudy = lPerformMultStudy;
 325 }
 326
 327 void AliV0Module::SetLowMultValue ( Int_t lLoMultBound       ){
 328   //Lower boundary (inclusive) in integer number for mult selection.
 329   //Note: If in PbPb and you want, say, 10-20% centrality, set this
 330   //to 10.
 331   fLoMultBound = lLoMultBound;
 332 }
 333
 334 void AliV0Module::SetHighMultValue ( Int_t lHiMultBound       ){
 335   //Lower boundary (inclusive) in integer number for mult selection.
 336   //Note: If in PbPb and you want, say, 10-20% centrality, set this
 337   //to 10.
 338   fHiMultBound = lHiMultBound;
 339 }
 340
 341 void AliV0Module::SetSpecialArmenterosCutK0s ( Bool_t lSpecialArmenterosCutK0s ){
 342   //Special armenteros cut: |alpha|<5*pt_{arm}
 343   fSpecialArmenterosCutK0s = lSpecialArmenterosCutK0s;
 344 }
 345
 346 void AliV0Module::SetDefaultCuts(){
 347    //Sets Default cuts for analysis. (adjusted for adequate pp analysis)
 348   cout<<"[AliV0Module] Setting default cuts for particle species: "<<fWhichParticle<<endl;
 349   //Set Cuts - topological
 350   SetCutV0Radius       (0.500);
 351   SetCutDCANegToPV     (0.060);
 352   SetCutDCAPosToPV     (0.060);
 353   SetCutDCAV0Daughters (1.000);
 354   if ( fWhichParticle == "K0Short")
 355     SetCutV0CosPA      (0.970);
 356   if ( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" )
 357     SetCutV0CosPA      (0.995);
 358
 359
 360   //Set Cuts - other
 361   if ( fWhichParticle == "K0Short")
 362     SetCutProperLifetime                          ( 20);
 363   if ( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" )
 364     SetCutProperLifetime                          ( 30);
 365
 366   SetCutTPCPIDNSigmas                           (  5);
 367   SetCutSigmaForSignalExtraction                (  6);
 368   SetCutLeastNumberOfCrossedRows                ( 70);
 369   SetCutLeastNumberOfCrossedRowsOverFindable    (0.8);
 370   SetCutDaughterEta                             (0.8);
 371   if ( fWhichParticle == "K0Short")
 372     SetCutCompetingV0Rejection      (0.005);
 373   if ( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" )
 374     SetCutCompetingV0Rejection      (0.010);
 375 }
 376
 377 TString AliV0Module::IntToString(int input){
 378   //Integer to TString Converter
 379         char dummyChar[50];
 380         sprintf(dummyChar, "%d", input);
 381         TString outputstring = dummyChar;
 382         return outputstring;
 383 }
 384
 385 TString AliV0Module::DoubleToString(double input){
 386   //Double to TString Converter
 387         char dummyChar[50];
 388         sprintf(dummyChar, "%.3f", input);
 389         TString outputstring = dummyChar;
 390         return outputstring;
 391 }
 392
 393 Double_t AliV0Module::ErrorInRatio ( Double_t A, Double_t Aerr, Double_t B, Double_t Berr ){
 394   //Error in a Ratio
 395         if(B!=0){
 396                 Double_t errorfromtop = Aerr*Aerr / (B*B) ;
 397                 Double_t errorfrombottom = ((A*A)/(B*B*B*B)) * Berr * Berr;
 398                 return TMath::Sqrt( errorfromtop + errorfrombottom );
 399         }
 400         return 1;
 401 }
 402
 403 Double_t AliV0Module::MyGeant3FlukaCorrectionForProtons(const Double_t *x, const Double_t *par){
 404         //Parametrization Used for Geant3/Fluka Correction for protons
 405         //Credit: Antonin Maire
 406         return 1 - par[0]*TMath::Exp(par[1]*x[0]) + par[2];
 407 }
 408
 409
 410 Double_t AliV0Module::MyGeant3FlukaCorrectionForAntiProtons(const Double_t *x, const Double_t *par){
 411         // Parametrization Used for Geant3/Fluka Correction for antiprotons
 412         // Credit: Antonin Maire
 413         //
 414         // La fonction A*TMath::Erf( B.x ) ne marche pas à bas pt.
 415         // Différentes fonctions ont été testées.
 416         // La fonction suivante semble donner de bons résultats
 417         // On peut jouer sur la puissance n du terme "1/x^n" pour repousser le pt auquel la fonction prend la valeur 1.0
 418         // Ici, pour n = 0.2, la fonction prend la valeur 0.9990 en pt = 10 GeV/c
 419
 420         return 1 - par[0]*TMath::Exp(par[1]*x[0]) + par[2] + par[3]*1/TMath::Power(x[0], 0.2)*TMath::Log(x[0]);
 421 }
 422
 423 Double_t AliV0Module::MyLevyPtXi(const Double_t *pt, const Double_t *par)
 424 {
 425   //Levy Fit Function
 426   Double_t lMass  = 1.32171; //Xi Mass
 427   Double_t ldNdy  = par[0];
 428   Double_t lTemp = par[1];
 429   Double_t lPower = par[2];
 430
 431   Double_t lBigCoef = ((lPower-1)*(lPower-2)) / (lPower*lTemp*(lPower*lTemp+lMass*(lPower-2)));
 432   Double_t lInPower = 1 + (TMath::Sqrt(pt[0]*pt[0]+lMass*lMass)-lMass) / (lPower*lTemp);
 433
 434   return ldNdy * pt[0] * lBigCoef * TMath::Power(lInPower,(-1)*lPower);
 435   //return ldNdy * lBigCoef * TMath::Power(lInPower,(-1)*lPower);
 436 }
 437
 438 Double_t AliV0Module::MyBgPol1(const Double_t *x, const Double_t *par)
 439 {
 440   //Function for background fitting, rejects peak region
 441         if ( x[0] > par[2] && x[0] < par[3]) {
 442                 TF1::RejectPoint();
 443                 return 0;
 444         }
 445         return par[0] + par[1]*x[0];
 446 }
 447
 448 Double_t AliV0Module::MyBgPolToEval1(const Double_t *x, const Double_t *par)
 449 {
 450   //Just a plain linear function.
 451         return par[0] + par[1]*x[0];
 452 }
 453
 454 Double_t AliV0Module::RoundToThousandth( const Double_t lToRound ){
 455   //Round any number to a hundredth...
 456   //Well, within machine precision...
 457   return TMath::Nint( 1000 * lToRound ) * 0.001;
 458 }
 459
 460 void AliV0Module::DoAnalysis(){
 461   //----------------------
 462   // V0 Analysis Function
 463   //----------------------
 464   //
 465   //Consists of the following steps:
 466   //
 467   // (1) Loop over Real data candidates, acquire peak position and widths.
 468   // (2) Loop over Real data candidates, extract signal with variable extraction
 469   //     areas. Two loops: totally avoids binning data, allowing for sub-MeV/c^2
 470   //     granularity in signal extraction areas.
 471   // (3) Loop over MC data reconstructed candidates, find associated-to-MC primary
 472   //     candidates for efficiency numerator.
 473   // (4) (if Lambda or AntiLambda) Perform Geant3/Fluka correction.
 474   // (5) Get generated primary V0 histograms from MC file for efficiency denominator.
 475   // (6) (if Lambda or AntiLambda + FD correction enabled) Open MC file for feeddown
 476   //     subtraction. Loop over candidates, find Lambda associated with primary Xi
 477   //     to fill feeddown matrix. Scale the feeddown contribution according to real-life
 478   //     measured Xi- production under specified feeddown subtraction scheme (see
 479   //     function SetFeeddownTreatment). Perform subtraction of raw Lambda or AntiLambda
 480   //     count estimated to be coming from charged or neutral Xi.
 481   // (7) Perform detection efficiency correction and compute final corrected spectra.
 482   //
 483   //
 484   // Normalization:
 485   //  --- Number of Inelastic Events estimated to be:
 486   //      <Number of CINT1B triggers> / <CINT1B/INEL ratio>
 487   //
 488   //  --- Efficiency denominator: filled at pre-physics selection stage. All signal
 489   //      losses are therefore estimated from Monte Carlo.
 490   //
 491   //  --- Pileup: not included in MC. Fraction of events removed by IsPileupFromSPD
 492   //      is artificially removed from Number of Inelastic Events as well.
 493   //
 494   // Output: written to specified file (with SetOutputFile(...)).
 495   //   --- includes spectra and efficiencies in base directory.
 496   //   --- includes a number of different subdirectories storing plots such as
 497   //       such as invariant mass histos, signal extraction, pt resolution, G3/F
 498   //       related distributions, feeddown related distributions, and so on.
 499
 500   //Set Batch Mode: Ignore all canvases in this section
 501   gROOT->SetBatch (kTRUE);
 502
 503   cout<<"======================================"<<endl;
 504   cout<<" -------- Spectra Extraction -------- "<<endl;
 505   cout<<"======================================"<<endl;
 506   cout<<endl;
 507   if(fptbinnumb == -1){
 508     cout<<"[AliV0Module] It's unclear what kind of pt binning you want."<<endl;
 509     cout<<"[AliV0Module] Most likely you forgot to set it using SetPtBinLimits..."<<endl;
 510     cout<<"[AliV0Module] Analysis will NOT be done. Returning."<<endl;
 511     return;
 512   }
 513   cout<<"--------------- Configuration --------------------------"<<endl;
 514   cout<<" Analysed Particle.............: "<<fWhichParticle<<endl;
 515   cout<<" Rapidity Window...............: "<<fRapidityBoundary<<endl;
 516   cout<<" CINT1B/INEL Ratio used........: "<<fCINT1BoverINELratio<<endl;
 517   cout<<" V0 Decay Radius...............: "<<fCutV0Radius<<endl;
 518   cout<<" DCA Negative track to PV......: "<<fCutDCANegToPV<<endl;
 519   cout<<" DCA Positive track to PV......: "<<fCutDCAPosToPV<<endl;
 520   cout<<" DCA V0 Daughters..............: "<<fCutDCAV0Daughters<<endl;
 521   cout<<" Cosine of Pointing Angle V0...: "<<fCutV0CosPA<<endl;
 522   cout<<" Proper Lifetime cut (cm)......: "<<fCutProperLifetime<<endl;
 523   cout<<" TPC dE/dx sigmas cut (Real)...: "<<fCutTPCPIDNSigmas<<endl;
 524   cout<<" CutNSigmasForSignalExtraction.: "<<fCutNSigmasForSignalExtraction<<endl;
 525   cout<<" Least # of Crossed Rows.......: "<<fCutLeastNumberOfCrossedRows<<endl;
 526   cout<<" Least # of C. R. over findable: "<<fCutLeastNumberOfCrossedRowsOverFindable<<endl;
 527   cout<<" Daughter Track |eta| < .......: "<<fCutDaughterEta<<endl;
 528   cout<<" Competing V0 Reject. (GeV/c^2): "<<fCutCompetingV0Rejection<<endl;
 529   cout<<"--------------- File Names -----------------------------"<<endl;
 530   cout<<" Real Data File................: "<<fRealDataFile<<endl;
 531   cout<<" MC File.......................: "<<fMCDataFile<<endl;
 532   cout<<" MC File for feeddown..........: "<<fFeedDownDataFile<<endl;
 533   cout<<" Analysis output file..........: "<<fOutputDataFile<<endl;
 534   cout<<"--------------------------------------------------------"<<endl;
 535   cout<<endl;
 536
 537   //=== Real data loop 1: Acquire peak positions, widths ===
 538   //--- Preparing... ---
 539   //defining helping histogram - only used to FindBin Index!========
 540         TH1F* fHistPt           = new TH1F("fHistPt","Dummy;p_{T} (GeV/c);Counts",fptbinnumb,fptbinlimits);
 541   //Peak Position Histograms
 542         TH1F* fHistPeakPosition = new TH1F("fHistPeakPosition","Peak Position (Real);p_{T} (GeV/c);Peak Position",fptbinnumb,fptbinlimits);
 543         TH1F* fHistPeakPositionMC       = new TH1F("fHistPeakPositionMC","Peak Position (MC);p_{T} (GeV/c);Peak Position",fptbinnumb,fptbinlimits);
 544   //Signal to Noise Histograms
 545         TH1F* fHistSigToNoise     = new TH1F("fHistSigToNoise","Signal to Noise Ratio (real);p_{T} (GeV/c);Sig / Bg",fptbinnumb,fptbinlimits);
 546         TH1F* fHistSigToNoiseMC = new TH1F("fHistSigToNoiseMC","Signal to Noise Ratio (MC);p_{T} (GeV/c);Sig / Bg",fptbinnumb,fptbinlimits);
 547
 548   //Signal Extraction Range Histogram
 549   TH1F* fHistSignalExtractionRange      = new TH1F("fHistSignalExtractionRange","Sig. Ext. Range;p_{T} (GeV/c);Range",fptbinnumb,fptbinlimits);
 550
 551    //Resolution Histogram (filled with MC)
 552         TH2F* f2dHistPtResolution               = new TH2F("f2dHistPtResolution","p_{t} Resolution;p_{t} (reco);p_{t} (mc)",fptbinnumb,fptbinlimits, fptbinnumb, fptbinlimits);
 553
 554         TH2F* f2dHistPtBlur             = new TH2F("f2dHistPtBlur","p_{t} blurring matrix;p_{t} (reco);p_{t} (mc)",fptbinnumb,fptbinlimits, fptbinnumb, fptbinlimits);
 555         TH2F* f2dHistPtSharpen          = new TH2F("f2dHistPtSharpen","p_{t} sharpening matrix;p_{t} (reco);p_{t} (mc)",fptbinnumb,fptbinlimits, fptbinnumb, fptbinlimits);
 556
 557   /********************************************************
 558
 559     ---> Let's Remember the limits of the data we're analyzing!
 560     ---> Important so that we don't try signal extraction outside
 561     ---> the bundaries of available data in the Tree object.
 562
 563     From AliAnalysisTaskExtractV0.cxx
 564
 565     //Second Selection: rough 20-sigma band, parametric.
 566     //K0Short: Enough to parametrize peak broadening with linear function.
 567     Double_t UpperLimitK0Short = (5.63707e-01) + (1.14979e-02)*tree_lPt;
 568     Double_t LowerLimitK0Short = (4.30006e-01) - (1.10029e-02)*tree_lPt;
 569
 570     //Lambda: Linear (for higher pt) plus exponential (for low-pt broadening)
 571     //[0]+[1]*x+[2]*TMath::Exp(-[3]*x)
 572     Double_t UpperLimitLambda = (1.13688e+00) + (5.27838e-03)*tree_lPt + (8.42220e-02)*TMath::Exp(-(3.80595e+00)*tree_lPt);
 573     Double_t LowerLimitLambda = (1.09501e+00) - (5.23272e-03)*tree_lPt - (7.52690e-02)*TMath::Exp(-(3.46339e+00)*tree_lPt);
 574
 575   ********************************************************/
 576
 577   TF1 *fKDataUpper = new TF1("fKDataUpper","[0]+[1]*x",0,20);
 578   TF1 *fKDataLower = new TF1("fKDataLower","[0]+[1]*x",0,20);
 579   TF1 *fLDataUpper = new TF1("fLDataUpper","[0]+[1]*x+[2]*TMath::Exp([3]*x)",0,20);
 580   TF1 *fLDataLower = new TF1("fLDataLower","[0]+[1]*x+[2]*TMath::Exp([3]*x)",0,20);
 581
 582   fKDataUpper->SetParameter(0, 5.63707e-01);
 583   fKDataUpper->SetParameter(1, 1.14979e-02);
 584   fKDataLower->SetParameter(0, 4.30006e-01);
 585   fKDataLower->SetParameter(1,-1.10029e-02);
 586
 587   fLDataUpper->SetParameter(0, 1.13688e+00);
 588   fLDataUpper->SetParameter(1, 5.27838e-03);
 589   fLDataUpper->SetParameter(2, 8.42220e-02);
 590   fLDataUpper->SetParameter(3,-3.80595e+00);
 591
 592   fLDataLower->SetParameter(0, 1.09501e+00);
 593   fLDataLower->SetParameter(1,-5.23272e-03);
 594   fLDataLower->SetParameter(2,-7.52690e-02);
 595   fLDataLower->SetParameter(3,-3.46339e+00);
 596
 597   Int_t lWeAreAtBin = 0;
 598   Double_t lParticleMass = -1;          //needed for proper lifetime selection
 599   Double_t lCompetingParticleMass = -1; //needed for Competing V0 Rejection
 600   Double_t lHistoLowerBoundary = -1;
 601   Double_t lHistoUpperBoundary = -1;
 602   Long_t lHistoNBins = -1;
 603   if ( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ){
 604     lParticleMass          = 1.115683;
 605     lCompetingParticleMass = 0.4976;
 606     lHistoLowerBoundary = 1.116 - 0.1;
 607     lHistoUpperBoundary = 1.116 + 0.1;
 608     lHistoNBins = 200; //1MeV/c^2 binning
 609   }
 610   if ( fWhichParticle == "K0Short" ){
 611     lParticleMass          = 0.4976;
 612     lCompetingParticleMass = 1.115683;
 613     lHistoLowerBoundary = 0.498 - 0.15;
 614     lHistoUpperBoundary = 0.498 + 0.15;
 615     lHistoNBins = 300; //1MeV/c^2 binning
 616   }
 617         //Setting Up Base Histograms to use===============================
 618   TH1F* lHistoFullV0[100];
 619   TH1F* lHistoFullV0MC[100];
 620   TCanvas* lCanvasHistoFullV0[100];
 621   TCanvas* lCanvasHistoFullV0MC[100];
 622   TH1F *lHistResolution[100];
 623   TF1 *lHistResolutionGaussian[100];
 624   char lHistResolutionName[100];
 625         TH1F* fHistResolutionVsPt               = new TH1F("fHistResolutionVsPt","Resolution vs p_{t};p_{t} (GeV/c);#Delta(p_{t}) = p_{t}-p_{t}^{true} (GeV/c)",fptbinnumb,fptbinlimits);
 626         TH1F* fHistResolutionVsPtDivByBinWidth          = new TH1F("fHistResolutionVsPtDivByBinWidth","Resolution vs p_{t};p_{t} (GeV/c);#Delta(p_{t})/#Delta^{bin}(p_{t}) = (p_{t}-p_{t}^{true})/#Delta^{bin}(p_{t}) (GeV/c)",fptbinnumb,fptbinlimits);
 627         TH1F* fHistResolutionVsPtWithGaussians          = new TH1F("fHistResolutionVsPtWithGaussians","Resolution vs p_{t};p_{t} (GeV/c);#Delta(p_{t}) = p_{t}-p_{t}^{true} (GeV/c)",fptbinnumb,fptbinlimits);
 628         TH1F* fHistResolutionVsPtDivByBinWidthWithGaussians             = new TH1F("fHistResolutionVsPtDivByBinWidthWithGaussians","Resolution vs p_{t};p_{t} (GeV/c);#Delta(p_{t})/#Delta^{bin}(p_{t}) = (p_{t}-p_{t}^{true})/#Delta^{bin}(p_{t}) (GeV/c)",fptbinnumb,fptbinlimits);
 629   for(Int_t ibin=0;ibin<fptbinnumb;ibin++){
 630   }
 631   char histname[80];    TString bindescription = "";
 632   for(Int_t ihist=0;ihist<100;ihist++){
 633     //Histo For Real Data
 634     sprintf(histname,"lHistoFullV0%i",ihist);
 635     if(fWhichParticle == "Lambda")     bindescription="#Lambda, bin #";
 636     if(fWhichParticle == "AntiLambda") bindescription="#bar{#Lambda}, bin #";
 637     if(fWhichParticle == "K0Short")    bindescription="K^{0}_{S}, bin #";
 638                 bindescription.Append(IntToString( ihist ));
 639     if ( ihist < fptbinnumb ){
 640                   bindescription.Append(", ");
 641                   bindescription.Append(DoubleToString(fptbinlimits[ ihist ]));
 642                   bindescription.Append("-");
 643                   bindescription.Append(DoubleToString(fptbinlimits[ ihist+1 ]));
 644                   bindescription.Append("GeV/c");
 645     }
 646     lHistoFullV0[ihist]   =     new TH1F(histname,"Candidates;M(appropriate) (GeV/c^{2});Counts", lHistoNBins,lHistoLowerBoundary,lHistoUpperBoundary);
 647     lHistoFullV0[ihist]->SetTitle(bindescription);
 648     sprintf(histname,"lCanvasHistoFullV0%i",ihist);
 649     lCanvasHistoFullV0[ihist] = new TCanvas(histname, bindescription, 800, 600);
 650     lCanvasHistoFullV0[ihist] -> SetFillColor(kWhite);
 651     lCanvasHistoFullV0[ihist] -> SetLeftMargin( 0.175 );
 652     lCanvasHistoFullV0[ihist] -> SetBottomMargin( 0.175 );
 653
 654     //Histo for MC
 655     sprintf(histname,"lHistoFullV0MC%i",ihist);
 656     if(fWhichParticle == "Lambda")     bindescription="#Lambda, MC, bin #";
 657     if(fWhichParticle == "AntiLambda") bindescription="#bar{#Lambda}, MC, bin #";
 658     if(fWhichParticle == "K0Short")    bindescription="K^{0}_{S}, MC, bin #";
 659                 bindescription.Append(IntToString( ihist ));
 660     if ( ihist < fptbinnumb ){
 661                   bindescription.Append(", ");
 662                   bindescription.Append(DoubleToString(fptbinlimits[ ihist ]));
 663                   bindescription.Append("-");
 664                   bindescription.Append(DoubleToString(fptbinlimits[ ihist+1 ]));
 665                   bindescription.Append("GeV/c");
 666     }
 667     lHistoFullV0MC[ihist]   =   new TH1F(histname,"Candidates;M(appropriate) (GeV/c^{2});Counts", lHistoNBins,lHistoLowerBoundary,lHistoUpperBoundary);
 668     lHistoFullV0MC[ihist]->SetTitle(bindescription);
 669     sprintf(histname,"lCanvasHistoFullV0MC%i",ihist);
 670     lCanvasHistoFullV0MC[ihist] = new TCanvas(histname, bindescription, 800, 600);
 671     lCanvasHistoFullV0MC[ihist] -> SetFillColor(kWhite);
 672     lCanvasHistoFullV0MC[ihist] -> SetLeftMargin( 0.175 );
 673     lCanvasHistoFullV0MC[ihist] -> SetBottomMargin( 0.175 );
 674
 675     //Histo for resolution tests
 676     sprintf(lHistResolutionName,"lHistResolution%i",ihist);
 677     Long_t lNumberOfBinsResolution = 5000;
 678     if ( (fptbinlimits[ihist+1]+fptbinlimits[ihist])*0.5 > 5 )
 679       lNumberOfBinsResolution = 500;
 680     if ( (fptbinlimits[ihist+1]+fptbinlimits[ihist])*0.5 > 10 )
 681       lNumberOfBinsResolution = 50;
 682
 683     lHistResolution[ihist] = new TH1F ( lHistResolutionName,bindescription,lNumberOfBinsResolution, -5, 5); //histo with 5MeV/c precision!
 684     sprintf(lHistResolutionName,"lHistResolutionGaussian%i",ihist);
 685     lHistResolutionGaussian[ihist] = new TF1(lHistResolutionName, "[0]*TMath::Gaus(x,[1],[2])",-5,5);
 686
 687   }
 688         //================================================================
 689   cout<<endl;
 690
 691
 692   cout<<"--------------- Open Real Data File --------------------"<<endl;
 693         TFile* file = TFile::Open(fRealDataFile, "READ");
 694         file->cd("PWG2CheckLambda_PP");
 695         TList* v0list  = (TList*)file->FindObjectAny("clistV0");
 696   TTree* lTree;
 697   TH1F* fHistV0MultiplicityForTrigEvt;
 698   TH1F* fHistV0MultiplicityForSelEvtNoTPCOnly;
 699   TH1F* fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup;
 700   lTree = (TTree*)file->FindObjectAny("fTree");
 701   fHistV0MultiplicityForTrigEvt = (TH1F*)v0list->FindObject("fHistV0MultiplicityForTrigEvt");
 702   fHistV0MultiplicityForSelEvtNoTPCOnly = (TH1F*)v0list->FindObject("fHistV0MultiplicityForSelEvtNoTPCOnly");
 703   fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup = (TH1F*)v0list->FindObject("fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup");
 704
 705         Long_t lNCandidates = lTree->GetEntries();
 706         Long_t lNEvents = fHistV0MultiplicityForTrigEvt->GetEntries();
 707   Double_t lEvtFracAfterPileup =
 708     ((double)(fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup->GetEntries()))/
 709     ((double)(fHistV0MultiplicityForSelEvtNoTPCOnly->GetEntries()));
 710
 711   Double_t lNInelasticEvents = lEvtFracAfterPileup * (lNEvents / fCINT1BoverINELratio);
 712
 713   cout<<" Number of CINT1B triggers.....: "<<lNEvents<<endl;
 714         cout<<" Pileup Fraction...............: "<<1 - lEvtFracAfterPileup<<endl;
 715   cout<<" Estimated N_{INEL}............: "<<lNInelasticEvents<<endl;
 716         cout<<" V0 Candidates.................: "<<lNCandidates<<endl;
 717   cout<<"--------------------------------------------------------"<<endl;
 718
 719   Double_t lNEventsBeforePileup = 0;
 720   Double_t lNEventsAfterPileup = 0;
 721   if( fPerformMultiplicityStudy == kTRUE ){
 722     lNInelasticEvents = 0;
 723
 724     TH1F* fHistMultiplicity             = (TH1F*)v0list->FindObject("fHistMultiplicityForTrigEvt");
 725     TH1F* fHistMultiplicityBeforePileup = (TH1F*)v0list->FindObject("fHistMultiplicityNoTPCOnly");
 726     TH1F* fHistMultiplicityAfterPileup  = (TH1F*)v0list->FindObject("fHistMultiplicityNoTPCOnlyNoPileup");
 727
 728
 729
 730     for( Int_t ibin = fHistMultiplicity->FindBin(fLoMultBound); ibin<fHistMultiplicity->FindBin(fHiMultBound)+1; ibin++){
 731       lNInelasticEvents += fHistMultiplicity->GetBinContent(ibin);
 732       lNEventsBeforePileup += fHistMultiplicityBeforePileup->GetBinContent(ibin);
 733       lNEventsAfterPileup  += fHistMultiplicityAfterPileup ->GetBinContent(ibin);
 734     }
 735     cout<<" Number of events, no pileup rejection..: "<<lNInelasticEvents<<endl;
 736         cout<<" Number of events, this multiplicity....: "<<lNInelasticEvents * ( lNEventsAfterPileup / lNEventsBeforePileup ) <<endl;
 737     lNInelasticEvents = lNInelasticEvents * ( lNEventsAfterPileup / lNEventsBeforePileup );
 738     cout<<"--------------------------------------------------------"<<endl;
 739
 740   }
 741
 742         //Variable Definition=============================================
 743         //Kinematic
 744         Float_t lPt, lRap, lPtMC, lNegEta, lPosEta;
 745   //Invariant Masses
 746   Float_t lInvariantMass, lInvariantMassCompetingOne, lInvariantMassCompetingTwo; //wildcard
 747         //DCA Variables
 748         Float_t lDcaV0Daughters;
 749         Float_t lDcaPosToPrimVertex,  lDcaNegToPrimVertex;
 750         //Cosine of Pointing Angle variable
 751         Float_t lV0CosinePointingAngle;
 752   //Decay Radius and distance over total momentum
 753   Float_t lV0Radius, lDistOverTotMom;
 754   //Least Number of TPC Clusters
 755   Int_t lLeastNbrCrossedRows;
 756   Float_t lLeastNbrCrossedRowsOverFindable;
 757   //TPC dE/dx acquired with AliPIDResponse class
 758   Float_t lNSigmasPosProton,lNSigmasNegProton,lNSigmasPosPion,lNSigmasNegPion;
 759   Float_t lArmPt,lArmAlpha;
 760   //Multiplicity Variable
 761   Int_t lMultiplicity = -1;
 762         //================================================================
 763
 764         //Linking to Tree=================================================
 765
 766   //--- Base Variables ----------------------------------------------
 767         lTree->SetBranchAddress("fTreeVariablePosEta",&lPosEta);
 768         lTree->SetBranchAddress("fTreeVariableNegEta",&lNegEta);
 769         lTree->SetBranchAddress("fTreeVariablePt",&lPt);
 770   if ( fWhichParticle == "Lambda"      )  lTree->SetBranchAddress("fTreeVariableInvMassLambda",&lInvariantMass);
 771   if ( fWhichParticle == "AntiLambda"  )        lTree->SetBranchAddress("fTreeVariableInvMassAntiLambda",&lInvariantMass);
 772   if ( fWhichParticle == "K0Short"     )        lTree->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMass);
 773   if ( fWhichParticle == "Lambda"      ){ //For symmetry of computation...
 774     lTree->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMassCompetingOne);
 775     lTree->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMassCompetingTwo);
 776   }
 777   if ( fWhichParticle == "AntiLambda"  ){ //For symmetry of computation...
 778     lTree->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMassCompetingOne);
 779     lTree->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMassCompetingTwo);
 780   }
 781   if ( fWhichParticle == "K0Short"     ){
 782     lTree->SetBranchAddress("fTreeVariableInvMassLambda"    ,&lInvariantMassCompetingOne);
 783     lTree->SetBranchAddress("fTreeVariableInvMassAntiLambda",&lInvariantMassCompetingTwo);
 784   }
 785   if ( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" )
 786         lTree->SetBranchAddress("fTreeVariableRapLambda",&lRap);
 787   if ( fWhichParticle == "K0Short" )
 788         lTree->SetBranchAddress("fTreeVariableRapK0Short",&lRap);
 789         lTree->SetBranchAddress("fTreeVariableDistOverTotMom",&lDistOverTotMom);
 790         lTree->SetBranchAddress("fTreeVariableLeastNbrCrossedRows",&lLeastNbrCrossedRows);
 791         lTree->SetBranchAddress("fTreeVariableLeastRatioCrossedRowsOverFindable",&lLeastNbrCrossedRowsOverFindable);
 792   //--- TPC dEdx Variables ------------------------------------------
 793         lTree->SetBranchAddress("fTreeVariableNSigmasPosProton",&lNSigmasPosProton);
 794         lTree->SetBranchAddress("fTreeVariableNSigmasNegProton",&lNSigmasNegProton);
 795         lTree->SetBranchAddress("fTreeVariableNSigmasPosPion",&lNSigmasPosPion);
 796         lTree->SetBranchAddress("fTreeVariableNSigmasNegPion",&lNSigmasNegPion);
 797   //--- Topological selection variables -----------------------------
 798         lTree->SetBranchAddress("fTreeVariableV0Radius",&lV0Radius);
 799         lTree->SetBranchAddress("fTreeVariableDcaNegToPrimVertex",&lDcaNegToPrimVertex);
 800         lTree->SetBranchAddress("fTreeVariableDcaPosToPrimVertex",&lDcaPosToPrimVertex);
 801         lTree->SetBranchAddress("fTreeVariableDcaV0Daughters",&lDcaV0Daughters);
 802         lTree->SetBranchAddress("fTreeVariableV0CosineOfPointingAngle",&lV0CosinePointingAngle);
 803   //--- Multiplicity Variable ----------------------------------------
 804   lTree->SetBranchAddress("fTreeVariableMultiplicity",&lMultiplicity);
 805   //--- Armenteros-Podolansky ----------------------------------------
 806   lTree->SetBranchAddress("fTreeVariableAlphaV0",&lArmAlpha);
 807   lTree->SetBranchAddress("fTreeVariablePtArmV0",&lArmPt);
 808         //================================================================
 809
 810
 811   cout<<"--------------- Real Data File Loop 1 ------------------"<<endl;
 812   Long_t lOneTenthOfNCandidates = ((double)(lNCandidates) / 10. );
 813   for(Long_t icand = 0;icand<lNCandidates;icand++){
 814                 lTree->GetEntry(icand);
 815     if( icand % lOneTenthOfNCandidates == 0 )
 816         cout<<" Currently at candidate........: "<<icand<<" / "<<lNCandidates<<" ( "<<(long)(((double)(icand)/(double)(lNCandidates))*(100.+1e-3))<<"% )"<<endl;
 817                 if(TMath::Abs(lRap)<fRapidityBoundary &&
 818             TMath::Abs(lNegEta)       <= fCutDaughterEta               &&
 819             TMath::Abs(lPosEta)       <= fCutDaughterEta               &&
 820             lV0Radius                 >= fCutV0Radius                  &&
 821             lDcaNegToPrimVertex       >= fCutDCANegToPV                &&
 822             lDcaPosToPrimVertex       >= fCutDCAPosToPV                &&
 823             lDcaV0Daughters           <= fCutDCAV0Daughters            &&
 824             lV0CosinePointingAngle    >= fCutV0CosPA                   &&
 825             lParticleMass*lDistOverTotMom    <= fCutProperLifetime     &&
 826             lLeastNbrCrossedRows             >= fCutLeastNumberOfCrossedRows             &&
 827             lLeastNbrCrossedRowsOverFindable >= fCutLeastNumberOfCrossedRowsOverFindable &&
 828             TMath::Abs(lInvariantMassCompetingOne - lCompetingParticleMass) > fCutCompetingV0Rejection &&
 829             TMath::Abs(lInvariantMassCompetingTwo - lCompetingParticleMass) > fCutCompetingV0Rejection &&
 830           ( fSpecialArmenterosCutK0s == kFALSE || ( fSpecialArmenterosCutK0s == kTRUE && lArmPt*5>TMath::Abs(lArmAlpha) ) ) &&
 831           ( //official response code
 832           ( fWhichParticle == "Lambda"
 833           && TMath::Abs(lNSigmasNegPion)   <= fCutTPCPIDNSigmas
 834           && TMath::Abs(lNSigmasPosProton) <= fCutTPCPIDNSigmas) ||
 835           ( fWhichParticle == "AntiLambda"
 836           && TMath::Abs(lNSigmasPosPion)   <= fCutTPCPIDNSigmas
 837           && TMath::Abs(lNSigmasNegProton) <= fCutTPCPIDNSigmas) ||
 838           ( fWhichParticle == "K0Short"
 839           && TMath::Abs(lNSigmasPosPion)   <= fCutTPCPIDNSigmas
 840           && TMath::Abs(lNSigmasNegPion)   <= fCutTPCPIDNSigmas)
 841           ) &&
 842           ( //Multiplicity Switch
 843           fPerformMultiplicityStudy == kFALSE || //No mult selection, or...
 844           (fPerformMultiplicityStudy == kTRUE &&  //inside mult bin
 845           lMultiplicity>=fLoMultBound &&
 846           lMultiplicity<=fHiMultBound)
 847           )
 848       ){
 849                   lWeAreAtBin = fHistPt->FindBin(lPt)-1;
 850                   if(lWeAreAtBin == -1) lWeAreAtBin = 99; //UnderFlow, special treatment
 851                   lHistoFullV0[lWeAreAtBin]->Fill(lInvariantMass); //fill with specific inv mass
 852                 }
 853         }
 854   cout<<"--------------- Loop Completed -------------------------"<<endl;
 855   cout<<endl;
 856
 857   cout<<"--------------- Peak Finding (gauss+linear) ------------"<<endl;
 858  //== Variables for holding peak position, width ==============
 859   Double_t lPeakPosition[100];
 860   Double_t lPeakWidth[100];
 861   Double_t lLeftBgLeftLimit[100];
 862   Double_t lLeftBgRightLimit[100];
 863   Double_t lPeakLeftLimit[100];
 864   Double_t lPeakRightLimit[100];
 865   Double_t lRightBgLeftLimit[100];
 866   Double_t lRightBgRightLimit[100];
 867   //May be needed if bg areas are different in the future
 868   //Double_t lScaleFactor[100];
 869
 870   TLine *lLineLeftMost[100];
 871   TLine *lLineLeft[100];
 872   TLine *lLineRight[100];
 873   TLine *lLineRightMost[100];
 874
 875   TLine *lLineLeftMostMC[100];
 876   TLine *lLineLeftMC[100];
 877   TLine *lLineRightMC[100];
 878   TLine *lLineRightMostMC[100];
 879
 880   Double_t lMiddle = 0;
 881   Double_t lUpperFit = 0;
 882   Double_t lLowerFit = 0;
 883
 884   char fgausname[100];
 885   TF1 *fgausPt[100];
 886   for(Int_t ibin = 0; ibin<fptbinnumb; ibin++){
 887     cout<<"---> Peak Finding, bin #"<<ibin<<"..."<<endl;
 888     sprintf(fgausname,"fGausPt%i",ibin);
 889         if ( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ){
 890       lMiddle = 0.5 * ( fLDataLower->Eval( fHistPt->GetBinCenter(ibin+1) ) + fLDataUpper->Eval( fHistPt->GetBinCenter(ibin+1) ) );
 891       lUpperFit = lMiddle + 0.7 * ( fLDataUpper->Eval( fHistPt->GetBinCenter(ibin+1) ) - lMiddle );
 892       lLowerFit = lMiddle - 0.7 * ( lMiddle - fLDataLower->Eval( fHistPt->GetBinCenter(ibin+1) ) );
 893       fgausPt[ibin]= new TF1(fgausname,"[0]*TMath::Gaus(x,[1],[2])+[3]*x+[4]", lLowerFit, lUpperFit );
 894           fgausPt[ibin]->SetParameter(1,1.116);
 895           fgausPt[ibin]->SetParameter(2,0.0025);
 896       fgausPt[ibin]->SetParLimits(1,1,1.2);
 897       fgausPt[ibin]->SetParLimits(2,0.001,0.01);
 898     }
 899         if ( fWhichParticle == "K0Short"){
 900       fgausPt[ibin]= new TF1(fgausname,"[0]*TMath::Gaus(x,[1],[2])+[3]*x+[4]",0.498-0.06,0.498+0.060);
 901           fgausPt[ibin]->SetParameter(1,0.498);
 902           fgausPt[ibin]->SetParameter(2,0.004);
 903     }
 904     fgausPt[ibin]->SetParameter(0,lHistoFullV0[ibin]->GetMaximum() * 0.9);
 905
 906
 907           fgausPt[ibin]->SetParameter(3,0);
 908           fgausPt[ibin]->SetParameter(4,lHistoFullV0[ibin]->GetMaximum() * 0.1);
 909     lHistoFullV0[ibin]->Fit(fgausname,"QREM0");
 910     lPeakPosition[ibin] = fgausPt[ibin]->GetParameter(1);
 911     lPeakWidth[ibin] = TMath::Abs( fgausPt[ibin]->GetParameter(2) );
 912     cout<<"---> ["<<fptbinlimits[ibin]<<" - "<<fptbinlimits[ibin+1]<<" GeV/c]\tPeak at: "<<lPeakPosition[ibin]<<", sigma = "<<lPeakWidth[ibin]<<endl;
 913     //Find Corresponding Limits In this bin
 914     lLeftBgLeftLimit[ibin]  = lPeakPosition[ibin] - 2.*fCutNSigmasForSignalExtraction*lPeakWidth[ibin];
 915     lLeftBgRightLimit[ibin] = lPeakPosition[ibin] - 1.*fCutNSigmasForSignalExtraction*lPeakWidth[ibin];
 916     lPeakLeftLimit[ibin]    = lPeakPosition[ibin] - 1.*fCutNSigmasForSignalExtraction*lPeakWidth[ibin];
 917     lPeakRightLimit[ibin]   = lPeakPosition[ibin] + 1.*fCutNSigmasForSignalExtraction*lPeakWidth[ibin];
 918     lRightBgLeftLimit[ibin] = lPeakPosition[ibin] + 1.*fCutNSigmasForSignalExtraction*lPeakWidth[ibin];
 919     lRightBgRightLimit[ibin]= lPeakPosition[ibin] + 2.*fCutNSigmasForSignalExtraction*lPeakWidth[ibin];
 920     //lScaleFactor[ibin] = (lPeakRightLimit[ibin] - lPeakLeftLimit[ibin]) / ( (lLeftBgRightLimit[ibin]-lLeftBgLeftLimit[ibin]) + (lRightBgRightLimit[ibin] - lRightBgLeftLimit[ibin]) );
 921     fHistPeakPosition->SetBinContent(ibin+1, lPeakPosition[ibin]);
 922     fHistPeakPosition->SetBinError(ibin+1, lPeakWidth[ibin]);
 923     //Create Signal Extraction Range Histogram
 924     fHistSignalExtractionRange->SetBinContent(ibin+1, lPeakPosition[ibin]);
 925     fHistSignalExtractionRange->SetBinError(ibin+1, 2.*fCutNSigmasForSignalExtraction*lPeakWidth[ibin] );
 926     //Create appropriate TLine Objects for Canvases
 927     lLineLeftMost[ibin]  = new TLine( lLeftBgLeftLimit[ibin],   0, lLeftBgLeftLimit[ibin],   lHistoFullV0[ibin]->GetMaximum() * 0.95 );
 928     lLineLeft[ibin]      = new TLine( lLeftBgRightLimit[ibin],  0, lLeftBgRightLimit[ibin],  lHistoFullV0[ibin]->GetMaximum() * 0.95 );
 929     lLineRight[ibin]     = new TLine( lRightBgLeftLimit[ibin],  0, lRightBgLeftLimit[ibin],  lHistoFullV0[ibin]->GetMaximum() * 0.95 );
 930     lLineRightMost[ibin] = new TLine( lRightBgRightLimit[ibin], 0, lRightBgRightLimit[ibin], lHistoFullV0[ibin]->GetMaximum() * 0.95 );
 931
 932     //Preparing Canvas for storing...
 933     lCanvasHistoFullV0[ibin]->cd();
 934     lHistoFullV0[ibin]->Draw();
 935     lLineLeftMost[ibin]->Draw();
 936     lLineLeft[ibin]->Draw();
 937     lLineRight[ibin]->Draw();
 938     lLineRightMost[ibin]->Draw();
 939     fgausPt[ibin]->Draw("same");
 940
 941   }
 942   cout<<"--------------- Peak Finding Finished ------------------"<<endl;
 943
 944   //Defining Signal holding variables===============================
 945   Double_t lSigRealV0[100];
 946   Double_t lSigErrRealV0[100];
 947   Double_t lSigMCV0[100];
 948   Double_t lSigErrMCV0[100];
 949   //================================================================
 950         Double_t lLeftPlusRightBgV0[100];
 951         Double_t lSigPlusCenterBgV0[100];
 952         Double_t lLeftPlusRightBgV0MC[100];
 953         Double_t lSigPlusCenterBgV0MC[100];
 954   for(Long_t n1 = 0; n1<100; n1++){
 955      lLeftPlusRightBgV0[n1]=0;
 956      lSigPlusCenterBgV0[n1]=0;
 957      lLeftPlusRightBgV0MC[n1]=0;
 958      lSigPlusCenterBgV0MC[n1]=0;
 959   }
 960   //================================================================
 961   cout<<endl;
 962   cout<<"--------------- Real Data File Loop 2 ------------------"<<endl;
 963   for(Long_t icand = 0;icand<lNCandidates;icand++){
 964                 lTree->GetEntry(icand);
 965     if( icand % lOneTenthOfNCandidates == 0 )
 966         cout<<" Currently at candidate........: "<<icand<<" / "<<lNCandidates<<" ( "<<(long)(((double)(icand)/(double)(lNCandidates))*(100.+1e-3))<<"% )"<<endl;
 967                 if(TMath::Abs(lRap)<fRapidityBoundary &&
 968             TMath::Abs(lNegEta)       <= fCutDaughterEta               &&
 969             TMath::Abs(lPosEta)       <= fCutDaughterEta               &&
 970             lV0Radius                 >= fCutV0Radius                  &&
 971             lDcaNegToPrimVertex       >= fCutDCANegToPV                &&
 972             lDcaPosToPrimVertex       >= fCutDCAPosToPV                &&
 973             lDcaV0Daughters           <= fCutDCAV0Daughters            &&
 974             lV0CosinePointingAngle    >= fCutV0CosPA                   &&
 975             lParticleMass*lDistOverTotMom    <= fCutProperLifetime     &&
 976             lLeastNbrCrossedRows             >= fCutLeastNumberOfCrossedRows             &&
 977             lLeastNbrCrossedRowsOverFindable >= fCutLeastNumberOfCrossedRowsOverFindable &&
 978             TMath::Abs(lInvariantMassCompetingOne - lCompetingParticleMass) > fCutCompetingV0Rejection &&
 979             TMath::Abs(lInvariantMassCompetingTwo - lCompetingParticleMass) > fCutCompetingV0Rejection &&
 980           ( fSpecialArmenterosCutK0s == kFALSE || ( fSpecialArmenterosCutK0s == kTRUE && lArmPt*5>TMath::Abs(lArmAlpha) ) ) &&
 981           ( //official response code
 982           ( fWhichParticle == "Lambda"
 983           && TMath::Abs(lNSigmasNegPion)   <= fCutTPCPIDNSigmas
 984           && TMath::Abs(lNSigmasPosProton) <= fCutTPCPIDNSigmas) ||
 985           ( fWhichParticle == "AntiLambda"
 986           && TMath::Abs(lNSigmasPosPion)   <= fCutTPCPIDNSigmas
 987           && TMath::Abs(lNSigmasNegProton) <= fCutTPCPIDNSigmas) ||
 988           ( fWhichParticle == "K0Short"
 989           && TMath::Abs(lNSigmasPosPion)   <= fCutTPCPIDNSigmas
 990           && TMath::Abs(lNSigmasNegPion)   <= fCutTPCPIDNSigmas)
 991           ) &&
 992           ( //Multiplicity Switch
 993           fPerformMultiplicityStudy == kFALSE || //No mult selection, or...
 994           (fPerformMultiplicityStudy == kTRUE &&  //inside mult bin
 995           lMultiplicity>=fLoMultBound &&
 996           lMultiplicity<=fHiMultBound)
 997           )
 998       ){ // Start Entry Loop
 999                   lWeAreAtBin = fHistPt->FindBin(lPt)-1;
1000                   if(lWeAreAtBin == -1) lWeAreAtBin = 99; //UnderFlow, special treatment
1001
1002       //Extract left and right background areas and peak
1003       //--- Left Background Sample
1004                         if(lInvariantMass>lLeftBgLeftLimit[lWeAreAtBin]    && lInvariantMass<lLeftBgRightLimit[lWeAreAtBin]   ){ lLeftPlusRightBgV0[lWeAreAtBin]++; }
1005                         if(lInvariantMass>lRightBgLeftLimit[lWeAreAtBin]   && lInvariantMass<lRightBgRightLimit[lWeAreAtBin]  ){ lLeftPlusRightBgV0[lWeAreAtBin]++; }
1006       //--- Peak Region
1007                         if(lInvariantMass>lPeakLeftLimit[lWeAreAtBin]      && lInvariantMass<lPeakRightLimit[lWeAreAtBin]     ){ lSigPlusCenterBgV0[lWeAreAtBin]++; }
1008                 } // End Entry Loop
1009         }
1010   cout<<"--------------- Loop Completed -------------------------"<<endl;
1011   cout<<endl;
1012   cout<<"--------------- Memory Cleanup -------------------------"<<endl;
1013   if ( v0list ) {
1014     v0list->Delete();
1015     delete v0list;
1016   }
1017
1018   file->Close("R");
1019   file->Delete();
1020   delete file;
1021   cout<<endl;
1022
1023   TF1 *lfitNoise[50];
1024   TF1 *lSampleNoise[50];
1025   char lFitNameOne[50];
1026
1027   if( fFitBackgroundSwitch ){
1028     cout<<"--------------- Backgrounds: fitting with linear -------"<<endl;
1029     for(long i=0; i<fptbinnumb; i++){
1030       //Define Function to Fit Background
1031       sprintf(lFitNameOne,"lfitNoise%i",((int)(i)));
1032       //cout<<"creating fitnoise, named "<<lFitNameOne<<endl;
1033       lfitNoise[i] = new TF1(lFitNameOne, this, &AliV0Module::MyBgPol1,
1034         RoundToThousandth ( lLeftBgLeftLimit[i]   ),
1035         RoundToThousandth ( lRightBgRightLimit[i] ), 4 , "AliV0Module", "MyBgPol1");
1036       lfitNoise[i] -> FixParameter(2,RoundToThousandth ( lLeftBgRightLimit[i] ) );
1037       lfitNoise[i] -> FixParameter(3,RoundToThousandth ( lRightBgLeftLimit[i] ) );
1038       lfitNoise[i] -> SetParameter(0,lLeftPlusRightBgV0[i] * lHistoFullV0[i]->GetBinWidth(5) / (lRightBgLeftLimit[i]-lLeftBgRightLimit[i] + 1e-6 ) );
1039       lfitNoise[i] -> SetParameter(1,0);
1040       cout<<"Guessed Parameter 0 fot "<<lFitNameOne<<" to be "<<lLeftPlusRightBgV0[i] * lHistoFullV0[i]->GetBinWidth(5) / (lRightBgLeftLimit[i]-lLeftBgRightLimit[i] + 1e-6 )<<endl;
1041       sprintf(lFitNameOne,"lSampleNoise%i",((int)(i)));
1042
1043       //Define Function to Sample Background
1044       //cout<<"creating sample "<<i<<endl;
1045       lSampleNoise[i] = new TF1(lFitNameOne, this, &AliV0Module::MyBgPolToEval1,
1046         RoundToThousandth ( lLeftBgLeftLimit[i]   ),
1047         RoundToThousandth ( lRightBgRightLimit[i] ), 2 , "AliV0Module", "MyBgPolToEval1");
1048     }
1049     for(long i=0; i<fptbinnumb; i++){
1050       //cout<<"Fitting function for bin "<<i<<", get name = "<<lfitNoise[i]->GetName()<<endl;
1051       sprintf(lFitNameOne,"lfitNoise%i",((int)(i)));
1052       lHistoFullV0[i] -> Fit(lFitNameOne,"LLrie+0");
1053       lSampleNoise[i]->SetParameter(0, lfitNoise[i]->GetParameter(0) );
1054       lSampleNoise[i]->SetParameter(1, lfitNoise[i]->GetParameter(1) );
1055     }
1056     for(long i=0; i<fptbinnumb; i++){
1057       cout<<"Overriding Background info: Was "<<lLeftPlusRightBgV0[i]<<", is now "<<lSampleNoise[i]->Integral( lPeakLeftLimit[i], lPeakRightLimit[i] )/lHistoFullV0[i]->GetBinWidth(5)<<endl;
1058       lLeftPlusRightBgV0[i] = lSampleNoise[i]->Integral( lPeakLeftLimit[i], lPeakRightLimit[i] )/lHistoFullV0[i]->GetBinWidth(5);
1059     }
1060     cout<<"--------------- Fitting Finished! ----------------------"<<endl;
1061   }
1062
1063   //=============================================================
1064   // Compute Signal + Sig to noise
1065   for(Int_t ipoint = 0; ipoint<fptbinnumb; ipoint++){
1066       //Signal computation
1067       lSigRealV0[ipoint] = lSigPlusCenterBgV0[ipoint] - lLeftPlusRightBgV0[ipoint];
1068       lSigErrRealV0[ipoint] = TMath::Sqrt(lSigPlusCenterBgV0[ipoint]+lLeftPlusRightBgV0[ipoint]);
1069       //Error Formula: Equivalent to Sqrt(S+B+B) = Sqrt(S+2B)
1070
1071       //Signal-to-noise computation
1072       if( lLeftPlusRightBgV0[ipoint] != 0 ){
1073         fHistSigToNoise->SetBinContent(ipoint+1, (lSigPlusCenterBgV0[ipoint] - lLeftPlusRightBgV0[ipoint])/lLeftPlusRightBgV0[ipoint] );
1074       }else{
1075         fHistSigToNoise->SetBinContent(ipoint+1, -1) ; //-1 means: no background
1076       }
1077   }
1078   //=============================================================
1079
1080   cout<<"--------------- Extracted Signal -----------------------"<<endl;
1081   for(Int_t ipoint = 0; ipoint<fptbinnumb; ipoint++){
1082       cout<<"---> ["<<fptbinlimits[ipoint]<<" - "<<fptbinlimits[ipoint+1]<<" GeV/c]\tSignal: "<<lSigRealV0[ipoint]<<" +/- "<<lSigErrRealV0[ipoint]<<endl;
1083   }
1084   cout<<"--------------------------------------------------------"<<endl;
1085
1086   //=============================================================
1087   // Preparations for MC loop
1088   //Extra info available only in MC
1089   Int_t lPID = 0;
1090   Float_t lNegTransvMomentumMC, lPosTransvMomentumMC;
1091   Int_t lPIDPositive = 0;
1092   Int_t lPIDNegative = 0;
1093   Int_t lPIDMother = -1; Float_t lPtMother = -1;
1094   Int_t lPrimaryStatus = 0;
1095   Int_t lPrimaryStatusMother = 0;
1096   //=============================================================
1097   cout<<endl;
1098   //Prepare containers for Geant3/fluka
1099   //Will only be used and saved if we're dealing with Lambdas
1100   TH1F *lProtonMomentum[100];
1101   char lNameOne[100];
1102   for(Int_t ibin=0;ibin<100;ibin++){
1103     sprintf(lNameOne,"lProtonMomentumBin%i",ibin);
1104     lProtonMomentum[ibin] = new TH1F(lNameOne,"",800,0,20);
1105     if(fWhichParticle == "Lambda")     bindescription="#Lambda, bin #";
1106     if(fWhichParticle == "AntiLambda") bindescription="#bar{#Lambda}, bin #";
1107     if(fWhichParticle == "K0Short")    bindescription="K^{0}_{S}, bin #";
1108                 bindescription.Append(IntToString( ibin ));
1109     if ( ibin < fptbinnumb ){
1110                   bindescription.Append(", ");
1111                   bindescription.Append(DoubleToString(fptbinlimits[ ibin ]));
1112                   bindescription.Append("-");
1113                   bindescription.Append(DoubleToString(fptbinlimits[ ibin+1 ]));
1114                   bindescription.Append("GeV/c");
1115     }
1116     lProtonMomentum[ibin]->SetTitle(bindescription);
1117   }
1118
1119         //MC File Acquisition=============================================
1120   cout<<"--------------- Opening MC file ------------------------"<<endl;
1121         TFile* fileMC = TFile::Open(fMCDataFile, "READ");
1122         fileMC->cd("PWG2CheckPerformanceLambda_PP_MC");
1123         TList* v0listMC  = (TList*)fileMC->FindObjectAny("clistV0MC");
1124   TTree *lTreeMC;
1125   TH1F* fHistMultiplicityBeforeTrigSelMC;
1126   lTreeMC = (TTree*)fileMC->FindObjectAny("fTree");
1127   fHistMultiplicityBeforeTrigSelMC = (TH1F*)v0listMC->FindObject("fHistMultiplicityBeforeTrigSel");
1128
1129         Long_t lNCandidatesMC = lTreeMC->GetEntries();
1130         Long_t lNEventsMC = fHistMultiplicityBeforeTrigSelMC->GetEntries();
1131         cout<<" MC Events (before trig sel)...: "<<lNEventsMC<<endl;
1132         cout<<" Cascade MC Candidates.........: "<<lNCandidatesMC<<endl;
1133   cout<<"--------------------------------------------------------"<<endl;
1134         //================================================================
1135
1136   //Linking to Tree=================================================
1137   //--- Base Variables ----------------------------------------------
1138         lTreeMC->SetBranchAddress("fTreeVariablePosEta",&lPosEta);
1139         lTreeMC->SetBranchAddress("fTreeVariableNegEta",&lNegEta);
1140         lTreeMC->SetBranchAddress("fTreeVariablePrimaryStatus",&lPrimaryStatus);
1141         lTreeMC->SetBranchAddress("fTreeVariablePrimaryStatusMother",&lPrimaryStatusMother);
1142         lTreeMC->SetBranchAddress("fTreeVariablePt",&lPt);
1143         lTreeMC->SetBranchAddress("fTreeVariablePtXiMother",&lPtMother);
1144         lTreeMC->SetBranchAddress("fTreeVariablePtMC",&lPtMC);
1145   if ( fWhichParticle == "Lambda"      )  lTreeMC->SetBranchAddress("fTreeVariableInvMassLambda",&lInvariantMass);
1146   if ( fWhichParticle == "AntiLambda"  )        lTreeMC->SetBranchAddress("fTreeVariableInvMassAntiLambda",&lInvariantMass);
1147   if ( fWhichParticle == "K0Short"     )        lTreeMC->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMass);
1148   if ( fWhichParticle == "Lambda"      ){ //For symmetry of computation...
1149     lTreeMC->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMassCompetingOne);
1150     lTreeMC->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMassCompetingTwo);
1151   }
1152   if ( fWhichParticle == "AntiLambda"  ){ //For symmetry of computation...
1153     lTreeMC->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMassCompetingOne);
1154     lTreeMC->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMassCompetingTwo);
1155   }
1156   if ( fWhichParticle == "K0Short"     ){
1157     lTreeMC->SetBranchAddress("fTreeVariableInvMassLambda"    ,&lInvariantMassCompetingOne);
1158     lTreeMC->SetBranchAddress("fTreeVariableInvMassAntiLambda",&lInvariantMassCompetingTwo);
1159   }
1160   //if ( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" )
1161   //    lTreeMC->SetBranchAddress("fTreeVariableRapLambda",&lRap);
1162   //if ( fWhichParticle == "K0Short" )
1163   //    lTreeMC->SetBranchAddress("fTreeVariableRapK0Short",&lRap);
1164
1165   lTreeMC->SetBranchAddress("fTreeVariableRapMC",&lRap);
1166         lTreeMC->SetBranchAddress("fTreeVariableLeastNbrCrossedRows",&lLeastNbrCrossedRows);
1167         lTreeMC->SetBranchAddress("fTreeVariableLeastRatioCrossedRowsOverFindable",&lLeastNbrCrossedRowsOverFindable);
1168   //--- Topological selection variables -----------------------------
1169         lTreeMC->SetBranchAddress("fTreeVariableV0Radius",&lV0Radius);
1170         lTreeMC->SetBranchAddress("fTreeVariableDcaNegToPrimVertex",&lDcaNegToPrimVertex);
1171         lTreeMC->SetBranchAddress("fTreeVariableDcaPosToPrimVertex",&lDcaPosToPrimVertex);
1172         lTreeMC->SetBranchAddress("fTreeVariableDcaV0Daughters",&lDcaV0Daughters);
1173         lTreeMC->SetBranchAddress("fTreeVariableV0CosineOfPointingAngle",&lV0CosinePointingAngle);
1174         lTreeMC->SetBranchAddress("fTreeVariableDistOverTotMom",&lDistOverTotMom);
1175   //---- Extra  Information Available in MC only --------------------
1176         lTreeMC->SetBranchAddress("fTreeVariableNegTransvMomentumMC",&lNegTransvMomentumMC);
1177         lTreeMC->SetBranchAddress("fTreeVariablePosTransvMomentumMC",&lPosTransvMomentumMC);
1178         lTreeMC->SetBranchAddress("fTreeVariablePID",&lPID);
1179         lTreeMC->SetBranchAddress("fTreeVariablePIDMother",&lPIDMother);
1180         lTreeMC->SetBranchAddress("fTreeVariablePIDPositive",&lPIDPositive);
1181         lTreeMC->SetBranchAddress("fTreeVariablePIDNegative",&lPIDNegative);
1182   //--- Multiplicity ------------------------------------------------
1183   lTreeMC->SetBranchAddress("fTreeVariableMultiplicity",&lMultiplicity);
1184   //--- Armenteros-Podolansky ----------------------------------------
1185   lTreeMC->SetBranchAddress("fTreeVariableAlphaV0",&lArmAlpha);
1186   lTreeMC->SetBranchAddress("fTreeVariablePtArmV0",&lArmPt);
1187         //================================================================
1188
1189   //================================================================
1190   cout<<endl;
1191   cout<<"--------------- MC Data File Loop  ---------------------"<<endl;
1192   Long_t lOneTenthOfNCandidatesMC = ((double)(lNCandidatesMC) / 10. );
1193   for(Long_t icand = 0;icand<lNCandidatesMC;icand++){
1194                 lTreeMC->GetEntry(icand);
1195     if( icand % lOneTenthOfNCandidatesMC == 0 )
1196         cout<<" Currently at candidate........: "<<icand<<" / "<<lNCandidatesMC<<" ( "<<(long)(((double)(icand)/(double)(lNCandidatesMC))*(100.+1e-3))<<"% )"<<endl;
1197                 if(TMath::Abs(lRap)<fRapidityBoundary &&
1198             TMath::Abs(lNegEta)       <= fCutDaughterEta               &&
1199             TMath::Abs(lPosEta)       <= fCutDaughterEta               &&
1200             lV0Radius                 >= fCutV0Radius                  &&
1201             lDcaNegToPrimVertex       >= fCutDCANegToPV                &&
1202             lDcaPosToPrimVertex       >= fCutDCAPosToPV                &&
1203             lDcaV0Daughters           <= fCutDCAV0Daughters            &&
1204             lV0CosinePointingAngle    >= fCutV0CosPA                   &&
1205             lParticleMass*lDistOverTotMom    <= fCutProperLifetime     &&
1206             lLeastNbrCrossedRows             >= fCutLeastNumberOfCrossedRows             &&
1207             lLeastNbrCrossedRowsOverFindable >= fCutLeastNumberOfCrossedRowsOverFindable &&
1208             TMath::Abs(lInvariantMassCompetingOne - lCompetingParticleMass) > fCutCompetingV0Rejection &&
1209             TMath::Abs(lInvariantMassCompetingTwo - lCompetingParticleMass) > fCutCompetingV0Rejection &&
1210           ( fSpecialArmenterosCutK0s == kFALSE || ( fSpecialArmenterosCutK0s == kTRUE && lArmPt*5>TMath::Abs(lArmAlpha) ) ) &&
1211           ( //perfect PID association, IsPhysicalPrimary association
1212           ( fWhichParticle == "Lambda"
1213           && lPID           == 3122 //V0 is a Lambda
1214           && lPIDPositive   == 2212 //Pos Daughter is p
1215           && lPIDNegative   == -211 //Neg Daughter is pi-
1216           && lPrimaryStatus == 1
1217           ) ||
1218           ( fWhichParticle == "AntiLambda"
1219           && lPID           == -3122  //V0 is an AntiLambda
1220           && lPIDPositive   == 211    //Pos Daughter is pi+
1221           && lPIDNegative   == -2212  //Neg Daughter is antiproton
1222           && lPrimaryStatus == 1
1223           ) ||
1224           ( fWhichParticle == "K0Short"
1225           && lPID           == 310   //V0 is an AntiLambda
1226           && lPIDPositive   == 211   //Pos Daughter is pi+
1227           && lPIDNegative   == -211  //Neg Daughter is pi-
1228           && lPrimaryStatus == 1     //K0Short is a primary
1229           )
1230           ) &&
1231           ( //Multiplicity Switch
1232           fPerformMultiplicityStudy == kFALSE || //No mult selection, or...
1233           (fPerformMultiplicityStudy == kTRUE &&  //inside mult bin
1234           lMultiplicity>=fLoMultBound &&
1235           lMultiplicity<=fHiMultBound)
1236           )
1237       ){ // Start Entry Loop
1238                   lWeAreAtBin = fHistPt->FindBin(lPtMC)-1;
1239                   if(lWeAreAtBin == -1) lWeAreAtBin = 99; //UnderFlow, special treatment
1240                   lHistoFullV0MC[lWeAreAtBin]->Fill(lInvariantMass); //fill with specific inv mass
1241       f2dHistPtResolution->Fill(lPt,lPtMC);
1242       //Extract left and right background areas and peak
1243       //--- Left Background Sample
1244                         if(lInvariantMass>lLeftBgLeftLimit[lWeAreAtBin]    && lInvariantMass<lLeftBgRightLimit[lWeAreAtBin]   ){ lLeftPlusRightBgV0MC[lWeAreAtBin]++; }
1245                         if(lInvariantMass>lRightBgLeftLimit[lWeAreAtBin]   && lInvariantMass<lRightBgRightLimit[lWeAreAtBin]  ){ lLeftPlusRightBgV0MC[lWeAreAtBin]++; }
1246       //--- Peak Region
1247                         if(lInvariantMass>lPeakLeftLimit[lWeAreAtBin]      && lInvariantMass<lPeakRightLimit[lWeAreAtBin]     ){ lSigPlusCenterBgV0MC[lWeAreAtBin]++; }
1248       //--- Info may be needed for geant3/fluka
1249       if ( fWhichParticle == "Lambda"     ) lProtonMomentum[lWeAreAtBin]->Fill( lPosTransvMomentumMC );
1250       if ( fWhichParticle == "AntiLambda" ) lProtonMomentum[lWeAreAtBin]->Fill( lNegTransvMomentumMC );
1251       //--- Resolution tests
1252       lHistResolution[lWeAreAtBin]->Fill(lPt - lPtMC);
1253                 } // End Entry Loop
1254         }
1255   cout<<"--------------- Loop Completed -------------------------"<<endl;
1256   cout<<endl;
1257
1258   cout<<"--------------- X-Check Peak Finding (gauss+linear) ----"<<endl;
1259  //== Variables for holding peak position, width ==============
1260   Double_t lPeakPositionMC[100];
1261   Double_t lPeakWidthMC[100];
1262
1263   char fgausnameMC[100];
1264   TF1 *fgausPtMC[100];
1265   for(Int_t ibin = 0; ibin<fptbinnumb; ibin++){
1266     cout<<"---> Peak Finding, bin #"<<ibin<<" (perfect count = "<<lHistoFullV0MC[ibin]->GetEntries()<<")..."<<endl;
1267     sprintf(fgausnameMC,"fGausPtMC%i",ibin);
1268         if ( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ){
1269       fgausPtMC[ibin]= new TF1(fgausnameMC,"[0]*TMath::Gaus(x,[1],[2])",1.116-0.040,1.116+0.040);
1270           fgausPtMC[ibin]->SetParameter(1,1.116);
1271           fgausPtMC[ibin]->SetParameter(2,0.0025);
1272           fgausPtMC[ibin]->SetParLimits(1,1.1,1.2);
1273           fgausPtMC[ibin]->SetParLimits(2,0.001,0.02);
1274     }
1275         if ( fWhichParticle == "K0Short"){
1276       fgausPtMC[ibin]= new TF1(fgausnameMC,"[0]*TMath::Gaus(x,[1],[2])",0.498-0.06,0.498+0.060);
1277           fgausPtMC[ibin]->SetParameter(1,0.498);
1278           fgausPtMC[ibin]->SetParameter(2,0.004);
1279     }
1280     fgausPtMC[ibin]->SetParameter(0,lHistoFullV0MC[ibin]->GetMaximum() * 0.9);
1281           //fgausPtMC[ibin]->SetParameter(3,0);
1282           //fgausPtMC[ibin]->SetParameter(4,HistoFullV0MC[ibin]->GetMaximum() * 0.1);
1283     lHistoFullV0MC[ibin]->Fit(fgausnameMC,"QREM0");
1284     lPeakPositionMC[ibin] = TMath::Abs( fgausPtMC[ibin]->GetParameter(1) );
1285     lPeakWidthMC[ibin] = fgausPtMC[ibin]->GetParameter(2);
1286     cout<<"---> ["<<fptbinlimits[ibin]<<" - "<<fptbinlimits[ibin+1]<<" GeV/c]\tPeak at: "<<lPeakPositionMC[ibin]<<", sigma = "<<lPeakWidthMC[ibin]<<endl;
1287     fHistPeakPositionMC->SetBinContent(ibin+1, lPeakPositionMC[ibin]);
1288     fHistPeakPositionMC->SetBinError(ibin+1, lPeakWidthMC[ibin]);
1289
1290     //Create appropriate TLine Objects for Canvases
1291     lLineLeftMostMC[ibin]  = new TLine( lLeftBgLeftLimit[ibin],   0, lLeftBgLeftLimit[ibin],   lHistoFullV0MC[ibin]->GetMaximum() * 0.95 );
1292     lLineLeftMC[ibin]      = new TLine( lLeftBgRightLimit[ibin],  0, lLeftBgRightLimit[ibin],  lHistoFullV0MC[ibin]->GetMaximum() * 0.95 );
1293     lLineRightMC[ibin]     = new TLine( lRightBgLeftLimit[ibin],  0, lRightBgLeftLimit[ibin],  lHistoFullV0MC[ibin]->GetMaximum() * 0.95 );
1294     lLineRightMostMC[ibin] = new TLine( lRightBgRightLimit[ibin], 0, lRightBgRightLimit[ibin], lHistoFullV0MC[ibin]->GetMaximum() * 0.95 );
1295
1296     //Preparing Canvas for storing...
1297     lCanvasHistoFullV0MC[ibin]->cd();
1298     lHistoFullV0MC[ibin]->Draw();
1299     lLineLeftMostMC[ibin]->Draw();
1300     lLineLeftMC[ibin]->Draw();
1301     lLineRightMC[ibin]->Draw();
1302     lLineRightMostMC[ibin]->Draw();
1303     fgausPtMC[ibin]->Draw("same");
1304
1305   }
1306   cout<<"--------------- Peak Finding Finished (in MC) ----------"<<endl;
1307   cout<<endl;
1308
1309   TF1 *lfitNoiseMC[50];
1310   TF1 *lSampleNoiseMC[50];
1311   char lFitNameOneMC[50];
1312
1313   if( fFitBackgroundSwitch ){
1314     cout<<"--------------- Backgrounds: fitting with linear -------"<<endl;
1315     for(long i=0; i<fptbinnumb; i++){
1316       //Define Function to Fit Background
1317       sprintf(lFitNameOneMC,"lfitNoiseMC%i",(int)i);
1318       cout<<"creating fitnoise, named "<<lFitNameOneMC<<endl;
1319       lfitNoiseMC[i] = new TF1(lFitNameOneMC, this, &AliV0Module::MyBgPol1,
1320         RoundToThousandth ( lLeftBgLeftLimit[i]   ),
1321         RoundToThousandth ( lRightBgRightLimit[i] ), 4 , "AliV0Module", "MyBgPol1");
1322       lfitNoiseMC[i] -> FixParameter(2,RoundToThousandth ( lLeftBgRightLimit[i] ) );
1323       lfitNoiseMC[i] -> FixParameter(3,RoundToThousandth ( lRightBgLeftLimit[i] ) );
1324       lfitNoiseMC[i] -> SetParameter(0,lLeftPlusRightBgV0MC[i]*lHistoFullV0MC[i]->GetMaximum() / (lSigPlusCenterBgV0MC[i]+1e-6) );
1325       lfitNoiseMC[i] -> SetParameter(1,0);
1326       sprintf(lFitNameOneMC,"lSampleNoiseMC%i",(int)i);
1327
1328       //Define Function to Sample Background
1329       cout<<"creating sample "<<i<<endl;
1330       lSampleNoiseMC[i] = new TF1(lFitNameOneMC, this, &AliV0Module::MyBgPolToEval1,
1331         RoundToThousandth ( lLeftBgLeftLimit[i]   ),
1332         RoundToThousandth ( lRightBgRightLimit[i] ), 2, "AliV0Module", "MyBgPolToEval1" );
1333     }
1334     for(long i=0; i<fptbinnumb; i++){
1335       cout<<"Fitting function for bin "<<i<<", get name = "<<lfitNoiseMC[i]->GetName()<<endl;
1336       sprintf(lFitNameOneMC,"lfitNoiseMC%i",(int)i);
1337       lHistoFullV0MC[i] -> Fit(lFitNameOneMC,"LLrie+0");
1338       lSampleNoiseMC[i]->SetParameter(0, lfitNoiseMC[i]->GetParameter(0) );
1339       lSampleNoiseMC[i]->SetParameter(1, lfitNoiseMC[i]->GetParameter(1) );
1340     }
1341     for(long i=0; i<fptbinnumb; i++){
1342       cout<<"Overriding Background info: Was "<<lLeftPlusRightBgV0MC[i]<<", is now "<<lSampleNoiseMC[i]->Integral( lPeakLeftLimit[i], lPeakRightLimit[i] )/lHistoFullV0MC[i]->GetBinWidth(5)<<endl;
1343       lLeftPlusRightBgV0MC[i] = lSampleNoiseMC[i]->Integral( lPeakLeftLimit[i], lPeakRightLimit[i] )/lHistoFullV0MC[i]->GetBinWidth(5);
1344     }
1345     cout<<"--------------- Fitting Finished! ----------------------"<<endl;
1346   }
1347
1348
1349   //=============================================================
1350   // Compute Signal
1351   for(Int_t ipoint = 0; ipoint<fptbinnumb; ipoint++){
1352       lSigMCV0[ipoint] = lSigPlusCenterBgV0MC[ipoint] - lLeftPlusRightBgV0MC[ipoint];
1353       lSigErrMCV0[ipoint] = TMath::Sqrt(lSigPlusCenterBgV0MC[ipoint]+lLeftPlusRightBgV0MC[ipoint]);
1354       //Error Formula: Equivalent to Sqrt(S+B+B) = Sqrt(S+2B)
1355
1356       //Signal-to-noise computation
1357       if( lLeftPlusRightBgV0MC[ipoint] != 0 ){
1358         fHistSigToNoiseMC->SetBinContent(ipoint+1, (lSigPlusCenterBgV0MC[ipoint] - lLeftPlusRightBgV0MC[ipoint])/lLeftPlusRightBgV0MC[ipoint] );
1359       }else{
1360         fHistSigToNoiseMC->SetBinContent(ipoint+1, -1) ; //-1 means: no background
1361       }
1362   }
1363   //=============================================================
1364   cout<<"--------------- Extracted Signal (MC) ------------------"<<endl;
1365   for(Int_t ipoint = 0; ipoint<fptbinnumb; ipoint++){
1366       cout<<"---> ["<<fptbinlimits[ipoint]<<" - "<<fptbinlimits[ipoint+1]<<" GeV/c]\tSignal: "<<lSigMCV0[ipoint]<<" +/- "<<lSigErrMCV0[ipoint]<<endl;
1367   }
1368   cout<<"--------------------------------------------------------"<<endl;
1369   cout<<endl;
1370
1371
1372   cout<<"--------------- Process Generated V0s ------------------"<<endl;
1373         //Filling histogram with original MC particles====================
1374   TH3F* f3dHistGenPtVsYVsMultV0 = 0x0;
1375           if(fWhichParticle == "Lambda")
1376       f3dHistGenPtVsYVsMultV0                   = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultLambda");
1377           if(fWhichParticle == "AntiLambda")
1378       f3dHistGenPtVsYVsMultV0                   = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultAntiLambda");
1379           if(fWhichParticle == "K0Short")
1380       f3dHistGenPtVsYVsMultV0                   = (TH3F*)v0listMC->FindObject("f3dHistPrimRawPtVsYVsMultK0Short");
1381
1382         TH1D *fHistDummyV0 = 0x0;
1383         if( fPerformMultiplicityStudy == kFALSE ){
1384           fHistDummyV0 = f3dHistGenPtVsYVsMultV0->ProjectionX("fHistDummyV0",
1385                   f3dHistGenPtVsYVsMultV0->GetYaxis()->FindBin(-fRapidityBoundary+1e-2),  //avoid taking previous bin
1386                   f3dHistGenPtVsYVsMultV0->GetYaxis()->FindBin(+fRapidityBoundary-1e-2), //avoid taking next bin
1387                   f3dHistGenPtVsYVsMultV0->GetZaxis()->FindBin(-1), //Not interested in Multiplicity so far, integrate
1388                   f3dHistGenPtVsYVsMultV0->GetZaxis()->FindBin(101) //Not interested in Multiplicity so far, integrate
1389           );
1390   }
1391         if( fPerformMultiplicityStudy == kTRUE ){
1392           fHistDummyV0 = f3dHistGenPtVsYVsMultV0->ProjectionX("fHistDummyV0",
1393                   f3dHistGenPtVsYVsMultV0->GetYaxis()->FindBin(-fRapidityBoundary+1e-2),  //avoid taking previous bin
1394                   f3dHistGenPtVsYVsMultV0->GetYaxis()->FindBin(+fRapidityBoundary-1e-2), //avoid taking next bin
1395                   f3dHistGenPtVsYVsMultV0->GetZaxis()->FindBin(fLoMultBound), //Interested in Multiplicity!
1396                   f3dHistGenPtVsYVsMultV0->GetZaxis()->FindBin(fHiMultBound)  //Interested in Multiplicity!
1397           );
1398   }
1399
1400         TH1F *fHistMCCountbyptV0        = new TH1F("fHistMCCountbyptV0","V0 MC count;p_{T} (GeV/c);Counts",fptbinnumb,fptbinlimits);
1401
1402         Double_t temppt;
1403         for(long i = 1; i<fHistDummyV0->GetNbinsX()+1;i++){
1404                 temppt = fHistDummyV0->GetXaxis()->GetBinCenter(i);
1405                 for(long filling = 0; filling<fHistDummyV0->GetBinContent(i);filling++){
1406                         fHistMCCountbyptV0->Fill(temppt);
1407                 }
1408         }
1409   cout<<"--------------- Generated V0 Dump ----------------------"<<endl;
1410   for(Int_t ipoint = 0; ipoint<fptbinnumb; ipoint++){
1411       cout<<"---> ["<<fptbinlimits[ipoint]<<" - "<<fptbinlimits[ipoint+1]<<" GeV/c]\tSignal: "<<fHistMCCountbyptV0->GetBinContent(ipoint+1)<<" +/- "<<TMath::Sqrt(fHistMCCountbyptV0->GetBinContent(ipoint+1))<<endl;
1412   }
1413   cout<<"--------------------------------------------------------"<<endl;
1414
1415   //---------------------------------------------------------------
1416   //extremely curious test!!! experimental only, FIXME
1417   cout<<"Resolution tests ongoing - Experimental! "<<endl;
1418   cout<<"--------------------------------------------------------"<<endl;
1419   cout<<" Usual tests: (Pt - PtMC) histograms "<<endl;
1420   cout<<" Fitting gaussians..."<<endl;
1421   for(Int_t ibin=0;ibin<fptbinnumb;ibin++){
1422     lHistResolutionGaussian[ibin]->SetParameter(0, lHistResolution[ibin]->GetMaximum() );
1423     lHistResolutionGaussian[ibin]->SetParameter(1, lHistResolution[ibin]->GetMean()    );
1424     lHistResolutionGaussian[ibin]->SetParameter(2, lHistResolution[ibin]->GetRMS()     );
1425     sprintf(lHistResolutionName,"lHistResolutionGaussian%i",ibin);
1426     lHistResolution[ibin]->Fit(lHistResolutionName,"IREM0");
1427   }
1428
1429   Double_t lTotSigMCPtMC[100];
1430   Double_t lTotSigMCPtMCXCheck[100];
1431   Double_t lTotSigMCPtReco[100];
1432   Double_t lTotSigMCPtRecoXCheck[100];
1433   for(Int_t ibin=0;ibin<100;ibin++){
1434     lTotSigMCPtMC[ibin] = 0;
1435     lTotSigMCPtReco[ibin] = 0;
1436     lTotSigMCPtRecoXCheck[ibin] = 0;
1437   }
1438   for(Int_t ibin=0; ibin<fptbinnumb; ibin++){   //reco pt
1439     for(Int_t jbin=0; jbin<fptbinnumb; jbin++){ //mc pt
1440       lTotSigMCPtMC[jbin] += ((double)f2dHistPtResolution->GetBinContent(ibin+1,jbin+1));
1441       lTotSigMCPtReco[ibin] += ((double)f2dHistPtResolution->GetBinContent(ibin+1,jbin+1));
1442     }
1443   }
1444
1445
1446   TMatrixD fResolutionMatrix(fptbinnumb,fptbinnumb);
1447   for(Int_t ibin=0; ibin<fptbinnumb; ibin++){   //reco pt
1448     for(Int_t jbin=0; jbin<fptbinnumb; jbin++){ //mc pt
1449       fResolutionMatrix[ibin][jbin] = ((double)f2dHistPtResolution->GetBinContent(ibin+1,jbin+1)) / lTotSigMCPtMC[jbin] ;
1450       f2dHistPtBlur->SetBinContent(ibin+1,jbin+1,(float)fResolutionMatrix[ibin][jbin]);
1451     }
1452   }
1453   fResolutionMatrix.Print();
1454   cout<<"---> Debug Test 1: Check if fResolutionMatrix * lTotSigMCPtMC is lTotSigMCPtReco"<<endl;
1455   cout<<"---> Multiplying... "<<endl;
1456   for(Int_t ibin=0; ibin<fptbinnumb; ibin++){   //reco pt
1457     for(Int_t jbin=0; jbin<fptbinnumb; jbin++){ //mc pt
1458       //if resolution is perfect: lTempSigReal[ibin] = lSigRealV0[ibin]
1459       // (as fResolution is delta[ibin][jbin])
1460       lTotSigMCPtRecoXCheck[ibin] += lTotSigMCPtMC[jbin]*fResolutionMatrix[ibin][jbin];
1461     }
1462   }
1463
1464   cout<<"---> Compare for Cross-check!..."<<endl;
1465   for(Int_t ibin=0; ibin<fptbinnumb; ibin++){
1466     cout<<"---> bin #"<<ibin<<"\t"<<lTotSigMCPtRecoXCheck[ibin]<<"\t"<<lTotSigMCPtReco[ibin]<<endl;
1467   }
1468
1469
1470   cout<<"Resolution matrix has determinant: "<<fResolutionMatrix.Determinant()<<endl;
1471   if(TMath::Abs(fResolutionMatrix.Determinant()) > 1e-5){
1472     cout<<"Determinant is different from zero enough to do inversion."<<endl;
1473     cout<<"Will, however, suppress data too far from diagonal!"<<endl;
1474     cout<<"Truncation: meant to ensure better numerical stability"<<endl;
1475
1476     for(Int_t ibin=0; ibin<fptbinnumb; ibin++){   //reco pt
1477       for(Int_t jbin=0; jbin<fptbinnumb; jbin++){ //mc pt
1478         if ( TMath::Abs(ibin-jbin) > 1 ){
1479           fResolutionMatrix[ibin][jbin] = 0;
1480         }
1481       }
1482     }
1483
1484     fResolutionMatrix.Invert();
1485     fResolutionMatrix.Print();
1486   }
1487
1488   fResolutionMatrix.Print();
1489   cout<<"---> Debug Test 2: Check if fResolutionMatrix^-1 * lTotSigMCPtReco is lTotSigMCPtMC"<<endl;
1490   cout<<"---> AKA \"the big test\"!"<<endl;
1491   cout<<"---> Multiplying... "<<endl;
1492   for(Int_t ibin=0; ibin<fptbinnumb; ibin++){   //reco pt
1493     for(Int_t jbin=0; jbin<fptbinnumb; jbin++){ //mc pt
1494       //if resolution is perfect: lTempSigReal[ibin] = lSigRealV0[ibin]
1495       // (as fResolution is delta[ibin][jbin])
1496       lTotSigMCPtMCXCheck[ibin] += lTotSigMCPtReco[jbin]*fResolutionMatrix[ibin][jbin];
1497     }
1498   }
1499   cout<<"---> Compare for Cross-check!..."<<endl;
1500   for(Int_t ibin=0; ibin<fptbinnumb; ibin++){
1501     cout<<"---> bin #"<<ibin<<"\t"<<lTotSigMCPtMCXCheck[ibin]<<"\t"<<lTotSigMCPtMC[ibin]<<endl;
1502   }
1503
1504
1505   for(Int_t ibin=0; ibin<fptbinnumb; ibin++){   //reco pt
1506     for(Int_t jbin=0; jbin<fptbinnumb; jbin++){ //mc pt
1507       f2dHistPtSharpen->SetBinContent(ibin+1,jbin+1,(float)fResolutionMatrix[ibin][jbin]);
1508     }
1509   }
1510
1511   Double_t lTempSigRealV0[100];
1512   for(Int_t ibin=0; ibin<fptbinnumb; ibin++){
1513     lTempSigRealV0[ibin] = 0;
1514   }
1515   //To use: de-convolute by doing matrix multiplication with lSigRealV0[]...
1516   for(Int_t ibin=0; ibin<fptbinnumb; ibin++){   //reco pt
1517     for(Int_t jbin=0; jbin<fptbinnumb; jbin++){ //mc pt
1518       //if resolution is perfect: lTempSigReal[ibin] = lSigRealV0[ibin]
1519       // (as fResolution is delta[ibin][jbin])
1520       lTempSigRealV0[ibin] += lSigRealV0[jbin]*fResolutionMatrix[ibin][jbin];
1521     }
1522   }
1523
1524   cout<<"Compare considering resolutions"<<endl;
1525   for(Int_t ibin=0; ibin<fptbinnumb; ibin++){
1526     cout<<"bin #"<<ibin<<"\t"<<lTempSigRealV0[ibin]<<"\t"<<lSigRealV0[ibin]<<endl;
1527   }
1528   cout<<endl;
1529   cout<<"inspect ratio for trends..."<<endl;
1530   for(Int_t ibin=0; ibin<fptbinnumb; ibin++){
1531     cout<<"bin #"<<ibin<<"\t"<<lTempSigRealV0[ibin]/lSigRealV0[ibin]<<endl;
1532   }
1533   //---------------------------------------------------------------
1534
1535   //=============================================================
1536   // Compute Efficiency
1537   Double_t lEfficiency[100];
1538   Double_t lEfficiencyError[100];
1539   for(Int_t ipoint = 0; ipoint<fptbinnumb; ipoint++){
1540     lEfficiency[ipoint] = lSigMCV0[ipoint] / fHistMCCountbyptV0->GetBinContent(ipoint+1);
1541     lEfficiencyError[ipoint] = ErrorInRatio(
1542       lSigMCV0[ipoint],
1543       lSigErrMCV0[ipoint],
1544       fHistMCCountbyptV0->GetBinContent(ipoint+1),
1545       TMath::Sqrt(fHistMCCountbyptV0->GetBinContent(ipoint+1))
1546     );
1547   }
1548   //=============================================================
1549   cout<<endl;
1550   cout<<"--------------- Memory Cleanup -------------------------"<<endl;
1551   fHistDummyV0->Delete();
1552   fHistMCCountbyptV0->Delete();
1553   v0listMC->Delete();
1554   delete v0listMC;
1555   fileMC->Close("R");
1556   fileMC->Delete();
1557   delete fileMC;
1558   cout<<endl;
1559   cout<<"--------------- Pure Efficiency Numbers ----------------"<<endl;
1560         TH1F* fHistPureEfficiency               = new TH1F("fHistPureEfficiency","Pure Efficiency;p_{T} (GeV/c);Pure Efficiency",fptbinnumb,fptbinlimits);
1561         TH1F* fHistEfficiency               = new TH1F("fHistEfficiency","Efficiency;p_{T} (GeV/c);Efficiency",fptbinnumb,fptbinlimits);
1562   if( fWhichParticle == "K0Short"      ) fHistPureEfficiency->SetTitle("K^{0}_{S} Efficiency");
1563   if( fWhichParticle == "Lambda"       ) fHistPureEfficiency->SetTitle("#Lambda Efficiency (no G3/F corr.)");
1564   if( fWhichParticle == "AntiLambda"   ) fHistPureEfficiency->SetTitle("#bar{#Lambda} Efficiency (no G3/F corr.)");
1565
1566   if( fWhichParticle == "Lambda"       ) fHistEfficiency->SetTitle("#Lambda Efficiency (with G3/F corr.)");
1567   if( fWhichParticle == "AntiLambda"   ) fHistEfficiency->SetTitle("#bar{#Lambda} Efficiency (with G3/F corr.)");
1568
1569   for(Int_t ipoint = 0; ipoint<fptbinnumb; ipoint++){
1570     cout<<"---> ["<<fptbinlimits[ipoint]<<" - "<<fptbinlimits[ipoint+1]<<" GeV/c]\tEff.: "<<lEfficiency[ipoint]<<" +/- "<<lEfficiencyError[ipoint]<<endl;
1571     fHistPureEfficiency->SetBinContent(ipoint+1, lEfficiency[ipoint]);
1572     fHistPureEfficiency->SetBinError(ipoint+1, lEfficiencyError[ipoint]);
1573   }
1574   cout<<"--------------------------------------------------------"<<endl;
1575
1576
1577   cout<<endl;
1578   //If it's Lambda or AntiLambda, do Geant3/fluka correction
1579   //Keep fit function outside 'if' scope
1580   TF1 *fitGeant3FlukaCorr = 0x0;
1581   TH1D* h1dPanosCorrections =0x0;
1582   if( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ){
1583     cout<<"--------------- Geant3/Fluka Correction ----------------"<<endl;
1584
1585     //Credit for the original version of this code: Antonin Maire (thanks!)
1586
1587                         TFile *fileCorrGeanT3FlukaPanos = new TFile( "aCorrectionForPPbarCrossSections_FromPanos.root", "READ");
1588       TH2D* h2dCorrForCrossSectionProtons      = (TH2D*)  ( fileCorrGeanT3FlukaPanos->Get("gHistCorrectionForCrossSectionProtons")      );
1589       TH2D* h2dCorrCrossSectionAntiProtons     = (TH2D*)  ( fileCorrGeanT3FlukaPanos->Get("gHistCorrectionForCrossSectionAntiProtons")  );
1590                         //Find bins...
1591
1592       //TH2D* h2dAsMCPtBaryonVsPtCasc  = 0x0;
1593       TH2D* h2dCorrCrossSection      = 0x0;
1594
1595                         if( fWhichParticle == "Lambda" )     h2dCorrCrossSection   = h2dCorrForCrossSectionProtons;
1596       if( fWhichParticle == "AntiLambda" ) h2dCorrCrossSection   = h2dCorrCrossSectionAntiProtons;
1597
1598       Int_t myFirstYBinPanos = h2dCorrCrossSection ->GetXaxis()->FindBin(-fRapidityBoundary);
1599       Int_t myLastYBinPanos  = h2dCorrCrossSection ->GetXaxis()->FindBin(+fRapidityBoundary-1e-3);
1600
1601                         cout<<"Check limiting bins: "<<endl;
1602         Printf(" - Bin %d = %f to %f", myFirstYBinPanos, h2dCorrCrossSection ->GetXaxis()->GetBinLowEdge(myFirstYBinPanos), h2dCorrCrossSection ->GetXaxis()->GetBinUpEdge(myFirstYBinPanos) );
1603         Printf(" - Bin %d = %f to %f", myLastYBinPanos,  h2dCorrCrossSection ->GetXaxis()->GetBinLowEdge(myLastYBinPanos ), h2dCorrCrossSection ->GetXaxis()->GetBinUpEdge(myLastYBinPanos) );
1604
1605         // B.1.a.
1606         h1dPanosCorrections = h2dCorrCrossSection->ProjectionY("h1dPanosCorrections", myFirstYBinPanos, myLastYBinPanos, "e" ); // ~mid-rapidity
1607         h1dPanosCorrections->SetDirectory(0); //means manual cleanup later
1608         h1dPanosCorrections->Scale( (Double_t) (1.0/(myLastYBinPanos - myFirstYBinPanos +1.0)) );
1609         h1dPanosCorrections->SetYTitle("#epsilon_{GEANT3} / #epsilon_{FLUKA}");
1610
1611         //TF1 *fitGeant3FlukaCorr = 0x0; declared outside this scope
1612         if(fWhichParticle == "Lambda") fitGeant3FlukaCorr = new TF1("fitGeant3FlukaCorr", this, &AliV0Module::MyGeant3FlukaCorrectionForProtons,     0.25, 1.1, 3, "AliV0Module","MyGeant3FlukaCorrectionForProtons");
1613         if(fWhichParticle == "AntiLambda") fitGeant3FlukaCorr = new TF1("fitGeant3FlukaCorr", this, &AliV0Module::MyGeant3FlukaCorrectionForAntiProtons, 0.25, 1.1, 4, "AliV0Module", "MyGeant3FlukaCorrectionForAntiProtons");
1614
1615         h1dPanosCorrections->Fit("fitGeant3FlukaCorr","rime+0"); // Chi2
1616
1617         Printf("Test fit function...");
1618         Printf(" - fitGeant3FlukaCorr for pt =  .25 GeV/c : %f", fitGeant3FlukaCorr->Eval( .25) );
1619         Printf(" - fitGeant3FlukaCorr for pt =  .5  GeV/c : %f", fitGeant3FlukaCorr->Eval( .5) );
1620         Printf(" - fitGeant3FlukaCorr for pt =  1   GeV/c : %f", fitGeant3FlukaCorr->Eval( 1.0) );
1621         Printf(" - fitGeant3FlukaCorr for pt =  2   GeV/c : %f", fitGeant3FlukaCorr->Eval( 2.0) );
1622         Printf(" - fitGeant3FlukaCorr for pt =  5   GeV/c : %f", fitGeant3FlukaCorr->Eval( 5.0) );
1623         Printf(" - fitGeant3FlukaCorr for pt =  7   GeV/c : %f", fitGeant3FlukaCorr->Eval( 7.0) );
1624         Printf(" - fitGeant3FlukaCorr for pt =  8   GeV/c : %f", fitGeant3FlukaCorr->Eval( 8.0) );
1625         Printf(" - fitGeant3FlukaCorr for pt = 10   GeV/c : %f", fitGeant3FlukaCorr->Eval(10.0) );
1626     cout<<"--------------- Embedding G3/F in Efficiencies ---------"<<endl;
1627     for(Int_t ipoint = 0; ipoint<fptbinnumb; ipoint++){
1628                   lEfficiency[ipoint]         /= fitGeant3FlukaCorr->Eval( lProtonMomentum[ipoint]->GetMean() ) ;
1629       lEfficiencyError[ipoint]    /= fitGeant3FlukaCorr->Eval( lProtonMomentum[ipoint]->GetMean() ) ;
1630     }
1631     cout<<"--------------- G3/F Corrected Efficiencies ------------"<<endl;
1632     for(Int_t ipoint = 0; ipoint<fptbinnumb; ipoint++){
1633       cout<<"---> ["<<fptbinlimits[ipoint]<<" - "<<fptbinlimits[ipoint+1]<<" GeV/c]\tEff.: "<<lEfficiency[ipoint]<<" +/- "<<lEfficiencyError[ipoint]<<endl;
1634       fHistEfficiency->SetBinContent(ipoint+1, lEfficiency[ipoint]);
1635       fHistEfficiency->SetBinError(ipoint+1, lEfficiencyError[ipoint]);
1636     }
1637     cout<<"--------------------------------------------------------"<<endl;
1638     cout<<endl;
1639
1640     fileCorrGeanT3FlukaPanos->Close();
1641     fileCorrGeanT3FlukaPanos->Delete();
1642     delete fileCorrGeanT3FlukaPanos;
1643
1644   } //end Geant3/fluka correction if
1645
1646   //Declaring outside Feeddown scope to keep in memory
1647   TH2F *f2dFeedDownMatrix = 0x0;
1648   TH2F *f2dFeedDownEfficiency = 0x0;
1649   TH2F *f2dFeedDownEfficiencyGFCorrected = 0x0;
1650   TH1F *fHistFeeddownSubtraction = 0x0;
1651
1652   //=========================================================================
1653   // --- Feeddown correction section
1654   if( fFDSwitch != "NoFD" && fWhichParticle != "K0Short"){
1655     cout<<"--------------- Feeddown Correction --------------------"<<endl;
1656           //MC File Acquisition=============================================
1657     cout<<"--------------- Opening MC file for Feeddown -----------"<<endl;
1658           TFile* fileMCFD = TFile::Open(fFeedDownDataFile, "READ");
1659           fileMCFD->cd("PWG2CheckPerformanceLambda_PP_MC");
1660           TList* v0listMCFD  = (TList*)fileMCFD->FindObjectAny("clistV0MC");
1661     TTree* lTreeMCFD;
1662     TH1F* fHistMultiplicityBeforeTrigSelMCFD;
1663     lTreeMCFD = (TTree*)fileMCFD->FindObjectAny("fTree");
1664     fHistMultiplicityBeforeTrigSelMCFD = (TH1F*)v0listMCFD->FindObject("fHistMultiplicityBeforeTrigSel");
1665
1666           Long_t lNCandidatesMCFD = lTreeMCFD->GetEntries();
1667           Long_t lNEventsMCFD = fHistMultiplicityBeforeTrigSelMCFD->GetEntries();
1668           cout<<" MC Events (before trig sel)...: "<<lNEventsMCFD<<endl;
1669           cout<<" Cascade MC Candidates.........: "<<lNCandidatesMCFD<<endl;
1670     cout<<"--------------------------------------------------------"<<endl;
1671           //================================================================
1672
1673     //Linking to Tree=================================================
1674     //--- Base Variables ----------------------------------------------
1675           lTreeMCFD->SetBranchAddress("fTreeVariablePosEta",&lPosEta);
1676           lTreeMCFD->SetBranchAddress("fTreeVariableNegEta",&lNegEta);
1677           lTreeMCFD->SetBranchAddress("fTreeVariablePrimaryStatus",&lPrimaryStatus);
1678           lTreeMCFD->SetBranchAddress("fTreeVariablePrimaryStatusMother",&lPrimaryStatusMother);
1679           lTreeMCFD->SetBranchAddress("fTreeVariablePt",&lPt);
1680           lTreeMCFD->SetBranchAddress("fTreeVariablePtMC",&lPtMC);
1681         lTreeMCFD->SetBranchAddress("fTreeVariablePtXiMother",&lPtMother);
1682     if ( fWhichParticle == "Lambda"      )  lTreeMCFD->SetBranchAddress("fTreeVariableInvMassLambda",&lInvariantMass);
1683     if ( fWhichParticle == "AntiLambda"  )      lTreeMCFD->SetBranchAddress("fTreeVariableInvMassAntiLambda",&lInvariantMass);
1684     if ( fWhichParticle == "K0Short"     )      lTreeMCFD->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMass);
1685     if ( fWhichParticle == "Lambda"      ){ //For symmetry of computation...
1686       lTreeMCFD->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMassCompetingOne);
1687       lTreeMCFD->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMassCompetingTwo);
1688     }
1689     if ( fWhichParticle == "AntiLambda"  ){ //For symmetry of computation...
1690       lTreeMCFD->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMassCompetingOne);
1691       lTreeMCFD->SetBranchAddress("fTreeVariableInvMassK0s",&lInvariantMassCompetingTwo);
1692     }
1693     if ( fWhichParticle == "K0Short"     ){
1694       lTreeMCFD->SetBranchAddress("fTreeVariableInvMassLambda"    ,&lInvariantMassCompetingOne);
1695       lTreeMCFD->SetBranchAddress("fTreeVariableInvMassAntiLambda",&lInvariantMassCompetingTwo);
1696     }
1697     //if ( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" )
1698     //  lTreeMCFD->SetBranchAddress("fTreeVariableRapLambda",&lRap);
1699     //if ( fWhichParticle == "K0Short" )
1700     //  lTreeMCFD->SetBranchAddress("fTreeVariableRapK0Short",&lRap);
1701
1702     lTreeMCFD->SetBranchAddress("fTreeVariableRapMC",&lRap);
1703           lTreeMCFD->SetBranchAddress("fTreeVariableLeastNbrCrossedRows",&lLeastNbrCrossedRows);
1704           lTreeMCFD->SetBranchAddress("fTreeVariableLeastRatioCrossedRowsOverFindable",&lLeastNbrCrossedRowsOverFindable);
1705     //--- Topological selection variables -----------------------------
1706           lTreeMCFD->SetBranchAddress("fTreeVariableV0Radius",&lV0Radius);
1707           lTreeMCFD->SetBranchAddress("fTreeVariableDcaNegToPrimVertex",&lDcaNegToPrimVertex);
1708           lTreeMCFD->SetBranchAddress("fTreeVariableDcaPosToPrimVertex",&lDcaPosToPrimVertex);
1709           lTreeMCFD->SetBranchAddress("fTreeVariableDcaV0Daughters",&lDcaV0Daughters);
1710           lTreeMCFD->SetBranchAddress("fTreeVariableV0CosineOfPointingAngle",&lV0CosinePointingAngle);
1711           lTreeMCFD->SetBranchAddress("fTreeVariableDistOverTotMom",&lDistOverTotMom);
1712     //---- Extra  Information Available in MC only --------------------
1713           lTreeMCFD->SetBranchAddress("fTreeVariableNegTransvMomentumMC",&lNegTransvMomentumMC);
1714           lTreeMCFD->SetBranchAddress("fTreeVariablePosTransvMomentumMC",&lPosTransvMomentumMC);
1715           lTreeMCFD->SetBranchAddress("fTreeVariablePID",&lPID);
1716           lTreeMCFD->SetBranchAddress("fTreeVariablePIDMother",&lPIDMother);
1717           lTreeMCFD->SetBranchAddress("fTreeVariablePIDPositive",&lPIDPositive);
1718           lTreeMCFD->SetBranchAddress("fTreeVariablePIDNegative",&lPIDNegative);
1719     //---- Multiplicity info -------------------------------------------
1720     lTreeMCFD->SetBranchAddress("fTreeVariableMultiplicity",&lMultiplicity);
1721     //--- Armenteros-Podolansky ----------------------------------------
1722     lTreeMCFD->SetBranchAddress("fTreeVariableAlphaV0",&lArmAlpha);
1723     lTreeMCFD->SetBranchAddress("fTreeVariablePtArmV0",&lArmPt);
1724           //================================================================
1725
1726     //================================================================
1727     // Define Feeddown matrix
1728     Double_t lFeedDownMatrix [100][100];
1729     //       FeedDownMatrix [Lambda Bin][Xi Bin];
1730     for(Int_t ilb = 0; ilb<100; ilb++){
1731       for(Int_t ixb = 0; ixb<100; ixb++){
1732         lFeedDownMatrix[ilb][ixb]=0;
1733       }
1734     }
1735     //Define Xi Binning
1736           Double_t xibinlimits[] = {
1737        0.00,  0.20,  0.40,  0.60,  0.80,  0.90,
1738        1.00,  1.10,  1.20,  1.30,  1.40,  1.50,
1739        1.70,  1.90,  2.20,  2.60,  3.10,  3.90,
1740        4.90,  6.00,  7.20,  8.50 ,10.00, 12.00 };
1741     Long_t xibinnumb = sizeof(xibinlimits)/sizeof(Double_t) - 1;
1742
1743         TH1F* fHistPtXiReference =
1744       new TH1F("fHistPtXiReference","#Xi candidates uncorrected p_{T};p_{T} (GeV/c);Counts",xibinnumb,xibinlimits);
1745
1746     //Feeddown: will use a double-coordinate system:
1747     // ( lambda bin, xi bin ) all the time!
1748     lWeAreAtBin = 0;         // Lambda bin
1749     Int_t lWeAreAtXiBin = 0; // Xi Bin
1750
1751     cout<<"--------------- MC Data File Loop, Feeddown  -----------"<<endl;
1752     Long_t lOneTenthOfNCandidatesMCFD = ((double)(lNCandidatesMCFD) / 10. );
1753     for(Long_t icand = 0;icand<lNCandidatesMCFD;icand++){
1754                   lTreeMCFD->GetEntry(icand);
1755       if( icand % lOneTenthOfNCandidatesMCFD == 0 )
1756         cout<<" Currently at candidate........: "<<icand<<" / "<<lNCandidatesMCFD<<" ( "<<(long)(((double)(icand)/(double)(lNCandidatesMCFD))*(100.+1e-3))<<"% )"<<endl;
1757                   if(TMath::Abs(lRap)<fRapidityBoundary &&
1758               TMath::Abs(lNegEta)       <= fCutDaughterEta               &&
1759               TMath::Abs(lPosEta)       <= fCutDaughterEta               &&
1760               lV0Radius                 >= fCutV0Radius                  &&
1761               lDcaNegToPrimVertex       >= fCutDCANegToPV                &&
1762               lDcaPosToPrimVertex       >= fCutDCAPosToPV                &&
1763               lDcaV0Daughters           <= fCutDCAV0Daughters            &&
1764               lV0CosinePointingAngle    >= fCutV0CosPA                   &&
1765               lParticleMass*lDistOverTotMom    <= fCutProperLifetime     &&
1766               lLeastNbrCrossedRows             >= fCutLeastNumberOfCrossedRows             &&
1767               lLeastNbrCrossedRowsOverFindable >= fCutLeastNumberOfCrossedRowsOverFindable &&
1768               TMath::Abs(lInvariantMassCompetingOne - lCompetingParticleMass) > fCutCompetingV0Rejection &&
1769               TMath::Abs(lInvariantMassCompetingTwo - lCompetingParticleMass) > fCutCompetingV0Rejection &&
1770             ( fSpecialArmenterosCutK0s == kFALSE || ( fSpecialArmenterosCutK0s == kTRUE && lArmPt*5>TMath::Abs(lArmAlpha) ) ) &&
1771             ( //perfect PID association, IsPhysicalPrimary association
1772             ( fWhichParticle == "Lambda"
1773             && lPID           == 3122 //V0 is a Lambda
1774             && lPIDPositive   == 2212 //Pos Daughter is p
1775             && lPIDNegative   == -211 //Neg Daughter is pi-
1776             && (lPIDMother     == 3312 || (lPIDMother == 3322 && fFDSwitch=="UseMCRatio") )
1777             && lPrimaryStatusMother == 1 //Xi is actually a primary (should matter little)
1778             ) ||
1779             ( fWhichParticle == "AntiLambda"
1780             && lPID           == -3122  //V0 is an AntiLambda
1781             && lPIDPositive   == 211    //Pos Daughter is pi+
1782             && lPIDNegative   == -2212  //Neg Daughter is antiproton
1783             && (lPIDMother     ==-3312 || (lPIDMother ==-3322 && fFDSwitch=="UseMCRatio") )
1784             && lPrimaryStatusMother == 1 //Xi is actually a primary (should matter little)
1785             )
1786             ) &&
1787             ( //Multiplicity Switch
1788             fPerformMultiplicityStudy == kFALSE || //No mult selection, or...
1789             (fPerformMultiplicityStudy == kTRUE &&  //inside mult bin
1790             lMultiplicity>=fLoMultBound &&
1791             lMultiplicity<=fHiMultBound)
1792             )
1793         ){ // Start Entry Loop
1794                     lWeAreAtBin   = fHistPt->FindBin(lPtMC)-1;
1795                     if(lWeAreAtBin == -1) lWeAreAtBin   = 99; //UnderFlow, special treatment
1796         lWeAreAtXiBin = fHistPtXiReference->FindBin(lPtMother)-1;
1797                     if(lWeAreAtXiBin == -1) lWeAreAtXiBin = 99; //UnderFlow, special treatment
1798         //Populate Feeddown Matrix
1799         //cout<<" Populate at coordinates "<<WeAreAtBin<<" vs "<<WeAreAtXiBin<<" ....."<<endl;
1800         lFeedDownMatrix[lWeAreAtBin][lWeAreAtXiBin]++;
1801                   } // End Entry Loop
1802           }
1803     cout<<"--------------- Loop Completed -------------------------"<<endl;
1804     cout<<endl;
1805
1806
1807           //Filling histogram with original MC particles====================
1808     TH3F* f3dHistGenPtVsYVsMultV0FeedDown = 0x0;
1809
1810             if (fWhichParticle == "Lambda" )
1811         f3dHistGenPtVsYVsMultV0FeedDown = (TH3F*)v0listMCFD->FindObject("f3dHistGenPtVsYVsMultXiMinus");
1812             if (fWhichParticle == "AntiLambda" )
1813         f3dHistGenPtVsYVsMultV0FeedDown = (TH3F*)v0listMCFD->FindObject("f3dHistGenPtVsYVsMultXiPlus");
1814
1815           TH1D *fHistDummyV0FeedDown = f3dHistGenPtVsYVsMultV0FeedDown->ProjectionX("fHistDummyV0FeedDown",
1816                   f3dHistGenPtVsYVsMultV0FeedDown->GetYaxis()->FindBin(-fRapidityBoundary+1e-2),  //avoid taking previous bin
1817                   f3dHistGenPtVsYVsMultV0FeedDown->GetYaxis()->FindBin(+fRapidityBoundary-1e-2), //avoid taking next bin
1818                   f3dHistGenPtVsYVsMultV0FeedDown->GetZaxis()->FindBin(-1),
1819                   f3dHistGenPtVsYVsMultV0FeedDown->GetZaxis()->FindBin(101)
1820           );
1821
1822           TH1F *fHistMCCountbyptXiFeedDown      = new TH1F("fHistMCCountbyptXiFeedDown","#Xi MC count;p_{T} (GeV/c);Counts",xibinnumb,xibinlimits);
1823           //Double_t temppt, y; --- declared before
1824     for(long i = 1; i<fHistDummyV0FeedDown->GetNbinsX()+1;i++){
1825       temppt = fHistDummyV0FeedDown->GetXaxis()->GetBinCenter(i);
1826       for(long filling = 0; filling<fHistDummyV0FeedDown->GetBinContent(i);filling++){
1827         fHistMCCountbyptXiFeedDown->Fill(temppt);
1828       }
1829     }
1830     if(fWhichParticle == "Lambda")
1831       cout<<"--------------- Generated Xi- Dump ---------------------"<<endl;
1832     if(fWhichParticle == "AntiLambda")
1833       cout<<"--------------- Generated Xi+ Dump ---------------------"<<endl;
1834     for(Int_t ipoint = 0; ipoint<xibinnumb; ipoint++){
1835         cout<<"---> ["<<xibinlimits[ipoint]<<" - "<<xibinlimits[ipoint+1]<<" GeV/c]\tSignal: "<<fHistMCCountbyptXiFeedDown->GetBinContent(ipoint+1)<<" +/- "<<TMath::Sqrt(fHistMCCountbyptXiFeedDown->GetBinContent(ipoint+1))<<endl;
1836     }
1837     cout<<"--------------------------------------------------------"<<endl;
1838     cout<<endl;
1839     cout<<"Computing actual feeddown matrix"<<endl;
1840     Double_t lFeedDownEfficiency[100][100];
1841     Double_t lFeedDownEfficiencyGFCorrected[100][100];
1842     Double_t lFeedDownEfficiencyError[100][100];
1843     Double_t lFeedDownEfficiencyGFCorrectedError[100][100];
1844
1845     for(Int_t ilb = 0; ilb<fptbinnumb; ilb++){
1846       for(Int_t ixb = 0; ixb<xibinnumb; ixb++){
1847         if( fHistMCCountbyptXiFeedDown->GetBinContent(ixb+1)!= 0 ){ //avoid div by zero
1848           lFeedDownEfficiency[ilb][ixb]=((double)(lFeedDownMatrix[ilb][ixb])) / ((double)( fHistMCCountbyptXiFeedDown->GetBinContent(ixb+1) )) ;
1849           lFeedDownEfficiencyError[ilb][ixb]=ErrorInRatio(
1850             ((double)(lFeedDownMatrix[ilb][ixb])),
1851             TMath::Sqrt((double)(lFeedDownMatrix[ilb][ixb])),
1852             ((double)( fHistMCCountbyptXiFeedDown->GetBinContent(ixb+1) )),
1853             TMath::Sqrt((double)( fHistMCCountbyptXiFeedDown->GetBinContent(ixb+1) ))
1854           );
1855         }else{
1856           lFeedDownEfficiency[ilb][ixb] = 0;
1857           lFeedDownEfficiencyError[ilb][ixb] = 0;
1858         }
1859       }
1860     }
1861
1862     //Beware: the feeddown correction efficiency is computed with MC,
1863     //Which has the geant3/fluka problem. Thus, we actually geant3/fluka
1864     //correct the Feeddown efficiencies, too...
1865         for(Int_t ipoint = 0; ipoint<fptbinnumb; ipoint++){
1866       for(Int_t ixb = 0;ixb<xibinnumb; ixb++){
1867         lFeedDownEfficiencyGFCorrected[ipoint][ixb] = lFeedDownEfficiency[ipoint][ixb]/fitGeant3FlukaCorr->Eval( lProtonMomentum[ipoint]->GetMean() ) ;
1868         lFeedDownEfficiencyGFCorrectedError[ipoint][ixb] = lFeedDownEfficiencyError[ipoint][ixb]/fitGeant3FlukaCorr->Eval( lProtonMomentum[ipoint]->GetMean() ) ;
1869       }
1870     }
1871
1872     //Create FD TH2Fs for storing
1873     f2dFeedDownMatrix = new TH2F("f2dFeedDownMatrix","",fptbinnumb,fptbinlimits,xibinnumb,xibinlimits);
1874     f2dFeedDownEfficiency = new TH2F("f2dFeedDownEfficiency","",fptbinnumb,fptbinlimits,xibinnumb,xibinlimits);
1875     f2dFeedDownEfficiencyGFCorrected = new TH2F("f2dFeedDownEfficiencyGFCorrected","",fptbinnumb,fptbinlimits,xibinnumb,xibinlimits);
1876     f2dFeedDownMatrix->SetDirectory(0);
1877     f2dFeedDownEfficiency->SetDirectory(0);
1878     f2dFeedDownEfficiencyGFCorrected->SetDirectory(0);
1879     for(Int_t ilb = 0; ilb<fptbinnumb; ilb++){
1880       for(Int_t ixb = 0; ixb<xibinnumb; ixb++){
1881         f2dFeedDownMatrix->SetBinContent(ilb+1,ixb+1,lFeedDownMatrix[ilb][ixb]);
1882         f2dFeedDownEfficiency->SetBinContent(ilb+1,ixb+1,lFeedDownEfficiency[ilb][ixb]);
1883         f2dFeedDownEfficiency->SetBinError(ilb+1,ixb+1,lFeedDownEfficiencyError[ilb][ixb]);
1884         f2dFeedDownEfficiencyGFCorrected->SetBinContent(ilb+1,ixb+1,lFeedDownEfficiencyGFCorrected[ilb][ixb]);
1885         f2dFeedDownEfficiencyGFCorrected->SetBinError(ilb+1,ixb+1,lFeedDownEfficiencyGFCorrectedError[ilb][ixb]);
1886       }
1887     }
1888
1889     Double_t lProducedXi[100];
1890
1891           TF1 *lLevyFitXiFeedDown = new TF1("LevyFitXiFeedDown",this ,&AliV0Module::MyLevyPtXi, 0.0,15,3, "AliV0Module","MyLevyPXi");
1892
1893     //Set Parameters Adequately, as in paper
1894     //FIXME: These are the 7TeV Xi- parameters!!!
1895     //FIXME: Use points vs fit function (should matter little if Xi- fit is good)
1896
1897     if(fWhichParticle == "Lambda"){
1898       lLevyFitXiFeedDown->SetParameter(0, 7.98e-03);
1899       lLevyFitXiFeedDown->SetParameter(1, 3.4429e-1);
1900       lLevyFitXiFeedDown->SetParameter(2, 1.0787e+1);
1901     }
1902     if(fWhichParticle == "AntiLambda"){
1903       lLevyFitXiFeedDown->SetParameter(0, 7.79e-03);
1904       lLevyFitXiFeedDown->SetParameter(1, 3.3934e-1);
1905       lLevyFitXiFeedDown->SetParameter(2, 1.0428e+1);
1906     }
1907     //If you want me to double charged Xi feeddown, I'll do it here
1908     if( fFDSwitch == "DoubleChargedXi" ){
1909     lLevyFitXiFeedDown->SetParameter(0, lLevyFitXiFeedDown->GetParameter(0)*2 );
1910     }
1911     for(Int_t ixb = 0; ixb<xibinnumb; ixb++){
1912       lProducedXi[ixb] = lNInelasticEvents * lLevyFitXiFeedDown->Integral(xibinlimits[ixb],xibinlimits[ixb+1]);
1913     }
1914     if(fWhichParticle == "Lambda")
1915       cout<<"--------------- Generated Xi- Dump (real-corrected) ----"<<endl;
1916     if(fWhichParticle == "AntiLambda")
1917       cout<<"--------------- Generated Xi+ Dump (real-corrected) ----"<<endl;
1918     for(Int_t ixb = 0; ixb<xibinnumb; ixb++){
1919         cout<<"Xi bin "<<ixb<<"\t"<<lProducedXi[ixb]<<endl;
1920     }
1921     cout<<"--------------------------------------------------------"<<endl;
1922     cout<<endl;
1923
1924     //=== Computing actual feeddown numbers... ===
1925     Double_t lFeedDownToSubtract[100];
1926     Double_t lFeedDownToSubtractError[100];
1927     for(Int_t ilb = 0; ilb < fptbinnumb; ilb++){
1928       lFeedDownToSubtract[ilb] = 0;
1929       lFeedDownToSubtractError[ilb] = 0;
1930       for(Int_t ixb = 0; ixb < xibinnumb; ixb++){
1931         lFeedDownToSubtract[ilb] += lProducedXi[ixb]*lFeedDownEfficiencyGFCorrected[ilb][ixb];
1932         lFeedDownToSubtractError[ilb] += TMath::Power( (lProducedXi[ixb]*lFeedDownEfficiencyGFCorrectedError[ilb][ixb]), 2 );
1933       }
1934       lFeedDownToSubtractError[ilb] = TMath::Sqrt(lFeedDownToSubtractError[ilb]);
1935     }
1936         fHistFeeddownSubtraction                = new TH1F("fHistFeeddownSubtraction","FD Subtraction;p_{T} (GeV/c);Ratio removed",fptbinnumb,fptbinlimits);
1937     fHistFeeddownSubtraction->SetDirectory(0);
1938     for(Int_t ilb = 0; ilb<fptbinnumb; ilb++){
1939       fHistFeeddownSubtraction->SetBinContent( ilb+1 , ((double)(lFeedDownToSubtract[ilb])) / ((double)(lSigRealV0[ilb]) ) );
1940       fHistFeeddownSubtraction->SetBinError  ( ilb+1 ,
1941         ErrorInRatio(
1942           ((double)(lFeedDownToSubtract[ilb])),
1943           ((double)(lFeedDownToSubtractError[ilb])),
1944           ((double)(lSigRealV0[ilb]) ),
1945           ((double)(lSigErrRealV0[ilb]) ) )
1946         );
1947     }
1948
1949
1950
1951     cout<<"--------------- FD Subtraction Fraction ----------------"<<endl;
1952     for(Int_t ilb = 0; ilb<fptbinnumb; ilb++){
1953       cout<<"---> ["<<fptbinlimits[ilb]<<" - "<<fptbinlimits[ilb+1]<<" GeV/c]\t"<<fHistFeeddownSubtraction->GetBinContent(ilb+1)<<"\t+/-\t"<<fHistFeeddownSubtraction->GetBinError(ilb+1)<<endl;
1954     }
1955     cout<<"--------------------------------------------------------"<<endl;
1956     cout<<endl;
1957     cout<<"--------------------------------------------------------"<<endl;
1958     cout<<" Performing Actual Correction... "<<endl;
1959     for(Int_t ilb = 0; ilb<fptbinnumb; ilb++){
1960       lSigRealV0   [ilb] = lSigRealV0[ilb] - lFeedDownToSubtract[ilb];
1961       lSigErrRealV0[ilb] = TMath::Sqrt( TMath::Power(lSigErrRealV0[ilb],2) + TMath::Power(lFeedDownToSubtractError[ilb],2) );
1962     }
1963     cout<<"--------------------------------------------------------"<<endl;
1964     cout<<endl;
1965     cout<<"--------------- Memory Cleanup -------------------------"<<endl;
1966           lLevyFitXiFeedDown->Delete();
1967     v0listMCFD->Delete();
1968     delete v0listMCFD;
1969     fileMCFD->Close("R");
1970     fileMCFD->Delete();
1971     delete fileMCFD;
1972     cout<<endl;
1973   }//--- End Feeddown Correction section
1974   //=========================================================================
1975
1976   //=========================================================================
1977   //---> At this stage, everthing's just ready for the actual spectrum
1978   //---> computation to occur!
1979   Double_t lSpectrum[100];
1980   Double_t lSpectrumError[100];
1981
1982   for(Int_t ibin=0;ibin<fptbinnumb;ibin++){
1983     lSpectrum[ibin] = lSigRealV0[ibin] / lEfficiency [ibin];
1984     lSpectrumError[ibin] = ErrorInRatio(
1985       lSigRealV0    [ibin],
1986       lSigErrRealV0 [ibin],
1987       lEfficiency      [ibin],
1988       lEfficiencyError [ibin]);
1989   }
1990
1991   //Divide by: Bin Width, Rapidity Window, N_{Inel}
1992   for(Int_t ibin=0;ibin<fptbinnumb;ibin++){
1993     lSpectrum[ibin] /= (fptbinlimits[ibin+1]-fptbinlimits[ibin]);
1994     lSpectrum[ibin] /= lNInelasticEvents;
1995     lSpectrum[ibin] /= 2*fRapidityBoundary;
1996     lSpectrumError[ibin] /= (fptbinlimits[ibin+1]-fptbinlimits[ibin]);
1997     lSpectrumError[ibin] /= lNInelasticEvents;
1998     lSpectrumError[ibin] /= 2*fRapidityBoundary;
1999   }
2000
2001   TH1F* fHistPtLambda = new TH1F("fHistPtLambda","#Lambda Corrected Spectrum;p_{T} (GeV/c);1/N_{INEL} #frac{d^{2}N}{dydp_{t}}",fptbinnumb,fptbinlimits);
2002   TH1F* fHistPtAntiLambda = new TH1F("fHistPtAntiLambda","#bar{#Lambda} Corrected Spectrum;p_{T} (GeV/c);1/N_{INEL} #frac{d^{2}N}{dydp_{t}}",fptbinnumb,fptbinlimits);
2003   TH1F* fHistPtK0Short = new TH1F("fHistPtK0Short","K^{0}_{S} Corrected Spectrum;p_{T} (GeV/c);1/N_{INEL} #frac{d^{2}N}{dydp_{t}}",fptbinnumb,fptbinlimits);
2004
2005   //Copy to Histogram
2006   for(Int_t ibin=0;ibin<fptbinnumb;ibin++){
2007     if(fWhichParticle == "Lambda" ){
2008       fHistPtLambda->SetBinContent( ibin+1, lSpectrum[ibin]      );
2009       fHistPtLambda->SetBinError  ( ibin+1, lSpectrumError[ibin] );
2010     }
2011     if(fWhichParticle == "AntiLambda" ){
2012       fHistPtAntiLambda->SetBinContent( ibin+1, lSpectrum[ibin]      );
2013       fHistPtAntiLambda->SetBinError  ( ibin+1, lSpectrumError[ibin] );
2014     }
2015     if(fWhichParticle == "K0Short" ){
2016       fHistPtK0Short->SetBinContent( ibin+1, lSpectrum[ibin]      );
2017       fHistPtK0Short->SetBinError  ( ibin+1, lSpectrumError[ibin] );
2018     }
2019   }
2020
2021
2022   //=========================================================================
2023
2024   cout<<"--------------- Result Output --------------------------"<<endl;
2025   cout<<" ---> Writing information to "<<fOutputDataFile<<endl;
2026   // Open an output file
2027   TFile* lResultsFile = TFile::Open(fOutputDataFile, "RECREATE");
2028   if (!lResultsFile || !lResultsFile->IsOpen()){
2029     cout<<"Error! Couldn't open file!"<<endl;
2030     return;
2031   }
2032
2033   //Preparing Signal Extraction Range Canvas
2034   TCanvas *cSigExtRange = new TCanvas("cSigExtRange","Extraction Range",900,600);
2035   cSigExtRange->SetFillColor(kWhite);
2036   cSigExtRange->SetLeftMargin(0.17);
2037   cSigExtRange->SetRightMargin(0.17);
2038   cSigExtRange->cd();
2039   fHistSignalExtractionRange->SetFillColor(18);
2040   fHistSignalExtractionRange->SetMarkerStyle(20);
2041   if( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ){
2042     fHistSignalExtractionRange->GetYaxis()->SetRangeUser(1.115-0.08,1.115+0.08);
2043   }
2044   if( fWhichParticle == "K0Short" ){
2045     fHistSignalExtractionRange->GetYaxis()->SetRangeUser(0.498-0.15,0.498+0.15);
2046   }
2047   fHistSignalExtractionRange->Draw("E2");
2048   if( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ){
2049     fLDataUpper->Draw("same");
2050     fLDataLower->Draw("same");
2051   }
2052   if( fWhichParticle == "K0Short" ){
2053     fKDataUpper->Draw("same");
2054     fKDataLower->Draw("same");
2055   }
2056   //Saving Invariant Mass Plots (real data)
2057   lResultsFile->cd();
2058   TDirectoryFile *lInvMassReal = new TDirectoryFile("lInvMassReal","Invariant Mass Plots (Real Data)");
2059   lInvMassReal->cd();
2060   cSigExtRange->Write();
2061   for(Int_t ibin = 0; ibin<fptbinnumb; ibin++) {
2062     lCanvasHistoFullV0[ibin] -> Write();
2063   }
2064
2065   TDirectoryFile *lInvMassRealRawData = new TDirectoryFile("lInvMassRealRawData","Objects for Inv Mass Plots (Real Data)");
2066   lInvMassRealRawData->cd();
2067   fHistSignalExtractionRange->Write();
2068   fHistPeakPosition->Write();
2069   fHistSigToNoise->Write();
2070   for(Int_t ibin = 0; ibin<fptbinnumb; ibin++) {
2071     lHistoFullV0[ibin] -> Write();
2072     fgausPt[ibin]         -> Write();
2073     if( fFitBackgroundSwitch ) lfitNoise[ibin] -> Write();
2074   }
2075   if( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ){
2076     fLDataUpper->Write();
2077     fLDataLower->Write();
2078   }
2079   if( fWhichParticle == "K0Short" ){
2080     fKDataUpper->Write();
2081     fKDataLower->Write();
2082   }
2083
2084   //Saving Invariant Mass Plots (MC)
2085   lResultsFile->cd();
2086   TDirectoryFile *lInvMassMC = new TDirectoryFile("lInvMassMC","Invariant Mass Plots (Monte Carlo)");
2087   lInvMassMC->cd();
2088
2089   for(Int_t ibin = 0; ibin<fptbinnumb; ibin++) {
2090     lCanvasHistoFullV0MC[ibin] -> Write();
2091   }
2092
2093   TDirectoryFile *lInvMassMCRawData = new TDirectoryFile("lInvMassMCRawData","Objects for Inv Mass Plots (MC)");
2094   lInvMassMCRawData->cd();
2095   fHistPeakPositionMC->Write();
2096   fHistSigToNoiseMC->Write();
2097    f2dHistPtResolution->Write();
2098    f2dHistPtBlur->Write();
2099    f2dHistPtSharpen->Write();
2100   for(Int_t ibin = 0; ibin<fptbinnumb; ibin++) {
2101     lHistoFullV0MC[ibin] ->Write();
2102     fgausPtMC[ibin]         ->Write();
2103     if( fFitBackgroundSwitch ) lfitNoiseMC[ibin] -> Write();
2104   }
2105
2106   //Saving Geant3/Fluka Correction Data (MC)
2107   if( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ){
2108     lResultsFile->cd();
2109     TDirectoryFile *lGFCorrection = new TDirectoryFile("lGFCorrection","Geant3/Fluka Correction Histograms");
2110     lGFCorrection->cd();
2111     h1dPanosCorrections->Write();
2112     fitGeant3FlukaCorr->Write();
2113     for(Int_t ibin = 0; ibin<fptbinnumb; ibin++) {
2114       lProtonMomentum[ibin]  ->Write();
2115     }
2116   }
2117
2118   //Saving Feeddown Correction information, if needed
2119   if( (fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda")&&fFDSwitch!="NoFD" ){
2120     lResultsFile->cd();
2121     TDirectoryFile *lFeeddown = new TDirectoryFile("lFeeddown","Feeddown Subtraction Information");
2122     lFeeddown->cd();
2123     f2dFeedDownMatrix->Write();
2124     f2dFeedDownEfficiency->Write();
2125     f2dFeedDownEfficiencyGFCorrected->Write();
2126     fHistFeeddownSubtraction->Write();
2127   }
2128
2129   //Saving Resolution Information
2130
2131   //preparing...
2132   for(Int_t ibin=0; ibin<fptbinnumb; ibin++){
2133     fHistResolutionVsPt->SetBinContent(ibin+1,lHistResolution[ibin]->GetMean());
2134     fHistResolutionVsPt->SetBinError(ibin+1,lHistResolution[ibin]->GetRMS());
2135     fHistResolutionVsPtDivByBinWidth->SetBinContent(ibin+1,lHistResolution[ibin]->GetMean() / (fptbinlimits[ibin+1]-fptbinlimits[ibin]) );
2136     fHistResolutionVsPtDivByBinWidth->SetBinError(ibin+1,lHistResolution[ibin]->GetRMS() / (fptbinlimits[ibin+1]-fptbinlimits[ibin]) );
2137
2138     fHistResolutionVsPtWithGaussians->SetBinContent(ibin+1,lHistResolutionGaussian[ibin]->GetParameter(1));
2139     fHistResolutionVsPtWithGaussians->SetBinError(ibin+1,lHistResolutionGaussian[ibin]->GetParameter(2));
2140     fHistResolutionVsPtDivByBinWidthWithGaussians->SetBinContent(ibin+1,lHistResolutionGaussian[ibin]->GetParameter(1) / (fptbinlimits[ibin+1]-fptbinlimits[ibin]) );
2141     fHistResolutionVsPtDivByBinWidthWithGaussians->SetBinError(ibin+1,lHistResolutionGaussian[ibin]->GetParameter(2) / (fptbinlimits[ibin+1]-fptbinlimits[ibin]) );
2142   }
2143
2144
2145   lResultsFile->cd();
2146   TDirectoryFile *lDirResolution = new TDirectoryFile("lInfoResolution","Resolution Information");
2147   lDirResolution->cd();
2148   fHistResolutionVsPt->Write();
2149   fHistResolutionVsPtDivByBinWidth->Write();
2150   fHistResolutionVsPtWithGaussians->Write();
2151   fHistResolutionVsPtDivByBinWidthWithGaussians->Write();
2152   for(Int_t ibin = 0; ibin<fptbinnumb; ibin++) {
2153     lHistResolution[ibin] ->Write();
2154     lHistResolutionGaussian[ibin]->Write();
2155   }
2156
2157
2158   lResultsFile->cd();
2159
2160   if( fWhichParticle == "K0Short") fHistPureEfficiency->Write();
2161   if( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ){
2162     fHistPureEfficiency->Write();
2163     fHistEfficiency->Write();
2164   }
2165   if(fWhichParticle == "Lambda" ){
2166     fHistPtLambda     ->  Write();
2167   }
2168   if(fWhichParticle == "AntiLambda" ){
2169     fHistPtAntiLambda -> Write();
2170   }
2171   if(fWhichParticle == "K0Short" ){
2172     fHistPtK0Short    -> Write();
2173   }
2174
2175   lResultsFile->Write();
2176   lResultsFile->Close();
2177   delete lResultsFile;
2178
2179   //================================================
2180   //Manual Cleanup of all created Histograms
2181   //================================================
2182   // needed if you want to re-run the whole thing without
2183   // memory leaks (systematics, etc)
2184
2185   //switch on if you want large amounts of printout
2186   Bool_t lDebugCleaningProcess = kFALSE;
2187
2188   if (lDebugCleaningProcess) cout<<"fHistPt->Delete()"<<endl;
2189   fHistPt->Delete();
2190   if (lDebugCleaningProcess) cout<<"fHistPeakPosition->Delete()"<<endl;
2191   fHistPeakPosition->Delete();
2192   if (lDebugCleaningProcess) cout<<"fHistPeakPositionMC->Delete()"<<endl;
2193   fHistPeakPositionMC->Delete();
2194   if (lDebugCleaningProcess) cout<<"fHistSigToNoise->Delete()"<<endl;
2195   fHistSigToNoise->Delete();
2196   if (lDebugCleaningProcess) cout<<"fHistSigToNoiseMC->Delete()"<<endl;
2197   fHistSigToNoiseMC->Delete();
2198   if (lDebugCleaningProcess) cout<<"fHistSignalExtractionRange->Delete()"<<endl;
2199   fHistSignalExtractionRange->Delete();
2200   if (lDebugCleaningProcess) cout<<"fKData*->Delete()"<<endl;
2201   fKDataUpper->Delete();
2202   fKDataLower->Delete();
2203   fLDataUpper->Delete();
2204   fLDataLower->Delete();
2205    if(lDebugCleaningProcess) cout<<"f2dHistPtResolution->Delete()"<<endl;
2206    f2dHistPtResolution->Delete();
2207
2208   if(lDebugCleaningProcess) cout<<"lfitNoise*[*]->Delete(); lSampleNoise*[*]->Delete()"<<endl;
2209   if( fFitBackgroundSwitch ){
2210     for(long i=0; i<fptbinnumb; i++){
2211       lfitNoise[i]    -> Delete();
2212       lSampleNoise[i] -> Delete();
2213       lfitNoiseMC[i]    -> Delete();
2214       lSampleNoiseMC[i] -> Delete();
2215     }
2216   }
2217
2218
2219   //pt-by-pt histos
2220   fHistResolutionVsPt->Delete();
2221   fHistResolutionVsPtDivByBinWidth->Delete();
2222   fHistResolutionVsPtWithGaussians->Delete();
2223   fHistResolutionVsPtDivByBinWidthWithGaussians->Delete();
2224   f2dHistPtBlur->Delete();
2225   f2dHistPtSharpen->Delete();
2226   if (lDebugCleaningProcess) cout<<"lHistoFullV0*[*]->Delete()"<<endl;
2227   for(Int_t ihist=0;ihist<100;ihist++){
2228     lHistoFullV0[ihist]->Delete();
2229     lHistoFullV0MC[ihist]->Delete();
2230     lHistResolution[ihist]->Delete();
2231     lHistResolutionGaussian[ihist]->Delete();
2232   }
2233
2234   if (lDebugCleaningProcess) cout<<"lLine*[*]->Delete()"<<endl;
2235   //gaussian fit functions, drawing lines
2236   for(Int_t ibin = 0; ibin<fptbinnumb; ibin++){
2237     lLineLeftMost[ibin] ->Delete();
2238     lLineLeft[ibin]     ->Delete();
2239     lLineRight[ibin]    ->Delete();
2240     lLineRightMost[ibin]->Delete();
2241     lLineLeftMostMC[ibin] ->Delete();
2242     lLineLeftMC[ibin]     ->Delete();
2243     lLineRightMC[ibin]    ->Delete();
2244     lLineRightMostMC[ibin]->Delete();
2245
2246         if ( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ){
2247       fgausPt[ibin]->Delete();
2248       fgausPtMC[ibin]->Delete();
2249     }
2250         if ( fWhichParticle == "K0Short"){
2251       fgausPt[ibin]->Delete();
2252       fgausPtMC[ibin]->Delete();
2253     }
2254   }
2255   if (lDebugCleaningProcess) cout<<"lCanvasHistoFullV0*[*]->Delete()"<<endl;
2256   for(Int_t ihist=0;ihist<100;ihist++){
2257     lCanvasHistoFullV0[ihist]->Close();
2258     lCanvasHistoFullV0MC[ihist]->Close();
2259     delete lCanvasHistoFullV0[ihist];
2260     delete lCanvasHistoFullV0MC[ihist];
2261   }
2262
2263   if (lDebugCleaningProcess) cout<<"cSigExtRange->Delete()"<<endl;
2264   cSigExtRange->Close();
2265   delete cSigExtRange;
2266
2267   if (lDebugCleaningProcess) cout<<"lProtonMomentum[*]->Delete()"<<endl;
2268   //histograms for G3/F correction
2269   for(Int_t ibin=0;ibin<100;ibin++){
2270     lProtonMomentum[ibin]->Delete();
2271   }
2272
2273   if (lDebugCleaningProcess) cout<<"fHist*Efficiency->Delete()"<<endl;
2274   //MC info, efficiencies
2275         fHistPureEfficiency->Delete();
2276         fHistEfficiency->Delete();
2277   //data for G3/F, continued
2278   if( fWhichParticle == "Lambda" || fWhichParticle == "AntiLambda" ){
2279       h1dPanosCorrections->Delete();
2280       fitGeant3FlukaCorr->Delete();
2281   }
2282
2283   //histograms, feeddown
2284   if (lDebugCleaningProcess) cout<<"f2dFeedDown*->Delete()"<<endl;
2285   if( fFDSwitch != "NoFD" && fWhichParticle != "K0Short"){
2286     //fHistPtXiReference->Delete();
2287     //fHistMCCountbyptXiFeedDown->Delete();
2288     f2dFeedDownMatrix->Delete();
2289     f2dFeedDownEfficiency->Delete();
2290     f2dFeedDownEfficiencyGFCorrected->Delete();
2291     fHistFeeddownSubtraction->Delete();
2292   }
2293   if (lDebugCleaningProcess) cout<<"fHistPt*->Delete()"<<endl;
2294   //Corrected Spectra Histograms
2295   fHistPtLambda->Delete();
2296   fHistPtAntiLambda->Delete();
2297   fHistPtK0Short->Delete();
2298
2299   //Exit Batch Mode
2300   gROOT->SetBatch (kFALSE);
2301
2302   cout<<"--------------------------------------------------------"<<endl;
2303   cout<<" There, done! "<<endl;
2304   cout<<endl;
2305
2306
2307
2308
2309 }