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