PWGLF/totEt/macros/hadEt/CorrPIDLevyFit.C

   1 //Christine Nattrass, University of Tennessee at Knoxville
   2 //This macro is to calculate the contributions to Et from various particles based on Levy fits to ALICE data at 900 GeV
   3 //It uses d^2N/dpTdy from the papers, weighs it by Et, and integrates over all pT.
   4 //A=0 for mesons
   5 //A=1 for antinucleons since they would anihilate in the calorimeter if they were observed by a calorimeter
   6 //A=-1 for nucleons since they would not anihilate and their mass would not be measured
   7 //At the end this is used to calculate the pid calculation correction
   8 class AliAnalysisLevyPtModified{
   9 public:
  10   virtual ~AliAnalysisLevyPtModified(){;};
  11
  12
  13   Double_t Evaluate(Double_t *pt, Double_t *par)
  14   {
  15     Double_t ldNdy  = par[0];
  16     Double_t l2pi   = 2*TMath::Pi();
  17     Double_t lTemp = par[1];
  18     Double_t lPower = par[2];
  19     Double_t lMass  = par[3];
  20     Double_t A = par[4];
  21     Double_t lMassEt  = par[5];//only used for calculating Et
  22     //this is the Et we would calculate if we had identified the particle as having a mass lMassEt
  23     Double_t Et = TMath::Sqrt(pt[0]*pt[0]+lMassEt*lMassEt)+A*lMassEt;
  24
  25     Double_t lBigCoef = ((lPower-1)*(lPower-2)) / (l2pi*lPower*lTemp*(lPower*lTemp+lMass*(lPower-2)));
  26     Double_t lInPower = 1 + (TMath::Sqrt(pt[0]*pt[0]+lMass*lMass)-lMass) / (lPower*lTemp);
  27     //This is the density of particles times the Et weight
  28     return ldNdy *Et* pt[0] * lBigCoef * TMath::Power(lInPower,(-1)*lPower);
  29   }
  30   ClassDef(AliAnalysisLevyPtModified, 1);
  31 };
  32 class AliAnalysisLevyPtModifiedStrangeness{
  33 public:
  34   virtual ~AliAnalysisLevyPtModifiedStrangeness(){;};
  35
  36
  37   Double_t Evaluate(Double_t *pt, Double_t *par)
  38   {
  39     Double_t ldNdy  = par[0];
  40     Double_t l2pi   = 2*TMath::Pi();
  41     Double_t lTemp = par[1];
  42     Double_t lPower = par[2];
  43     Double_t lMass  = par[3];
  44     Double_t A = par[4];
  45     Double_t lMassEt  = par[5];//only used for calculating Et
  46     //this is the Et we would calculate if we had identified the particle as having a mass lMassEt
  47     Double_t Et = TMath::Sqrt(pt[0]*pt[0]+lMassEt*lMassEt)+A*lMassEt;
  48
  49     Double_t lBigCoef = ((lPower-1)*(lPower-2)) / (l2pi*lPower*lTemp*(lPower*lTemp+lMass*(lPower-2)));
  50     Double_t lInPower = 1 + (TMath::Sqrt(pt[0]*pt[0]+lMass*lMass)-lMass) / (lPower*lTemp);
  51     Double_t strangeness = (0.4333333*pt[0]-0.0666666)/(0.2333333333*pt[0]-0.01666666);
  52     if(strangeness<1.0 || pt[0]<0.1) strangeness = 1.0;
  53     //This is the density of particles times the Et weight
  54     return ldNdy *strangeness* Et* pt[0] * lBigCoef * TMath::Power(lInPower,(-1)*lPower);
  55   }
  56   ClassDef(AliAnalysisLevyPtModifiedStrangeness, 1);
  57 };
  58 class AliAnalysisLevyPtModifiedBaryonEnhanced{
  59 public:
  60   virtual ~AliAnalysisLevyPtModifiedBaryonEnhanced(){;};
  61
  62
  63   Double_t Evaluate(Double_t *pt, Double_t *par)
  64   {
  65     Double_t ldNdy  = par[0];
  66     Double_t l2pi   = 2*TMath::Pi();
  67     Double_t lTemp = par[1];
  68     Double_t lPower = par[2];
  69     Double_t lMass  = par[3];
  70     Double_t A = par[4];
  71     Double_t lMassEt  = par[5];//only used for calculating Et
  72     Double_t lBary0 = par[6];
  73     Double_t lBary1 = par[7];
  74     Double_t lBary2 = par[8];
  75     Double_t lBary3 = par[9];
  76     Double_t lBary4 = par[10];
  77     Double_t lBary5 = par[11];
  78     //this is the Et we would calculate if we had identified the particle as having a mass lMassEt
  79     Double_t Et = TMath::Sqrt(pt[0]*pt[0]+lMassEt*lMassEt)+A*lMassEt;
  80
  81     Double_t lBigCoef = ((lPower-1)*(lPower-2)) / (l2pi*lPower*lTemp*(lPower*lTemp+lMass*(lPower-2)));
  82     Double_t lInPower = 1 + (TMath::Sqrt(pt[0]*pt[0]+lMass*lMass)-lMass) / (lPower*lTemp);
  83     //this is the baryon enhancement factor
  84     Double_t baryonEnhancement = lBary0*pow(pt[0],lBary1)*exp(-pow(pt[0]/lBary2,lBary3))/(lBary4+lBary5*pt[0]);
  85     if(baryonEnhancement<1.0) baryonEnhancement = 1.0;
  86     Double_t strangeness = (0.4333333*pt[0]-0.0666666)/(0.2333333333*pt[0]-0.01666666);
  87     if(strangeness<1.0 || pt[0]<0.1) strangeness = 1.0;
  88     //This is the density of particles times the Et weight
  89     return ldNdy *Et*strangeness* baryonEnhancement*pt[0] * lBigCoef * TMath::Power(lInPower,(-1)*lPower);  }
  90   ClassDef(AliAnalysisLevyPtModifiedBaryonEnhanced, 1);
  91 };
  92 class AliAnalysisLevyPtModifiedStrangenessBaryonEnhanced{
  93 public:
  94   virtual ~AliAnalysisLevyPtModifiedStrangenessBaryonEnhanced(){;};
  95
  96
  97   Double_t Evaluate(Double_t *pt, Double_t *par)
  98   {
  99     Double_t ldNdy  = par[0];
 100     Double_t l2pi   = 2*TMath::Pi();
 101     Double_t lTemp = par[1];
 102     Double_t lPower = par[2];
 103     Double_t lMass  = par[3];
 104     Double_t A = par[4];
 105     Double_t lMassEt  = par[5];//only used for calculating Et
 106     Double_t lBary0 = par[6];
 107     Double_t lBary1 = par[7];
 108     Double_t lBary2 = par[8];
 109     Double_t lBary3 = par[9];
 110     Double_t lBary4 = par[10];
 111     Double_t lBary5 = par[11];
 112     //this is the Et we would calculate if we had identified the particle as having a mass lMassEt
 113     Double_t Et = TMath::Sqrt(pt[0]*pt[0]+lMassEt*lMassEt)+A*lMassEt;
 114
 115     Double_t lBigCoef = ((lPower-1)*(lPower-2)) / (l2pi*lPower*lTemp*(lPower*lTemp+lMass*(lPower-2)));
 116     Double_t lInPower = 1 + (TMath::Sqrt(pt[0]*pt[0]+lMass*lMass)-lMass) / (lPower*lTemp);
 117     //this is the baryon enhancement factor
 118     Double_t baryonEnhancement = lBary0*pow(pt[0],lBary1)*exp(-pow(pt[0]/lBary2,lBary3))/(lBary4+lBary5*pt[0]);
 119     if(baryonEnhancement<1.0) baryonEnhancement = 1.0;
 120     //This is the density of particles times the Et weight
 121     return ldNdy *Et* baryonEnhancement*pt[0] * lBigCoef * TMath::Power(lInPower,(-1)*lPower);  }
 122   ClassDef(AliAnalysisLevyPtModifiedStrangenessBaryonEnhanced, 1);
 123 };
 124
 125 void CorrPIDLevyFit(bool hadronic = false){
 126   //pt cuts where we can identify things
 127   float protoncut = 0.9+0.25;
 128   float protoncuterr = 0.25;
 129   float kaoncut = 0.45+0.25;
 130   float kaoncuterr = 0.25;
 131   //pt cut
 132   float ptcut = 0.1;
 133 // particle          & $\frac{dN}{dy}$    & T (GeV)           & n               & $\frac{dE_T}{dy}$& a &\ET \\ \hline
 134 // $\pi^{+}+\pi^{-}$ & 2.977  $\pm$ 0.15  & 0.126 $\pm$ 0.001 & 7.82 $\pm$ 0.1  &                  &   &     \\
 135   AliAnalysisLevyPtModified *function = new AliAnalysisLevyPtModified();
 136   fPion = new TF1("fPion",function, &AliAnalysisLevyPtModified::Evaluate,0,50,6,"AliAnalysisLevyPtModified","Evaluate");
 137   fPion->SetParameter(0,2.977);//dN/dy
 138   fPion->SetParameter(1,0.126);//T
 139   fPion->SetParameter(2,7.82);//n
 140   fPion->SetParameter(3,0.13957);//mass
 141   fPion->SetParameter(4,0);//A=0 for pions
 142   fPion->SetParameter(5,0.13957);//mass et
 143   float integralPion = fPion->Integral(ptcut,50);
 144   float integralErrPion = integralPion*0.007/2.977;
 145   float myerrorPionT = 0.0;
 146   float myerrorPionn = 0.0;
 147   float tmpPion;
 148   float T = 0.126;
 149   float Terr = 0.001;
 150   float n = 7.82;
 151   float nerr = 0.1;
 152   fPion->SetParameter(1,T+Terr);//T
 153   tmpPion = fPion->Integral(ptcut,50);
 154   myerrorPionT = TMath::Abs(integralPion-tmpPion);
 155   fPion->SetParameter(1,T-Terr);//T
 156   tmpPion = fPion->Integral(ptcut,50);
 157   if(TMath::Abs(integralPion-tmpPion)>myerrorPionT) myerrorPionT = TMath::Abs(integralPion-tmpPion);
 158   fPion->SetParameter(1,T);//T
 159   fPion->SetParameter(2,n+nerr);//n
 160   tmpPion = fPion->Integral(ptcut,50);
 161   myerrorPionn = TMath::Abs(integralPion-tmpPion);
 162   fPion->SetParameter(2,n-nerr);//n
 163   tmpPion = fPion->Integral(ptcut,50);
 164   if(TMath::Abs(integralPion-tmpPion)>myerrorPionn) myerrorPionn = TMath::Abs(integralPion-tmpPion);
 165   //This isn't strictly correct because the errors on the parameters should be correlated but it's close
 166   //To get the correct error one would have to fit the spectra data to get the covariance matrix...
 167   cout<<"Errors: dN/dy "<<integralErrPion<<" T "<<myerrorPionT<<" n "<<myerrorPionn<<endl;
 168   integralErrPion = TMath::Sqrt(TMath::Power(integralErrPion,2)+TMath::Power(myerrorPionT,2)+TMath::Power(myerrorPionn,2));
 169   cout<<"Pion Et = "<<integralPion<<"$\\pm$"<<integralErrPion<<endl;
 170
 171
 172
 173
 174 // particle          & $\frac{dN}{dy}$    & T (GeV)           & n               & $\frac{dE_T}{dy}$& a &\ET \\ \hline
 175 // $K^{+}+K^{-}$     & 0.366 $\pm$ 0.03   & 0.160 $\pm$ 0.006 & 6.087 $\pm$ 0.4 &                  &   &     \\
 176   fKaon = new TF1("fKaon",function, &AliAnalysisLevyPtModified::Evaluate,0,50,6,"AliAnalysisLevyPtModified","Evaluate");
 177   fKaon->SetParameter(0,0.366);//dN/dy
 178   fKaon->SetParameter(1,0.160);//T
 179   fKaon->SetParameter(2,6.087);//n
 180   fKaon->SetParameter(3,0.493677);//mass
 181   fKaon->SetParameter(4,0);//A=0 for kaons
 182   fKaon->SetParameter(5,0.493677);//mass et
 183   fKaonStrange = new TF1("fKaonStrange",function, &AliAnalysisLevyPtModifiedStrangeness::Evaluate,0,50,6,"AliAnalysisLevyPtModifiedStrangeness","Evaluate");
 184   for(int i=0;i<6;i++){fKaonStrange->SetParameter(i,fKaon->GetParameter(i));}
 185   float integralKaon = fKaon->Integral(ptcut,50);
 186   float integralKaonStrange = fKaonStrange->Integral(ptcut,50);
 187   float integralErrKaon = integralKaon*0.006/0.366;
 188   float myerrorKaonT = 0.0;
 189   float myerrorKaonn = 0.0;
 190   float tmpKaon;
 191   float T = 0.160;
 192   float Terr = 0.006;
 193   float n = 6.087;
 194   float nerr = 0.4;
 195   fKaon->SetParameter(1,T+Terr);//T
 196   tmpKaon = fKaon->Integral(ptcut,50);
 197   myerrorKaonT = TMath::Abs(integralKaon-tmpKaon);
 198   fKaon->SetParameter(1,T-Terr);//T
 199   tmpKaon = fKaon->Integral(ptcut,50);
 200   if(TMath::Abs(integralKaon-tmpKaon)>myerrorKaonT) myerrorKaonT = TMath::Abs(integralKaon-tmpKaon);
 201   fKaon->SetParameter(1,T);//T
 202   fKaon->SetParameter(2,n+nerr);//n
 203   tmpKaon = fKaon->Integral(ptcut,50);
 204   myerrorKaonn = TMath::Abs(integralKaon-tmpKaon);
 205   fKaon->SetParameter(2,n-nerr);//n
 206   tmpKaon = fKaon->Integral(ptcut,50);
 207   if(TMath::Abs(integralKaon-tmpKaon)>myerrorKaonn) myerrorKaonn = TMath::Abs(integralKaon-tmpKaon);
 208   //This isn't strictly correct because the errors on the parameters should be correlated but it's close
 209   cout<<"Errors: dN/dy "<<integralErrKaon<<" T "<<myerrorKaonT<<" n "<<myerrorKaonn<<endl;
 210   integralErrKaon = TMath::Sqrt(TMath::Power(integralErrKaon,2)+TMath::Power(myerrorKaonT,2)+TMath::Power(myerrorKaonn,2));
 211   cout<<"Kaon Et = "<<integralKaon<<"$\\pm$"<<integralErrKaon<<endl;
 212   cout<<"Kaon Et Strangeness Enhanced = "<<integralKaonStrange<<endl;
 213   //now calculate the integral for kaons we can identify:
 214   float integralKaonLowPtMean = fKaon->Integral(ptcut,kaoncut);
 215   float integralKaonLowPtPlus = fKaon->Integral(ptcut,kaoncut+kaoncuterr);
 216   float integralKaonLowPtMinus = fKaon->Integral(ptcut,kaoncut-kaoncuterr);
 217   float integralKaonStrangeLowPtMinus = fKaonStrange->Integral(ptcut,kaoncut-kaoncuterr);
 218   //now we switch the kaon mass to the pion mass
 219   fKaon->SetParameter(5,0.13957);//mass et
 220   fKaonStrange->SetParameter(5,0.13957);//mass et
 221   //and do the integrals
 222   float integralKaonHighPtMean = fKaon->Integral(kaoncut,50);
 223   float integralKaonHighPtPlus = fKaon->Integral(kaoncut+kaoncuterr,50);
 224   float integralKaonHighPtMinus = fKaon->Integral(kaoncut-kaoncuterr,50);
 225   float integralKaonStrangeHighPtMinus = fKaonStrange->Integral(kaoncut-kaoncuterr,50);
 226   float integralKaonNoID = fKaon->Integral(ptcut,50);
 227   cout<<"Kaon Et as pion = "<<integralKaonHighPtMean+integralKaonLowPtMean<<endl;
 228   cout<<"Kaon Et Strange as pion = "<<integralKaonStrangeHighPtMinus+integralKaonStrangeLowPtMinus<<endl;
 229   TF1 *fKaonPion = fKaon->Clone("fKaonPion");
 230   //and change it back just to avoid confusion
 231   fKaon->SetParameter(5,0.493677);//mass et
 232   fKaonStrange->SetParameter(5,0.493677);//mass et
 233
 234 // particle          & $\frac{dN}{dy}$    & T (GeV)           & n               & $\frac{dE_T}{dy}$& a &\ET \\ \hline
 235 // $p + \bar{p}$     &  0.157 $\pm$ 0.012 & 0.196 $\pm$ 0.009 & 8.6 $\pm$ 1.1   &                  &   &     \\
 236   fProton = new TF1("fProton",function, &AliAnalysisLevyPtModified::Evaluate,0,50,6,"AliAnalysisLevyPtModified","Evaluate");
 237   fProton->SetParameter(0,0.08);//dN/dy
 238   fProton->SetParameter(1,0.196);//T
 239   fProton->SetParameter(2,8.6);//n
 240   fProton->SetParameter(3,0.938272);//mass
 241   fProton->SetParameter(4,-1);//A=0 for protons
 242   fProton->SetParameter(5,0.938272);//mass et
 243   fProtonEnhanced = new TF1("fProtonEnhanced",function, &AliAnalysisLevyPtModifiedBaryonEnhanced::Evaluate,0,50,12,"AliAnalysisLevyPtModifiedBaryonEnhanced","Evaluate");
 244   for(int i=0;i<6;i++){fProtonEnhanced->SetParameter(i,fProton->GetParameter(i));}//set all of the spectra parameters to their normal values
 245   //and now set the baryon enhancement parameters
 246   fProtonEnhanced->SetParameter(6,0.900878);
 247   fProtonEnhanced->SetParameter(7,1.38882);
 248   fProtonEnhanced->SetParameter(8,2.6361);
 249   fProtonEnhanced->SetParameter(9,1.37751);
 250   fProtonEnhanced->SetParameter(10,0.5);
 251   fProtonEnhanced->SetParameter(11,-0.03);
 252   float integralProton = fProton->Integral(ptcut,50);
 253   float integralErrProton = integralProton*0.002/0.08;
 254   float myerrorProtonT = 0.0;
 255   float myerrorProtonn = 0.0;
 256   float tmpProton;
 257   float T = 0.196;
 258   float Terr = 0.009;
 259   float n = 8.6;
 260   float nerr = 1.1;
 261   fProton->SetParameter(1,T+Terr);//T
 262   tmpProton = fProton->Integral(ptcut,50);
 263   myerrorProtonT = TMath::Abs(integralProton-tmpProton);
 264   fProton->SetParameter(1,T-Terr);//T
 265   tmpProton = fProton->Integral(ptcut,50);
 266   if(TMath::Abs(integralProton-tmpProton)>myerrorProtonT) myerrorProtonT = TMath::Abs(integralProton-tmpProton);
 267   fProton->SetParameter(1,T);//T
 268   fProton->SetParameter(2,n+nerr);//n
 269   tmpProton = fProton->Integral(ptcut,50);
 270   myerrorProtonn = TMath::Abs(integralProton-tmpProton);
 271   fProton->SetParameter(2,n-nerr);//n
 272   tmpProton = fProton->Integral(ptcut,50);
 273   if(TMath::Abs(integralProton-tmpProton)>myerrorProtonn) myerrorProtonn = TMath::Abs(integralProton-tmpProton);
 274   //This isn't strictly correct because the errors on the parameters should be correlated but it's close
 275   cout<<"Errors: dN/dy "<<integralErrProton<<" T "<<myerrorProtonT<<" n "<<myerrorProtonn<<endl;
 276   integralErrProton = TMath::Sqrt(TMath::Power(integralErrProton,2)+TMath::Power(myerrorProtonT,2)+TMath::Power(myerrorProtonn,2));
 277   cout<<"Proton Et = "<<integralProton<<"$\\pm$"<<integralErrProton<<endl;
 278   //now calculate the integral for kaons we can identify:
 279   float integralProtonLowPtMean = fProton->Integral(ptcut,protoncut);
 280   float integralProtonLowPtPlus = fProton->Integral(ptcut,protoncut+protoncuterr);
 281   float integralProtonLowPtMinus = fProton->Integral(ptcut,protoncut-protoncuterr);
 282   float integralProtonEnhancedTrue =integralProtonLowPtMinus+ fProtonEnhanced->Integral(protoncut-protoncuterr,50);
 283   fProtonEnhanced->SetParameter(6,1.2*fProtonEnhanced->GetParameter(6));
 284   float integralProtonReallyEnhancedTrue = integralProtonLowPtMinus+ fProtonEnhanced->Integral(protoncut-protoncuterr,50);
 285   fProtonEnhanced->SetParameter(6,1.0/1.2*fProtonEnhanced->GetParameter(6));
 286   //now we switch the kaon mass to the pion mass
 287   fProton->SetParameter(5,0.13957);//mass et
 288   fProton->SetParameter(4,0);//A=0 for pions
 289   //and do the integrals
 290   float integralProtonHighPtMean = fProton->Integral(protoncut,50);
 291   float integralProtonHighPtPlus = fProton->Integral(protoncut+protoncuterr,50);
 292   float integralProtonHighPtMinus = fProton->Integral(protoncut-protoncuterr,50);
 293
 294   //now we switch the kaon mass to the pion mass
 295   fProtonEnhanced->SetParameter(5,0.13957);//mass et
 296   fProtonEnhanced->SetParameter(4,0);//A=0 for pions
 297   float integralProtonHighPtMinusEnhanced = fProtonEnhanced->Integral(protoncut-protoncuterr,50);
 298   fProtonEnhanced->SetParameter(6,1.2*fProtonEnhanced->GetParameter(6));
 299   float integralProtonHighPtMinusReallyEnhanced = fProtonEnhanced->Integral(protoncut-protoncuterr,50);
 300   fProtonEnhanced->SetParameter(6,1.0/1.2*fProtonEnhanced->GetParameter(6));
 301
 302   float integralProtonNoID = fProton->Integral(ptcut,50);
 303   float integralProtonEnhancedNoID = integralProtonHighPtMinusEnhanced+integralProtonLowPtMinus;
 304   float integralProtonReallyEnhancedNoID = integralProtonHighPtMinusReallyEnhanced+integralProtonLowPtMinus;
 305   cout<<"Proton Et as pion = "<<integralProtonHighPtMean+integralProtonLowPtMean<<endl;
 306   cout<<"Proton Et really enhanced = "<<integralProtonReallyEnhancedTrue<<"("<<integralProtonHighPtMinusReallyEnhanced+integralProtonLowPtMinus<<")"<<endl;
 307   TF1 *fProtonPion = fProton->Clone("fProtonPion");
 308   //and change it back just to avoid confusion
 309   fProton->SetParameter(5,0.938272);//mass et
 310   fProton->SetParameter(4,-1);//A=0 for protons
 311   fProtonEnhanced->SetParameter(5,0.938272);//mass et
 312   fProtonEnhanced->SetParameter(4,-1);//A=0 for protons
 313
 314
 315
 316   //Antiprotons...
 317   fProton->SetParameter(0,0.077);//dN/dy
 318   fProton->SetParameter(2,n);//n
 319   fProton->SetParameter(4,1);//A=0 for protons
 320   fAntiProtonEnhanced = new TF1("fAntiProtonEnhanced",function, &AliAnalysisLevyPtModifiedBaryonEnhanced::Evaluate,0,50,12,"AliAnalysisLevyPtModifiedBaryonEnhanced","Evaluate");
 321   for(int i=0;i<6;i++){fAntiProtonEnhanced->SetParameter(i,fProton->GetParameter(i));}//set all of the spectra parameters to their normal values
 322   //and now set the baryon enhancement parameters
 323   fAntiProtonEnhanced->SetParameter(6,0.900878);
 324   fAntiProtonEnhanced->SetParameter(7,1.38882);
 325   fAntiProtonEnhanced->SetParameter(8,2.6361);
 326   fAntiProtonEnhanced->SetParameter(9,1.37751);
 327   fAntiProtonEnhanced->SetParameter(10,0.5);
 328   fAntiProtonEnhanced->SetParameter(11,-0.03);
 329   float integralAntiProton = fProton->Integral(ptcut,50);
 330   float integralErrAntiProton = integralAntiProton*0.002/0.077;
 331   float myerrorAntiProtonT = 0.0;
 332   float myerrorAntiProtonn = 0.0;
 333   float tmpAntiProton;
 334   float T = 0.196;
 335   float Terr = 0.009;
 336   float n = 8.6;
 337   float nerr = 1.1;
 338   fProton->SetParameter(1,T+Terr);//T
 339   tmpAntiProton = fProton->Integral(ptcut,50);
 340   myerrorAntiProtonT = TMath::Abs(integralAntiProton-tmpAntiProton);
 341   fProton->SetParameter(1,T-Terr);//T
 342   tmpAntiProton = fProton->Integral(ptcut,50);
 343   if(TMath::Abs(integralAntiProton-tmpAntiProton)>myerrorAntiProtonT) myerrorAntiProtonT = TMath::Abs(integralAntiProton-tmpAntiProton);
 344   fProton->SetParameter(1,T);//T
 345   fProton->SetParameter(2,n+nerr);//n
 346   tmpAntiProton = fProton->Integral(ptcut,50);
 347   myerrorAntiProtonn = TMath::Abs(integralAntiProton-tmpAntiProton);
 348   fProton->SetParameter(2,n-nerr);//n
 349   tmpAntiProton = fProton->Integral(ptcut,50);
 350   if(TMath::Abs(integralAntiProton-tmpAntiProton)>myerrorAntiProtonn) myerrorAntiProtonn = TMath::Abs(integralAntiProton-tmpAntiProton);
 351   //This isn't strictly correct because the errors on the parameters should be correlated but it's close
 352   cout<<"Errors: dN/dy "<<integralErrAntiProton<<" T "<<myerrorAntiProtonT<<" n "<<myerrorAntiProtonn<<endl;
 353   integralErrAntiProton = TMath::Sqrt(TMath::Power(integralErrAntiProton,2)+TMath::Power(myerrorAntiProtonT,2)+TMath::Power(myerrorAntiProtonn,2));
 354   cout<<"AntiProton Et = "<<integralAntiProton<<"$\\pm$"<<integralErrAntiProton<<endl;
 355   TF1 *fAntiProton = fProton->Clone("fAntiProton");
 356   //now calculate the integral for kaons we can identify:
 357   float integralAntiProtonLowPtMean = fProton->Integral(ptcut,protoncut);
 358   float integralAntiProtonLowPtPlus = fProton->Integral(ptcut,protoncut+protoncuterr);
 359   float integralAntiProtonLowPtMinus = fProton->Integral(ptcut,protoncut-protoncuterr);
 360   float integralAntiProtonEnhancedTrue =integralAntiProtonLowPtMinus+ fAntiProtonEnhanced->Integral(protoncut-protoncuterr,50);
 361   fAntiProtonEnhanced->SetParameter(6,1.2*fAntiProtonEnhanced->GetParameter(6));
 362   float integralAntiProtonReallyEnhancedTrue = integralAntiProtonLowPtMinus+ fAntiProtonEnhanced->Integral(protoncut-protoncuterr,50);
 363   fAntiProtonEnhanced->SetParameter(6,1.0/1.2*fAntiProtonEnhanced->GetParameter(6));
 364   //now we switch the kaon mass to the pion mass
 365   fProton->SetParameter(5,0.13957);//mass et
 366   fProton->SetParameter(4,0);//A=0 for pions
 367   //and do the integrals
 368   float integralAntiProtonHighPtMean = fProton->Integral(protoncut,50);
 369   float integralAntiProtonHighPtPlus = fProton->Integral(protoncut+protoncuterr,50);
 370   float integralAntiProtonHighPtMinus = fProton->Integral(protoncut-protoncuterr,50);
 371
 372   //now we switch the kaon mass to the pion mass
 373   fAntiProtonEnhanced->SetParameter(5,0.13957);//mass et
 374   fAntiProtonEnhanced->SetParameter(4,0);//A=0 for pions
 375   float integralAntiProtonHighPtMinusEnhanced = fAntiProtonEnhanced->Integral(protoncut-protoncuterr,50);
 376   fAntiProtonEnhanced->SetParameter(6,1.2*fAntiProtonEnhanced->GetParameter(6));
 377   float integralAntiProtonHighPtMinusReallyEnhanced = fAntiProtonEnhanced->Integral(protoncut-protoncuterr,50);
 378   fAntiProtonEnhanced->SetParameter(6,1.0/1.2*fAntiProtonEnhanced->GetParameter(6));
 379
 380   float integralAntiProtonNoID = fProton->Integral(ptcut,50);
 381   cout<<"AntiProton Et as pion = "<<integralAntiProtonHighPtMean+integralAntiProtonLowPtMean<<endl;
 382   cout<<"AntiProton Et really enhanced = "<<integralAntiProtonReallyEnhancedTrue<<"("<<integralAntiProtonHighPtMinusReallyEnhanced+integralAntiProtonLowPtMinus<<")"<<endl;
 383   TF1 *fAntiProtonPion = fProton->Clone("fAntiProtonPion");
 384   //and change it back just to avoid confusion
 385   fProton->SetParameter(5,0.938272);//mass
 386   //cout<<"CHECK "<<fProton->GetParameter(5)<<endl;
 387   fProton->SetParameter(4,-1);//A=0 for protons
 388   fAntiProtonEnhanced->SetParameter(5,0.938272);//mass et
 389   fAntiProtonEnhanced->SetParameter(4,-1);//A=0 for protons
 390
 391   float totTrue = integralPion + integralKaon + integralProton + integralAntiProton;
 392   float totTrueErr = TMath::Sqrt(TMath::Power(integralErrAntiProton,2)+TMath::Power(integralErrProton,2)+TMath::Power(integralErrKaon,2)+TMath::Power(integralErrPion,2));
 393   float totTrueEnhanced = integralPion + integralKaon + integralProtonEnhancedTrue + integralAntiProtonEnhancedTrue;
 394   float totTrueReallyEnhanced = integralPion + integralKaon + integralProtonReallyEnhancedTrue + integralAntiProtonReallyEnhancedTrue;
 395   float totTrueStrange = integralPion + integralKaonStrange + integralProton + integralAntiProton;
 396   float totTrueStrangeEnhanced = integralPion + integralKaonStrange + integralProtonReallyEnhancedTrue + integralAntiProtonReallyEnhancedTrue;
 397   cout<<"totEt "<<totTrue<<"+/-"<<totTrueErr<<endl;
 398
 399   float measEt = integralPion+integralKaonLowPtMean+integralKaonHighPtMean+integralProtonLowPtMean+integralProtonHighPtMean+integralAntiProtonLowPtMean+integralAntiProtonHighPtMean;
 400   float measEtPlus = integralPion+integralKaonLowPtPlus+integralKaonHighPtPlus+integralProtonLowPtPlus+integralProtonHighPtPlus+integralAntiProtonLowPtPlus+integralAntiProtonHighPtPlus;
 401   float measEtMinus = integralPion+integralKaonLowPtMinus+integralKaonHighPtMinus+integralProtonLowPtMinus+integralProtonHighPtMinus+integralAntiProtonLowPtMinus+integralAntiProtonHighPtMinus;
 402   float measEtMinusStrange = integralPion+integralKaonStrangeLowPtMinus+integralKaonStrangeHighPtMinus+integralProtonLowPtMinus+integralProtonHighPtMinus+integralAntiProtonLowPtMinus+integralAntiProtonHighPtMinus;
 403   float measEtEnhanced = integralPion+integralKaonLowPtMinus+integralKaonHighPtMinus+integralProtonLowPtMinus+integralProtonHighPtMinusEnhanced+integralAntiProtonLowPtMinus+integralAntiProtonHighPtMinusEnhanced;
 404   float measEtReallyEnhanced = integralPion+integralKaonLowPtMinus+integralKaonHighPtMinus+integralProtonLowPtMinus+integralProtonHighPtMinusReallyEnhanced+integralAntiProtonLowPtMinus+integralAntiProtonHighPtMinusReallyEnhanced;
 405   cout<<"measEt "<<measEt<<" "<<measEtPlus<<" "<<measEtMinus<<endl;
 406   float measEtStrangeEnhanced = integralPion+integralKaonStrangeLowPtMinus+integralKaonStrangeHighPtMinus+integralProtonLowPtMinus+integralProtonHighPtMinusReallyEnhanced+integralAntiProtonLowPtMinus+integralAntiProtonHighPtMinusReallyEnhanced;
 407   cout<<"measEt "<<measEt<<" "<<measEtPlus<<" "<<measEtMinus<<endl;
 408   float measNoID = integralPion+integralKaonNoID+integralProtonNoID+integralAntiProtonNoID;
 409
 410   float fpid = measEt/totTrue;
 411   float fpiderr = fpid*totTrueErr/totTrue;
 412   float fpidsyserr = TMath::Abs(measEtMinus-measEt)/totTrue;
 413   float fnopid = measNoID/totTrue;
 414   float fnopiderr = fnopid*totTrueErr/totTrue;
 415   if(TMath::Abs(measEtPlus-measEt)> fpidsyserr) fpidsyserr = TMath::Abs(measEtPlus-measEt);
 416   cout<<"fpid "<<fpid<<" +/- "<<fpiderr<<" +/- "<<fpidsyserr<<endl;
 417   cout<<"fpid (no pt uncertainty) "<<measEtMinus/totTrue<<" +/- "<<(measEtMinus/totTrue)*totTrueErr/totTrue<<endl;
 418   cout<<"For no id fpid "<<fnopid<<" +/- "<<fnopiderr<<endl;
 419   cout<<"fpid enhanced "<<measEtEnhanced/totTrueEnhanced<<endl;
 420   cout<<"fpid strange "<<measEtMinusStrange/totTrueStrange<<endl;
 421   cout<<"fpid strange enhanced "<<measEtStrangeEnhanced/totTrueStrangeEnhanced<<endl;
 422   cout<<"fpid really enhanced "<<measEtReallyEnhanced/totTrueReallyEnhanced<<endl;
 423
 424
 425
 426   gStyle->SetOptTitle(0);
 427   gStyle->SetOptStat(0);
 428   gStyle->SetOptFit(0);
 429   TCanvas *c1 = new TCanvas("c1","c1",500,400);
 430   c1->SetTopMargin(0.02);
 431   c1->SetRightMargin(0.02);
 432   c1->SetBorderSize(0);
 433   c1->SetFillColor(0);
 434   c1->SetFillColor(0);
 435   c1->SetBorderMode(0);
 436   c1->SetFrameFillColor(0);
 437   c1->SetFrameBorderMode(0);
 438   c1->SetLogy();
 439   fPion->SetLineColor(1);
 440   fKaon->SetLineColor(2);
 441   fKaonPion->SetLineColor(2);
 442   fKaonPion->SetLineStyle(2);
 443   fProton->SetLineColor(4);
 444   fProtonPion->SetLineColor(4);
 445   fProtonPion->SetLineStyle(2);
 446   fAntiProton->SetLineColor(TColor::kCyan);
 447   fAntiProtonPion->SetLineColor(TColor::kCyan);
 448   fAntiProtonPion->SetLineStyle(2);
 449 //   cout<<"Kaon pion mass "<<Form("%2.4f",fKaonPion->GetParameter(5))<<" A "<<fKaonPion->GetParameter(4)<<endl;
 450 //   cout<<"Proton pion mass "<<fProtonPion->GetParameter(5)<<" A "<<fProtonPion->GetParameter(4)<<endl;
 451 //   cout<<"AntiProton pion mass "<<fAntiProtonPion->GetParameter(5)<<" A "<<fAntiProtonPion->GetParameter(4)<<endl;
 452 //   cout<<"Kaon mass "<<Form("%2.4f",fKaon->GetParameter(5))<<" A "<<fKaon->GetParameter(4)<<endl;
 453 //   cout<<"Proton mass "<<fProton->GetParameter(5)<<" A "<<fProton->GetParameter(4)<<endl;
 454 //   cout<<"AntiProton mass "<<fAntiProton->GetParameter(5)<<" A "<<fAntiProton->GetParameter(4)<<endl;
 455   TH1F *frame = new TH1F("frame","frame",1,0,2);
 456   frame->GetYaxis()->SetTitle("1/N_{eve}d^{2}NE_{T}/dydp_{T}");
 457   frame->GetXaxis()->SetTitle("p_{T}");
 458   //fPion->SetRange(0,2);
 459   frame->SetMinimum(1e-4);
 460   frame->SetMaximum(0.5);
 461   frame->Draw();
 462   fPion->Draw("same");
 463   fKaon->Draw("same");
 464   fProton->Draw("same");
 465   fAntiProton->Draw("same");
 466   fKaonPion->Draw("same");
 467   fProtonPion->Draw("same");
 468   fAntiProtonPion->Draw("same");
 469   TLegend *leg = new  TLegend(0.782258,0.774194,0.925403,0.962366);
 470   leg->AddEntry(fPion,"#pi");
 471   leg->AddEntry(fProton,"p");
 472   leg->AddEntry(fKaon,"K");
 473   leg->AddEntry(fAntiProton,"#bar{p}");
 474   leg->SetFillStyle(0);
 475   leg->SetFillColor(0);
 476   leg->SetBorderSize(0);
 477   leg->SetTextSize(0.06);
 478  leg->Draw();
 479  TLine *lineProton = new TLine(protoncut-protoncuterr,frame->GetMinimum(),protoncut-protoncuterr,frame->GetMaximum());
 480  lineProton->SetLineColor(fProton->GetLineColor());
 481  lineProton->SetLineStyle(3);
 482  lineProton->SetLineWidth(2);
 483  lineProton->Draw();
 484  TLine *lineKaon = new TLine(kaoncut-kaoncuterr,frame->GetMinimum(),kaoncut-kaoncuterr,frame->GetMaximum());
 485  lineKaon->SetLineColor(fKaon->GetLineColor());
 486  lineKaon->SetLineStyle(3);
 487  lineKaon->SetLineWidth(2);
 488  lineKaon->Draw();
 489
 490
 491   TCanvas *c2 = new TCanvas("c2","c2",500,400);
 492   c2->SetTopMargin(0.02);
 493   c2->SetRightMargin(0.02);
 494   c2->SetBorderSize(0);
 495   c2->SetFillColor(0);
 496   c2->SetFillColor(0);
 497   c2->SetBorderMode(0);
 498   c2->SetFrameFillColor(0);
 499   c2->SetFrameBorderMode(0);
 500   c2->SetLogy();
 501   fPionSpectrum = (TF1*) fPion->Clone("fPionSpectrum");
 502   fKaonSpectrum = (TF1*) fKaon->Clone("fKaonSpectrum");
 503   fProtonSpectrum = (TF1*) fProton->Clone("fProtonSpectrum");
 504   fAntiProtonSpectrum = (TF1*) fAntiProton->Clone("fAntiProtonSpectrum");
 505   fPionSpectrum->SetLineColor(1);
 506   fKaonSpectrum->SetLineColor(2);
 507   fProtonSpectrum->SetLineColor(4);
 508   fAntiProtonSpectrum->SetLineColor(TColor::kCyan);
 509   fKaonSpectrum->SetParameter(5,0.0);//mass et
 510   fProtonSpectrum->SetParameter(5,0.0);//mass et
 511   fPionSpectrum->SetParameter(5,0.0);//mass et
 512   fAntiProtonSpectrum->SetParameter(5,0.0);//mass et
 513   TH1F *frame2 = new TH1F("frame2","frame2",1,0,2);
 514   frame2->GetYaxis()->SetTitle("1/N_{eve}d^{2}N/dydp_{T}");
 515   frame2->GetXaxis()->SetTitle("p_{T}");
 516   //fPionSpectrum->SetRange(0,2);
 517   frame2->SetMinimum(1e-5);
 518   frame2->SetMaximum(1);
 519   frame2->Draw();
 520   fPionSpectrum->Draw("same");
 521   fKaonSpectrum->Draw("same");
 522   fProtonSpectrum->Draw("same");
 523   //fAntiProtonSpectrum->Draw("same");
 524   TLegend *leg2 = new  TLegend(0.782258,0.774194,0.925403,0.962366);
 525   leg2->AddEntry(fPionSpectrum,"#pi");
 526   leg2->AddEntry(fProtonSpectrum,"p,#bar{p}");
 527   leg2->AddEntry(fKaonSpectrum,"K");
 528   leg2->SetFillStyle(0);
 529   leg2->SetFillColor(0);
 530   leg2->SetBorderSize(0);
 531   leg2->SetTextSize(0.06);
 532  leg2->Draw();
 533
 534   TCanvas *c3 = new TCanvas("c3","c3",500,400);
 535   c3->SetTopMargin(0.02);
 536   c3->SetRightMargin(0.02);
 537   c3->SetBorderSize(0);
 538   c3->SetFillColor(0);
 539   c3->SetFillColor(0);
 540   c3->SetBorderMode(0);
 541   c3->SetFrameFillColor(0);
 542   c3->SetFrameBorderMode(0);
 543   //c3->SetLogy();
 544   TH1F *frame3 = new TH1F("frame3","frame3",1,0,2);
 545   frame3->GetYaxis()->SetTitle("E_{T}");
 546   frame3->GetXaxis()->SetTitle("p_{T}");
 547   //fPionSpectrum->SetRange(0,2);
 548   frame3->SetMinimum(0.0);
 549   frame3->SetMaximum(4);
 550   frame3->Draw();
 551
 552   TF1 *fPionEt = new TF1("fPionEt","sqrt(x*x+[0]*[0])",0,2);
 553   fPionEt->SetParameter(0,0.13957);
 554   fPionEt->Draw("same");
 555   TF1 *fKaonEt = new TF1("fKaonEt","sqrt(x*x+[0]*[0])",0,2);
 556   fKaonEt->SetParameter(0,0.493677);
 557   fKaonEt->Draw("same");
 558   TF1 *fProtonEt = new TF1("fProtonEt","sqrt(x*x+[0]*[0])-[0]",0,2);
 559   fProtonEt->SetParameter(0,0.938272);
 560   fProtonEt->Draw("same");
 561   TF1 *fAntiProtonEt = new TF1("fAntiProtonEt","sqrt(x*x+[0]*[0])+[0]",0,2);
 562   fAntiProtonEt->SetParameter(0,0.938272);
 563   fAntiProtonEt->Draw("same");
 564   fPionEt->SetLineColor(1);
 565   fKaonEt->SetLineColor(2);
 566   fProtonEt->SetLineColor(4);
 567   fAntiProtonEt->SetLineColor(TColor::kCyan);
 568   TLegend *leg3 = new  TLegend(0.782258,0.774194,0.925403,0.962366);
 569   leg3->AddEntry(fPionEt,"#pi");
 570   leg3->AddEntry(fProtonEt,"p");
 571   leg3->AddEntry(fAntiProtonEt,"#bar{p}");
 572   leg3->AddEntry(fKaonEt,"K");
 573   leg3->SetFillStyle(0);
 574   leg3->SetFillColor(0);
 575   leg3->SetBorderSize(0);
 576   leg3->SetTextSize(0.06);
 577   leg3->Draw();
 578
 579   c1->SaveAs("pics/PIDEtweight.eps");
 580   c2->SaveAs("pics/PIDSpectra.eps");
 581   c3->SaveAs("pics/PIDEt.eps");
 582
 583
 584 }
 585