calibration/ATRI/voltageCalib/getADCmVConversionTest.cxx
00001 00002 00003 00004 00005 00006 //AraRoot Includes 00007 #include "UsefulAtriStationEvent.h" 00008 #include "RawAtriStationEvent.h" 00009 #include "araSoft.h" 00010 00011 00012 //Root Includes 00013 #include "TTree.h" 00014 #include "TFile.h" 00015 #include "TH1.h" 00016 #include "TF1.h" 00017 #include "TTree.h" 00018 #include "TMath.h" 00019 #include "TCanvas.h" 00020 #include "TChain.h" 00021 #include "TGraph.h" 00022 #include "TMultiGraph.h" 00023 00024 00025 //Standard Includes 00026 #include <iostream> 00027 #include <fstream> 00028 00029 00030 #define MAX_ADC 4096 00031 00032 //Prototype Functions 00033 Int_t getADCmVConversion(char *baseName, Int_t runLow, Int_t runHigh, Int_t thisBlock, Int_t thisDda); 00034 TGraph *getZeroSubGraph(TGraph *gr, Int_t factor); 00035 TGraph *addGraphs(TGraph *grOne, TGraph *grTwo); 00036 00037 int main(int argc, char **argv) 00038 { 00039 Int_t runLow, runHigh, block, dda; 00040 char baseName[FILENAME_MAX]; 00041 00042 if(argc<6) { 00043 std::cerr << "Usage: " << argv[0] << " <baseDir> <runLow> <runHigh> <block> <dda> \n"; 00044 return 1; 00045 } 00046 sprintf(baseName, argv[1]); 00047 runLow=atoi(argv[2]); 00048 runHigh=atoi(argv[3]); 00049 block=atoi(argv[4]); 00050 dda=atoi(argv[5]); 00051 getADCmVConversion(baseName, runLow, runHigh, block, dda); 00052 00053 } 00054 00055 00056 00057 00058 00059 00060 Int_t getADCmVConversion(char *baseName, Int_t runLow, Int_t runHigh, Int_t thisBlock, Int_t thisDda){ 00061 00062 TChain *chain = new TChain("maxPosNegTree"); 00063 char runName[FILENAME_MAX]; 00064 for(int run=runLow; run<=runHigh;run++){ 00065 sprintf(runName, "%s/root/run%i/adcSampleBlock.root", baseName, run); 00066 chain->Add(runName); 00067 } 00068 char outFileName[FILENAME_MAX]; 00069 sprintf(outFileName, "%s/root/voltageCalibThursday/ADCmVConversionTest_block%idda%i.root", baseName, thisBlock, thisDda); 00070 00071 TFile *outFile = new TFile(outFileName, "RECREATE"); 00072 00073 Int_t dda=0, chan=0, sample=0, block=0, posNeg=0, ndf=0; 00074 Double_t maxPosValue=0, maxNegValue=0, voltage=0; 00075 Double_t p0=0, p1=0, p2=0, p3=0, p4=0, chiSq=0; 00076 00077 00078 00079 char exp[100]; 00080 char cut[100]; 00081 char grName[100]; 00082 TGraph *grTemp; 00083 TCanvas *tempCan; 00084 00085 TF1* funcPos = new TF1("myPosPol3", "pol3"); 00086 TF1* funcNeg = new TF1("myNegPol1", "pol1"); 00087 funcPos->SetParLimits(1, 0,2); 00088 funcPos->SetParLimits(2, -1,+1); 00089 funcPos->SetParLimits(3, -1,+1); 00090 funcNeg->SetParLimits(1, 0,+2); 00091 funcPos->FixParameter(0,0); 00092 funcNeg->FixParameter(0,0); 00093 00094 00095 00096 00097 TTree *fitTree = new TTree("fitTree", "Tree of fit parameters"); 00098 fitTree->Branch("dda", &dda, "dda/I"); 00099 fitTree->Branch("chan", &chan, "chan/I"); 00100 fitTree->Branch("sample", &sample, "sample/I"); 00101 fitTree->Branch("block", &block, "block/I"); 00102 fitTree->Branch("ndf", &ndf, "ndf/I"); 00103 fitTree->Branch("posNeg", &posNeg, "posNeg/I"); 00104 00105 fitTree->Branch("p0", &p0, "p0/D"); 00106 fitTree->Branch("p1", &p1, "p1/D"); 00107 fitTree->Branch("p2", &p2, "p2/D"); 00108 fitTree->Branch("p3", &p3, "p3/D"); 00109 fitTree->Branch("chiSq", &chiSq, "chiSq/D"); 00110 00111 00112 00113 block=thisBlock; 00114 dda=thisDda; 00115 // for(dda=0;dda<1;dda++){ 00116 for(chan=0;chan<1;chan++){ 00117 // fprintf(stderr, "dda %i chan %i ", dda, chan); 00118 for(sample=40;sample<45;sample++){ 00119 00120 if((dda==1 || dda==2)&&(chan>3)) continue; 00121 //Reset starting point of parameters 00122 funcPos->SetParameter(1, 1.0); 00123 funcPos->SetParameter(2, -0.001); 00124 funcPos->SetParameter(3, 0); 00125 funcNeg->SetParameter(1, 0.55); 00126 00127 // fprintf(stderr, "%i ", sample); 00128 sprintf(exp, "voltage/2:maxPosValue"); 00129 sprintf(cut, "dda==%i&&chan==%i&&block==%i&&sample==%i", dda, chan, block, sample); 00130 sprintf(grName, "grPos_%i_%i_%i_%i", dda, chan, block, sample); 00131 tempCan = new TCanvas(); 00132 chain->Draw(exp, cut); 00133 grTemp = (TGraph*)tempCan->GetPrimitive("Graph"); 00134 // TGraph *grPosTemp = (TGraph*) grTemp->Clone(); 00135 TGraph *grPosTemp = getZeroSubGraph(grTemp, +1); 00136 00137 grPosTemp->SetName(grName); 00138 grPosTemp->Fit("myPosPol3", "QBW"); 00139 p0=funcPos->GetParameter(0); 00140 p1=funcPos->GetParameter(1); 00141 p2=funcPos->GetParameter(2); 00142 p3=funcPos->GetParameter(3); 00143 chiSq =funcPos->GetChisquare(); 00144 ndf = funcPos->GetNDF(); 00145 p4=0; 00146 posNeg=+1; 00147 printf("%i %i %i %i %f %f\n", dda, chan, block, sample, p1, p2); 00148 fitTree->Fill(); 00149 delete tempCan; 00150 00151 00152 tempCan = new TCanvas(); 00153 sprintf(exp, "-1*voltage/2:maxNegValue"); 00154 sprintf(grName, "grNeg_%i_%i_%i_%i", dda, chan, block, sample); 00155 chain->Draw(exp, cut); 00156 grTemp = (TGraph*)tempCan->GetPrimitive("Graph"); 00157 // TGraph *grNegTemp = (TGraph*) grTemp->Clone(); 00158 TGraph *grNegTemp = getZeroSubGraph(grTemp, -1); 00159 00160 grNegTemp->SetName(grName); 00161 grNegTemp->Fit("myNegPol1", "QBW"); 00162 p0=funcNeg->GetParameter(0); 00163 p1=funcNeg->GetParameter(1); 00164 chiSq=funcNeg->GetChisquare(); 00165 p2=0; 00166 p3=0; 00167 p4=0; 00168 // p0=grNegTemp->GetFunction("myNegPol1")->GetParameter(0); 00169 // p1=grNegTemp->GetFunction("myNegPol1")->GetParameter(1); 00170 // chiSq = grNegTemp->GetFunction("myNegPol1")->GetChisquare(); 00171 00172 posNeg=-1; 00173 printf("%i %i %i %i %f %f\n", dda, chan, block, sample, p1, p2); 00174 fitTree->Fill(); 00175 delete tempCan; 00176 00177 00178 outFile->cd(); 00179 grPosTemp->Write(); 00180 grNegTemp->Write(); 00181 // grPosNeg->Write(); 00182 00183 // delete grPosNeg; 00184 delete grNegTemp; 00185 delete grPosTemp; 00186 00187 00188 } 00189 fprintf(stderr, "\n"); 00190 } 00191 // } 00192 00193 outFile->cd(); 00194 fitTree->Write(); 00195 00196 fprintf(stderr, "Done!\n"); 00197 00198 delete funcPos; 00199 delete funcNeg; 00200 00201 00202 return 0; 00203 00204 } 00205 00206 00207 00208 TGraph *getZeroSubGraph(TGraph *gr, Int_t factor){ 00209 Int_t numEntries = gr->GetN(); 00210 Double_t *yVals=gr->GetY(); 00211 Double_t *xVals=gr->GetX(); 00212 00213 if(numEntries<=0) return NULL; 00214 Double_t *xValsNew = new Double_t[numEntries]; 00215 Double_t xValZero=0; 00216 00217 for(int entry=0;entry<numEntries;entry++){ 00218 if(yVals[entry]==0){ 00219 xValZero=xVals[entry]; 00220 break; 00221 } 00222 } 00223 for(int entry=0;entry<numEntries;entry++){ 00224 Double_t newX = TMath::Sqrt( TMath::Power(xVals[entry],2) - TMath::Power(xValZero,2) ); 00225 xValsNew[entry] = factor*newX; 00226 } 00227 00228 TGraph *grNew = new TGraph(numEntries, xValsNew, yVals); 00229 00230 00231 delete [] xValsNew; 00232 00233 return grNew; 00234 00235 } 00236 00237 TGraph *addGraphs(TGraph *grOne, TGraph *grTwo){ 00238 Int_t numEntriesOne = grOne->GetN(); 00239 Int_t numEntriesTwo = grTwo->GetN(); 00240 Double_t *yVals1 = grOne->GetY(); 00241 Double_t *xVals1 = grOne->GetX(); 00242 Double_t *yVals2 = grTwo->GetY(); 00243 Double_t *xVals2 = grTwo->GetX(); 00244 00245 Double_t *yValsNew = new Double_t[numEntriesOne+numEntriesTwo]; 00246 Double_t *xValsNew = new Double_t[numEntriesOne+numEntriesTwo]; 00247 00248 for(int entry = 0;entry<numEntriesOne;entry++){ 00249 xValsNew[entry] = xVals1[entry]; 00250 yValsNew[entry] = yVals1[entry]; 00251 } 00252 for(int entry = 0;entry<numEntriesTwo;entry++){ 00253 xValsNew[entry+numEntriesOne] = xVals2[entry]; 00254 yValsNew[entry+numEntriesOne] = yVals2[entry]; 00255 } 00256 00257 TGraph *grNew = new TGraph(numEntriesOne+numEntriesTwo, xValsNew, yValsNew); 00258 00259 delete [] yValsNew; 00260 delete [] xValsNew; 00261 00262 return grNew; 00263 }
Generated on Tue Jul 16 16:58:02 2013 for ARA ROOT v3.10 Software by
