[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / charmFlow / AliEventPlaneResolutionHandler.cxx

/**************************************************************************
 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

#include <TMath.h>
#include <TH1F.h>
#include <TFile.h>
#include "AliVertexingHFUtils.h"
#include "AliEventPlaneResolutionHandler.h"
#include "AliLog.h"

/* $Id$ */

///////////////////////////////////////////////////////////////////
//                                                               //
// Class with functions useful for different D2H analyses        //
// - event plane resolution                                      //
// - <pt> calculation with side band subtraction                 //
// - tracklet multiplicity calculation                            //
// Origin: F.Prino, Torino, prino@to.infn.it                     //
//                                                               //
///////////////////////////////////////////////////////////////////

ClassImp(AliEventPlaneResolutionHandler)

//______________________________________________________________________
AliEventPlaneResolutionHandler::AliEventPlaneResolutionHandler():TObject(),
  fEventPlane(kTPCPosEta),
  fResolOption(kTwoEtaSub),
  fMinCent(30.),
  fMaxCent(50.),
  fCorrHistoName1(""),
  fCorrHistoName2(""),
  fCorrHistoName3(""),
  fNsubevents(2),
  fExtrapToFull(kFALSE),
  fWeight(0),
  fRootFileName("$ALICE_ROOT/PWGHF/vertexingHF/charmFlow/EventPlaneResolutionHistos.root")
{
  // Default contructor
  InitializeWeights();
}
//______________________________________________________________________
AliEventPlaneResolutionHandler::AliEventPlaneResolutionHandler(TString filename):TObject(),
  fEventPlane(kTPCPosEta),
  fResolOption(kTwoEtaSub),
  fMinCent(30.),
  fMaxCent(50.),
  fCorrHistoName1(""),
  fCorrHistoName2(""),
  fCorrHistoName3(""),
  fNsubevents(2),
  fExtrapToFull(kFALSE),
  fWeight(0),
  fRootFileName(filename.Data())
{
  // Standard contructor
  InitializeWeights();
}
//______________________________________________________________________
Double_t AliEventPlaneResolutionHandler::GetEventPlaneResolution(Double_t minCent, Double_t maxCent){
  // method to compute the event plane resolution starting from correlations

  TFile* inputFile=new TFile(fRootFileName.Data());
  if(!inputFile || (inputFile && !inputFile->IsOpen())){
    printf("Root file not opened --- Return dummy resolution value\n");
    return -1.;
  }
  Bool_t isOk=SetHistoNames();
  if(!isOk){
    printf("Return dummy resolution value\n");
    return -1.;
  }
  
  Int_t minCentrTimesTen=(Int_t)(minCent*10);
  Int_t maxCentrTimesTen=(Int_t)(maxCent*10);
  Int_t ncBin=minCentrTimesTen/25;
  Int_t whichWei=fWeight;
  if(minCent>50.5 || maxCent>50.5) fWeight=0;
  Double_t weight=-1.;  
  if(fWeight==1) weight=fNcoll[ncBin];
  else if(fWeight==2 && minCent>29.5 && maxCent<50.5)  weight=fYield24[ncBin];
  else if(fWeight==3 && minCent>29.5 && maxCent<50.5)  weight=fYield46[ncBin];
  else if(fWeight==4 && minCent>29.5 && maxCent<50.5)  weight=fYield612[ncBin];
  TH1F* hevpls2=0x0;
  TH1F* hevpls3=0x0;
  TH1F* hevpls1=(TH1F*)inputFile->Get(Form("%sCentr%d_%d",fCorrHistoName1.Data(),minCentrTimesTen,minCentrTimesTen+25));
  hevpls1->SetName(Form("%sCentr%d_%d",fCorrHistoName1.Data(),minCentrTimesTen,maxCentrTimesTen));
  if(weight>0.) hevpls1->Scale(weight);
  if(fNsubevents==3){
    hevpls2=(TH1F*)inputFile->Get(Form("%sCentr%d_%d",fCorrHistoName2.Data(),minCentrTimesTen,minCentrTimesTen+25));
    hevpls2->SetName(Form("%sCentr%d_%d",fCorrHistoName2.Data(),minCentrTimesTen,maxCentrTimesTen));
    hevpls3=(TH1F*)inputFile->Get(Form("%sCentr%d_%d",fCorrHistoName3.Data(),minCentrTimesTen,minCentrTimesTen+25));
    hevpls3->SetName(Form("%sCentr%d_%d",fCorrHistoName2.Data(),minCentrTimesTen,maxCentrTimesTen));
    if(weight>0.){
      hevpls2->Scale(weight);
      hevpls3->Scale(weight);
    }
  }
  for(Int_t iCentr=minCentrTimesTen+25; iCentr<maxCentrTimesTen; iCentr+=25){
    ncBin=iCentr/25;
    weight=-1.;  
    if(fWeight==1) weight=fNcoll[ncBin];
    else if(fWeight==2 && minCent>29.5 && maxCent<50.5)  weight=fYield24[ncBin];
    else if(fWeight==3 && minCent>29.5 && maxCent<50.5)  weight=fYield46[ncBin];
    else if(fWeight==4 && minCent>29.5 && maxCent<50.5)  weight=fYield612[ncBin];
    TH1F* htmp1=(TH1F*)inputFile->Get(Form("%sCentr%d_%d",fCorrHistoName1.Data(),iCentr,iCentr+25));
    if(weight>0.) htmp1->Scale(weight);
    hevpls1->Add(htmp1);
    if(fNsubevents==3){
      TH1F* htmp2=(TH1F*)inputFile->Get(Form("%sCentr%d_%d",fCorrHistoName2.Data(),iCentr,iCentr+25));
      TH1F* htmp3=(TH1F*)inputFile->Get(Form("%sCentr%d_%d",fCorrHistoName3.Data(),iCentr,iCentr+25));
      if(weight>0.){
	htmp2->Scale(weight);
	htmp3->Scale(weight);
      }
      hevpls2->Add(htmp2);
      hevpls3->Add(htmp3);
    }
  }
  fWeight=whichWei;
  Double_t resol=1.;
  if(fNsubevents==2){
    if(fExtrapToFull) resol=AliVertexingHFUtils::GetFullEvResol(hevpls1);
    else resol=AliVertexingHFUtils::GetSubEvResol(hevpls1);
  }else if(fNsubevents==3){
    Double_t correl1=hevpls1->GetMean();
    Double_t correl2=hevpls2->GetMean();
    Double_t correl3=hevpls3->GetMean();
    resol=TMath::Sqrt(correl1*correl2/correl3);
  }
  inputFile->Close();
  return resol;
}
//______________________________________________________________________
Bool_t AliEventPlaneResolutionHandler::SetHistoNames(){
  // set names for histograms to be read
  if(fEventPlane==kTPCPosEta && fResolOption==kTwoRandSub){
    fCorrHistoName1="hCorrelRandSubTPCPosEta";
    fNsubevents=2;
    fExtrapToFull=kTRUE;
  }else if(fEventPlane==kTPCPosEta && fResolOption==kTwoChargeSub){
    fCorrHistoName1="hCorrelChargeSubTPCPosEta";
    fNsubevents=2;
    fExtrapToFull=kTRUE;
  }else if(fEventPlane==kTPCPosEta && fResolOption==kTwoEtaSub){
    fCorrHistoName1="hCorrelTPCEtaPosTPCEtaNeg";
    fNsubevents=2;
    fExtrapToFull=kFALSE;
  }else if(fEventPlane==kTPCPosEta && fResolOption==kThreeSub){
    fCorrHistoName1="hCorrelTPCPosEtaV0A";
    fCorrHistoName2="hCorrelTPCPosEtaV0C";
    fCorrHistoName3="hCorrelV0AV0C";
    fNsubevents=3;
    fExtrapToFull=kFALSE;
  }else if(fEventPlane==kTPCFullEta && fResolOption==kTwoRandSub){
    fCorrHistoName1="hCorrelRandSubTPCFullEta";
    fNsubevents=2;
    fExtrapToFull=kTRUE;
  }else if(fEventPlane==kTPCFullEta && fResolOption==kTwoEtaSub){
    fCorrHistoName1="hCorrelTPCEtaPosTPCEtaNeg";
    fNsubevents=2;
    fExtrapToFull=kTRUE;
  }else if(fEventPlane==kTPCFullEta && fResolOption==kThreeSub){
    fCorrHistoName1="hCorrelTPCFullEtaV0A";
    fCorrHistoName2="hCorrelTPCFullEtaV0C";
    fCorrHistoName3="hCorrelV0AV0C";
    fNsubevents=3;
    fExtrapToFull=kFALSE;
  }else if(fEventPlane==kVZERO && fResolOption==kThreeSub){
    fCorrHistoName1="hCorrelTPCEtaPosV0";
    fCorrHistoName2="hCorrelTPCEtaNegV0";
    fCorrHistoName3="hCorrelTPCEtaPosTPCEtaNeg";
    fNsubevents=3;
    fExtrapToFull=kFALSE;
  }else if(fEventPlane==kVZERO && fResolOption==kThreeSubTPCGap){
    fCorrHistoName1="hCorrelTPCEtaGt02V0";
    fCorrHistoName2="hCorrelTPCEtaLt02V0";
    fCorrHistoName3="hCorrelTPCEtaGt02TPCEtaLt02";
    fNsubevents=3;
    fExtrapToFull=kFALSE;
  }else if(fEventPlane==kVZEROA && fResolOption==kThreeSub){
    fCorrHistoName1="hCorrelTPCFullEtaV0A";
    fCorrHistoName2="hCorrelV0AV0C";
    fCorrHistoName3="hCorrelTPCFullEtaV0C";
    fNsubevents=3;
    fExtrapToFull=kFALSE;
  }else if(fEventPlane==kVZEROC && fResolOption==kThreeSub){
    fCorrHistoName1="hCorrelTPCFullEtaV0C";
    fCorrHistoName2="hCorrelV0AV0C";
    fCorrHistoName3="hCorrelTPCFullEtaV0A";
    fNsubevents=3;
    fExtrapToFull=kFALSE;
  }else{
    printf("Not possible to compute the resolution with the required option\n");
    return kFALSE;
  }
  return kTRUE;
}
//______________________________________________________________________
void AliEventPlaneResolutionHandler::InitializeWeights(){
  // Ncoll values in 2.5% wide classes
  Double_t ncoll[20]={1790.77,1578.44,1394.82,1236.17,1095.08,
		      969.86,859.571,759.959,669.648,589.588,
		      516.039,451.409,392.853,340.493,294.426,
		      252.385,215.484,183.284,155.101,130.963};
  Double_t y24[20]={1,1,1,1,
		    1,1,1,1,1,1,1,1,
		    446,402,339,277,218,185,184,104};

  Double_t y46[20]={1,1,1,1,
		    1,1,1,1,1,1,1,1,
		    49,150,127,105,88,71,71,53};
  Double_t y612[20]={1,1,1,1,
		    1,1,1,1,1,1,1,1,
		    89,92,79,66,57,59,38,58};
  for(Int_t i=0; i<20; i++){
    fNcoll[i]=ncoll[i];
    fYield24[i]=y24[i];
    fYield46[i]=y46[i];
    fYield612[i]=y612[i];
  }
}
Commit	Line	Data
f5e6ce8f	1	/**************************************************************************
	2	* Copyright(c) 2007-2009, 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	#include <TMath.h>
	17	#include <TH1F.h>
	18	#include <TFile.h>
	19	#include "AliVertexingHFUtils.h"
	20	#include "AliEventPlaneResolutionHandler.h"
	21	#include "AliLog.h"
	22
	23	/* $Id$ */
	24
	25	///////////////////////////////////////////////////////////////////
	26	// //
	27	// Class with functions useful for different D2H analyses //
	28	// - event plane resolution //
	29	// - <pt> calculation with side band subtraction //
	30	// - tracklet multiplicity calculation //
	31	// Origin: F.Prino, Torino, prino@to.infn.it //
	32	// //
	33	///////////////////////////////////////////////////////////////////
	34
	35	ClassImp(AliEventPlaneResolutionHandler)
	36
	37	//______________________________________________________________________
	38	AliEventPlaneResolutionHandler::AliEventPlaneResolutionHandler():TObject(),
	39	fEventPlane(kTPCPosEta),
	40	fResolOption(kTwoEtaSub),
	41	fMinCent(30.),
	42	fMaxCent(50.),
	43	fCorrHistoName1(""),
	44	fCorrHistoName2(""),
	45	fCorrHistoName3(""),
	46	fNsubevents(2),
	47	fExtrapToFull(kFALSE),
c60a4ab5	48	fWeight(0),
f5e6ce8f	49	fRootFileName("$ALICE_ROOT/PWGHF/vertexingHF/charmFlow/EventPlaneResolutionHistos.root")
	50	{
	51	// Default contructor
c60a4ab5	52	InitializeWeights();
f5e6ce8f	53	}
	54	//______________________________________________________________________
	55	AliEventPlaneResolutionHandler::AliEventPlaneResolutionHandler(TString filename):TObject(),
	56	fEventPlane(kTPCPosEta),
	57	fResolOption(kTwoEtaSub),
	58	fMinCent(30.),
	59	fMaxCent(50.),
	60	fCorrHistoName1(""),
	61	fCorrHistoName2(""),
	62	fCorrHistoName3(""),
	63	fNsubevents(2),
	64	fExtrapToFull(kFALSE),
c60a4ab5	65	fWeight(0),
f5e6ce8f	66	fRootFileName(filename.Data())
	67	{
	68	// Standard contructor
c60a4ab5	69	InitializeWeights();
f5e6ce8f	70	}
	71	//______________________________________________________________________
	72	Double_t AliEventPlaneResolutionHandler::GetEventPlaneResolution(Double_t minCent, Double_t maxCent){
	73	// method to compute the event plane resolution starting from correlations
	74
	75	TFile* inputFile=new TFile(fRootFileName.Data());
	76	if(!inputFile \|\| (inputFile && !inputFile->IsOpen())){
	77	printf("Root file not opened --- Return dummy resolution value\n");
	78	return -1.;
	79	}
	80	Bool_t isOk=SetHistoNames();
	81	if(!isOk){
	82	printf("Return dummy resolution value\n");
	83	return -1.;
	84	}
	85
	86	Int_t minCentrTimesTen=(Int_t)(minCent*10);
	87	Int_t maxCentrTimesTen=(Int_t)(maxCent*10);
	88	Int_t ncBin=minCentrTimesTen/25;
c60a4ab5	89	Int_t whichWei=fWeight;
	90	if(minCent>50.5 \|\| maxCent>50.5) fWeight=0;
	91	Double_t weight=-1.;
	92	if(fWeight==1) weight=fNcoll[ncBin];
	93	else if(fWeight==2 && minCent>29.5 && maxCent<50.5) weight=fYield24[ncBin];
	94	else if(fWeight==3 && minCent>29.5 && maxCent<50.5) weight=fYield46[ncBin];
	95	else if(fWeight==4 && minCent>29.5 && maxCent<50.5) weight=fYield612[ncBin];
f5e6ce8f	96	TH1F* hevpls2=0x0;
	97	TH1F* hevpls3=0x0;
	98	TH1F* hevpls1=(TH1F*)inputFile->Get(Form("%sCentr%d_%d",fCorrHistoName1.Data(),minCentrTimesTen,minCentrTimesTen+25));
	99	hevpls1->SetName(Form("%sCentr%d_%d",fCorrHistoName1.Data(),minCentrTimesTen,maxCentrTimesTen));
c60a4ab5	100	if(weight>0.) hevpls1->Scale(weight);
f5e6ce8f	101	if(fNsubevents==3){
	102	hevpls2=(TH1F*)inputFile->Get(Form("%sCentr%d_%d",fCorrHistoName2.Data(),minCentrTimesTen,minCentrTimesTen+25));
	103	hevpls2->SetName(Form("%sCentr%d_%d",fCorrHistoName2.Data(),minCentrTimesTen,maxCentrTimesTen));
	104	hevpls3=(TH1F*)inputFile->Get(Form("%sCentr%d_%d",fCorrHistoName3.Data(),minCentrTimesTen,minCentrTimesTen+25));
	105	hevpls3->SetName(Form("%sCentr%d_%d",fCorrHistoName2.Data(),minCentrTimesTen,maxCentrTimesTen));
c60a4ab5	106	if(weight>0.){
	107	hevpls2->Scale(weight);
	108	hevpls3->Scale(weight);
f5e6ce8f	109	}
	110	}
	111	for(Int_t iCentr=minCentrTimesTen+25; iCentr<maxCentrTimesTen; iCentr+=25){
	112	ncBin=iCentr/25;
c60a4ab5	113	weight=-1.;
	114	if(fWeight==1) weight=fNcoll[ncBin];
	115	else if(fWeight==2 && minCent>29.5 && maxCent<50.5) weight=fYield24[ncBin];
	116	else if(fWeight==3 && minCent>29.5 && maxCent<50.5) weight=fYield46[ncBin];
	117	else if(fWeight==4 && minCent>29.5 && maxCent<50.5) weight=fYield612[ncBin];
f5e6ce8f	118	TH1F* htmp1=(TH1F*)inputFile->Get(Form("%sCentr%d_%d",fCorrHistoName1.Data(),iCentr,iCentr+25));
c60a4ab5	119	if(weight>0.) htmp1->Scale(weight);
f5e6ce8f	120	hevpls1->Add(htmp1);
	121	if(fNsubevents==3){
	122	TH1F* htmp2=(TH1F*)inputFile->Get(Form("%sCentr%d_%d",fCorrHistoName2.Data(),iCentr,iCentr+25));
	123	TH1F* htmp3=(TH1F*)inputFile->Get(Form("%sCentr%d_%d",fCorrHistoName3.Data(),iCentr,iCentr+25));
c60a4ab5	124	if(weight>0.){
	125	htmp2->Scale(weight);
	126	htmp3->Scale(weight);
f5e6ce8f	127	}
	128	hevpls2->Add(htmp2);
	129	hevpls3->Add(htmp3);
	130	}
	131	}
c60a4ab5	132	fWeight=whichWei;
f5e6ce8f	133	Double_t resol=1.;
	134	if(fNsubevents==2){
	135	if(fExtrapToFull) resol=AliVertexingHFUtils::GetFullEvResol(hevpls1);
	136	else resol=AliVertexingHFUtils::GetSubEvResol(hevpls1);
	137	}else if(fNsubevents==3){
	138	Double_t correl1=hevpls1->GetMean();
	139	Double_t correl2=hevpls2->GetMean();
	140	Double_t correl3=hevpls3->GetMean();
	141	resol=TMath::Sqrt(correl1*correl2/correl3);
	142	}
	143	inputFile->Close();
	144	return resol;
	145	}
	146	//______________________________________________________________________
	147	Bool_t AliEventPlaneResolutionHandler::SetHistoNames(){
	148	// set names for histograms to be read
	149	if(fEventPlane==kTPCPosEta && fResolOption==kTwoRandSub){
	150	fCorrHistoName1="hCorrelRandSubTPCPosEta";
	151	fNsubevents=2;
	152	fExtrapToFull=kTRUE;
	153	}else if(fEventPlane==kTPCPosEta && fResolOption==kTwoChargeSub){
	154	fCorrHistoName1="hCorrelChargeSubTPCPosEta";
	155	fNsubevents=2;
	156	fExtrapToFull=kTRUE;
	157	}else if(fEventPlane==kTPCPosEta && fResolOption==kTwoEtaSub){
	158	fCorrHistoName1="hCorrelTPCEtaPosTPCEtaNeg";
	159	fNsubevents=2;
	160	fExtrapToFull=kFALSE;
	161	}else if(fEventPlane==kTPCPosEta && fResolOption==kThreeSub){
	162	fCorrHistoName1="hCorrelTPCPosEtaV0A";
	163	fCorrHistoName2="hCorrelTPCPosEtaV0C";
	164	fCorrHistoName3="hCorrelV0AV0C";
	165	fNsubevents=3;
	166	fExtrapToFull=kFALSE;
	167	}else if(fEventPlane==kTPCFullEta && fResolOption==kTwoRandSub){
	168	fCorrHistoName1="hCorrelRandSubTPCFullEta";
	169	fNsubevents=2;
	170	fExtrapToFull=kTRUE;
	171	}else if(fEventPlane==kTPCFullEta && fResolOption==kTwoEtaSub){
	172	fCorrHistoName1="hCorrelTPCEtaPosTPCEtaNeg";
	173	fNsubevents=2;
	174	fExtrapToFull=kTRUE;
	175	}else if(fEventPlane==kTPCFullEta && fResolOption==kThreeSub){
	176	fCorrHistoName1="hCorrelTPCFullEtaV0A";
	177	fCorrHistoName2="hCorrelTPCFullEtaV0C";
	178	fCorrHistoName3="hCorrelV0AV0C";
	179	fNsubevents=3;
	180	fExtrapToFull=kFALSE;
	181	}else if(fEventPlane==kVZERO && fResolOption==kThreeSub){
	182	fCorrHistoName1="hCorrelTPCEtaPosV0";
	183	fCorrHistoName2="hCorrelTPCEtaNegV0";
	184	fCorrHistoName3="hCorrelTPCEtaPosTPCEtaNeg";
	185	fNsubevents=3;
	186	fExtrapToFull=kFALSE;
	187	}else if(fEventPlane==kVZERO && fResolOption==kThreeSubTPCGap){
	188	fCorrHistoName1="hCorrelTPCEtaGt02V0";
	189	fCorrHistoName2="hCorrelTPCEtaLt02V0";
	190	fCorrHistoName3="hCorrelTPCEtaGt02TPCEtaLt02";
	191	fNsubevents=3;
	192	fExtrapToFull=kFALSE;
	193	}else if(fEventPlane==kVZEROA && fResolOption==kThreeSub){
	194	fCorrHistoName1="hCorrelTPCFullEtaV0A";
	195	fCorrHistoName2="hCorrelV0AV0C";
	196	fCorrHistoName3="hCorrelTPCFullEtaV0C";
197	fNsubevents=3;
198	fExtrapToFull=kFALSE;
199	}else if(fEventPlane==kVZEROC && fResolOption==kThreeSub){
200	fCorrHistoName1="hCorrelTPCFullEtaV0C";
201	fCorrHistoName2="hCorrelV0AV0C";
202	fCorrHistoName3="hCorrelTPCFullEtaV0A";
203	fNsubevents=3;
204	fExtrapToFull=kFALSE;
205	}else{
206	printf("Not possible to compute the resolution with the required option\n");
207	return kFALSE;
208	}
209	return kTRUE;
210	}
211	//______________________________________________________________________
c60a4ab5	212	void AliEventPlaneResolutionHandler::InitializeWeights(){
f5e6ce8f	213	// Ncoll values in 2.5% wide classes
	214	Double_t ncoll[20]={1790.77,1578.44,1394.82,1236.17,1095.08,
	215	969.86,859.571,759.959,669.648,589.588,
	216	516.039,451.409,392.853,340.493,294.426,
	217	252.385,215.484,183.284,155.101,130.963};
c60a4ab5	218	Double_t y24[20]={1,1,1,1,
	219	1,1,1,1,1,1,1,1,
	220	446,402,339,277,218,185,184,104};
	221
	222	Double_t y46[20]={1,1,1,1,
	223	1,1,1,1,1,1,1,1,
	224	49,150,127,105,88,71,71,53};
	225	Double_t y612[20]={1,1,1,1,
	226	1,1,1,1,1,1,1,1,
	227	89,92,79,66,57,59,38,58};
	228	for(Int_t i=0; i<20; i++){
	229	fNcoll[i]=ncoll[i];
	230	fYield24[i]=y24[i];
	231	fYield46[i]=y46[i];
	232	fYield612[i]=y612[i];
	233	}
f5e6ce8f	234	}