ROCsoft: beamtime/cern-oct11/go4/AnalysisMacros/position_resolution.cxx (r4864/r2584) Source File

Go to the documentation of this file.
00001 #ifndef __CINT__
00002 
00003 #include <iostream>
00004 #include <string>
00005 #include <typeinfo>
00006 #include <iostream>
00007 #include <vector>
00008 
00009 
00010 #include "Riostream.h"
00011 #include "TRint.h"
00012 #include "TROOT.h"
00013 #include "TStyle.h"
00014 #include "Riostream.h"
00015 
00016 
00017 #include "TH2.h"
00018 #include "TH1.h"
00019 #include "TF1.h"
00020 
00021 #include "TGraph.h"
00022 #include "TMultiGraph.h"
00023 #include "TCanvas.h"
00024 #include "TProfile.h"
00025 #include "TStopwatch.h"
00026 
00027 #include "TFile.h"
00028 #include "TTree.h"
00029 #include "TBranch.h"
00030 #include "TMath.h"
00031 #include "TLegend.h"
00032 #include "TColor.h"
00033 #include "TText.h"
00034 #include "TKey.h"
00035 #include "TLatex.h"
00036 #include "TCanvas.h"
00037 
00038 #include "TRocEvent.h"
00039 #include "roc/Message.h"
00040 #include "TSpadicEvent.h"
00041 
00042 
00043 #include "TMbsCrateEvent.h"
00044 #include "TBeamMonitorEvent.h"
00045 #include "TCernOct11UnpackEvent.h"
00046 
00047 //#include "langaus.C"
00048 
00049 #endif
00050 
00051 #define NUM_SPADIC_CHA             8
00052 #define SPADIC_TRACE_SIZE         45
00053 #define SPADIC_ADC_MAX           255
00054 #define noBaselineTB               5 
00055 
00056 using namespace std;
00057 
00058 void Statusbar(Int_t i, Int_t n) {
00059   if (int(i * 100 / float(n)) - int((i-1) * 100 / float(n)) >= 1) {
00060     if (int(i * 100 / float(n)) == 1 || i == 1 || i == 0) 
00061       cout << "[" << flush;
00062     cout << "-" << flush;
00063     if (int(i * 10 / float(n)) - int((i-1) * 10 / float(n)) >= 1) 
00064       cout << "|";
00065     if (int(i * 100 / float(n)) >=99) 
00066       cout << "]" <<endl;
00067   }
00068 }
00069 
00070 void BaselineCorrection(Double_t pulse[NUM_SPADIC_CHA][SPADIC_TRACE_SIZE]) {
00071   for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {
00072     //cout << endl;
00073     Float_t baseline(0);        
00074     for (Int_t bin = 0; bin < noBaselineTB; bin++) {    
00075       baseline += pulse[ch][bin];
00076     }
00077     baseline /= Float_t(noBaselineTB);
00078     for (Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++) {       
00079       pulse[ch][bin] -= baseline;
00080     }
00081   }
00082 }
00083 Bool_t Distance2TrackFFM(Double_t x[8], Double_t distance[8]/*, Double_t mAlig[8][8], Double_t bAlig[8][8] /*, Bool_t MS = true*/) {
00084   //Alignment correction
00085  
00086   Double_t xPos[8] = {0.0};
00087   Int_t nPos = 0;
00088   for (Int_t i = 4; i < 8; i++) {
00089     xPos[i] = x[i];
00090     if (xPos[i] != 0.0) nPos++;
00091   }
00092   if (nPos < 4) return false;
00093   /*
00094     for (Int_t i = 5; i < 8; i++) {
00095     xPos[i] -= (mAlig[i][5] * xPos[i] + bAlig[i][5]); //best    // use suid = 0 as fixed reference     
00096     }
00097   
00098     for (Int_t i = 4; i < 8; i++) {
00099     xPos[i] -= (mAlig[i][i] * xPos[i] + bAlig[i][i]);// use suid = i as dynamic self reference 
00100     }
00101   */
00102   Double_t zMean(0.0), xMean(0.0), m(0.0), b(0.0), SUMzx(0.0), SUMzz(0.0);
00103   //Least squares fit: x = m * z + b
00104   //                  MS  Ms   Ms   MS  F    F    F    F
00105   Double_t zPos[8] = {1, 321, 639, 959, 1, 236, 634, 817};
00106   /*
00107   // exclude point of interest from fit
00108   for (Int_t j = 4; j < 8; j++) {
00109   Int_t count = 0;
00110   for (Int_t i = 4; i < 8; i++) {
00111   if (i != j) {
00112   zMean += zPos[i];
00113   xMean += xPos[i];
00114   count++;
00115   }
00116   }
00117   zMean /= Float_t(count);
00118   xMean /= Float_t(count);
00119   for (Int_t i = 4; i < 8; i++) {
00120   if (i != j) {
00121   SUMzx += (zPos[i] - zMean) * (xPos[i] - xMean);
00122   SUMzz += (zPos[i] - zMean) * (zPos[i] - zMean);
00123   }
00124   }
00125   m = SUMzx / SUMzz;
00126   b = xMean - m * zMean;
00127   distance[j] = xPos[j] - (m * zPos[j] + b);
00128   }
00129   */
00130   //include point of interest from fit
00131   Int_t count = 0;
00132   for (Int_t i = 4; i < 8; i++) {     
00133     zMean += zPos[i];
00134     xMean += xPos[i];
00135     count++;
00136   }
00137   zMean /= Float_t(count);
00138   xMean /= Float_t(count);
00139   for (Int_t i = 4; i < 8; i++) { 
00140     SUMzx += (zPos[i] - zMean) * (xPos[i] - xMean);
00141     SUMzz += (zPos[i] - zMean) * (zPos[i] - zMean);
00142   }
00143   m = SUMzx / SUMzz;
00144   b = xMean - m * zMean;
00145   //hmTrack->Fill(m);
00146   //hbTrack->Fill(b);
00147   for (Int_t i = 4; i < 8; i++) {
00148     distance[i] = xPos[i] - (m * zPos[i] + b);
00149   }
00150   return true;
00151 }
00152 Bool_t Distance2TrackMS(Double_t x[8], Double_t distance[8]/*, TH1I* hmTrack, TH1I* hbTrack/*, Double_t mAlig[8][8], Double_t bAlig[8][8] /*, Bool_t MS = true*/) {
00153   //Alignment correction
00154   // use suid = 0 as fixed reference
00155   Double_t xPos[4] = {0.0};
00156   Int_t nPos = 0;
00157   for (Int_t i = 0; i< 4; i++) {
00158     xPos[i] = x[i];
00159     if (xPos[i] != 0.0) nPos++;
00160   }
00161   if (nPos < 4) return false;
00162   /*
00163     for (Int_t i = 1; i < 4; i++) {
00164     xPos[i] -= (mAlig[i][0] * xPos[i] + bAlig[i][0]); //best    // use suid = 0 as fixed reference 
00165     } 
00166   
00167     for (Int_t i = 0; i < 4; i++) {
00168     xPos[i] -= (mAlig[i][i] * xPos[i] + bAlig[i][i]);// use suid = i as dynamic self reference 
00169     } 
00170   */
00171   Double_t zMean(0.0), xMean(0.0), m(0.0), b(0.0), SUMzx(0.0), SUMzz(0.0);
00172   //Least squares fit: x = m * z + b
00173   Double_t zPos[4] = {1, 321, 639, 959};
00174   // exclude point of interest from fit
00175   /*
00176     for (Int_t j = 0; j < 4; j++) {
00177     Int_t count = 0;
00178     for (Int_t i = 0; i < 4; i++) {
00179     if (i != j) {
00180     zMean += zPos[i];
00181     xMean += xPos[i];
00182     count++;
00183     }
00184     }   
00185     zMean /= Float_t(count);
00186     xMean /= Float_t(count);
00187     for (Int_t i = 0; i < 4; i++) {
00188     if (i != j) {
00189     SUMzx += (zPos[i] - zMean) * (xPos[i] - xMean);
00190     SUMzz += (zPos[i] - zMean) * (zPos[i] - zMean);
00191     }
00192     }
00193     m = SUMzx / SUMzz;
00194     b = xMean - m * zMean;
00195     mTrack->Fill(m);
00196     bTrack->Fill(b);
00197     distance[j] = xPos[j] - (m * zPos[j] + b);
00198     }
00199   */
00200   //include point of interest from fit
00201   Int_t count = 0;
00202   for (Int_t i = 0; i < 4; i++) {     
00203     zMean += zPos[i];
00204     xMean += xPos[i];
00205     count++;
00206   }
00207   zMean /= Float_t(count);
00208   xMean /= Float_t(count);
00209   for (Int_t i = 0; i < 4; i++) { 
00210     SUMzx += (zPos[i] - zMean) * (xPos[i] - xMean);
00211     SUMzz += (zPos[i] - zMean) * (zPos[i] - zMean);
00212   }
00213   m = SUMzx / SUMzz;
00214   b = xMean - m * zMean;
00215   //hmTrack->Fill(m);
00216   //hbTrack->Fill(b);
00217   for (Int_t i = 0; i < 4; i++) {
00218     distance[i] = xPos[i] - (m * zPos[i] + b);
00219   }
00220   return true;
00221 }
00222 
00223 void position_resolution(TString nameFileSet, Bool_t debug=0) {
00224 
00225   Bool_t msMerge = false;
00226 
00227   ifstream fileSet;
00228   fileSet.open(nameFileSet.Data());
00229 
00230   if (fileSet.is_open()) {
00231     TString ofname = "data2011/Analysis/PR_" + TString(nameFileSet(0,nameFileSet.Sizeof()-5)).Append(".root");
00232     TFile outputFile(ofname,"RECREATE","output of position_resolution.cxx");
00233     TString ifnamefit = "data2011/Analysis/" + TString(nameFileSet(0,nameFileSet.Sizeof()-5)).Append(".root");
00234     TFile inputFileFit(ifnamefit,"READ");
00235     TString mergeName[2] = {"NotMerged","Merged"};
00236     TString filterName[4] = {"wOwU","woOwU","woOwoU","wOwoU"}; // w == with; wo == without; O == overflow; U == undershoot
00237     TString recoSetName[3] = {"3PadInteg","3PadIntegWindow","3PadAmpl"}; 
00238     TString folder = "data2011/TTrees/";
00239     TString picEx = ".png";
00240     const Int_t usedSusibos = 4;
00241     const Int_t recoSets = 3;
00242     const Int_t nMerge = 2;
00243     const Int_t nFilter = 4;
00244     const Int_t nIteration = 5;//0;
00245 
00246     Double_t mTrack[usedSusibos][nIteration] = {{0.0}};
00247     Double_t bTrack[usedSusibos][nIteration] = {{0.0}};
00248     //const Int_t noBaselineTB = 5;
00249     Int_t topSuId[4]            = {         1,         11,        3,        15};
00250     Int_t usedSuId[8] = {         2,         12,        6,        18,         4,         5,        16,        17  };
00251     Float_t padWidth[8] = {       5,          5,        8,         8,         5,         5,         5,         5  };
00252   
00253     //Int_t allPion[usedSusibos][nMerge] = {{0}};
00254     //Int_t allElectron[usedSusibos][nMerge] = {{0}};
00255     //Int_t allGoodEvent[usedSusibos][nMerge] = {{0}};
00256     Int_t allEvents = 0;
00257 
00258     // Fit functions
00259     TF1 *gausFit = new TF1("gausFit","gaus",-10,10);
00260     gausFit->SetLineColor(2);
00261     TF1 *alignmentFit= new TF1("alignmentFit","[0]*x+[1]",0,60);
00262     alignmentFit->SetLineColor(2);
00263     TF1 *prFit = new TF1("prFit","gaus",-1,1);
00264     prFit->SetLineColor(2);
00265 
00266     // Histogramms
00267     TString name, title, newtitle, appending;
00268 
00269     TProfile *prResult[recoSets][nMerge][nFilter];
00270 
00271     TProfile *prwResult[recoSets][nMerge][nFilter];
00272 
00273     for (Int_t rs = 0; rs < recoSets; rs++) {
00274       for (Int_t imerge = 0; imerge < nMerge; imerge++) {
00275         for (Int_t ifilter = 0; ifilter < nFilter; ifilter++) {
00276           appending.Form("_%s_%s_%s",recoSetName[rs].Data(),mergeName[imerge].Data(),filterName[ifilter].Data());
00277           prResult[rs][imerge][ifilter] = new TProfile("prResult" + appending,"prResult" + appending,usedSusibos,-0.5,usedSusibos-0.5);
00278           prResult[rs][imerge][ifilter]->SetMarkerStyle(04);
00279           //prResult->SetXTitle("distance to track [mm]");
00280           prResult[rs][imerge][ifilter]->SetYTitle("position resolution [mm]");
00281           //prResult->SetErrorOption("g");
00282 
00283           prwResult[rs][imerge][ifilter] = new TProfile("prwResult" + appending,"prwResult" + appending,usedSusibos,-0.5,usedSusibos-0.5);
00284           prwResult[rs][imerge][ifilter]->SetMarkerStyle(04);
00285           //prwResult->SetXTitle("distance to track [mm]");
00286           prwResult[rs][imerge][ifilter]->SetYTitle("position resolution [mm]");
00287           //prwResult->SetErrorOption("g");
00288           for (Int_t suid = 0; suid < usedSusibos; suid++) {
00289             title.Form("Susibo %i",usedSuId[suid]);
00290             prResult[rs][imerge][ifilter]->GetXaxis()->SetBinLabel( suid+1, title);
00291             prwResult[rs][imerge][ifilter]->GetXaxis()->SetBinLabel(suid+1, title);
00292           }
00293         }
00294       }
00295     }
00296     TProfile *prf[usedSusibos]; 
00297     TProfile *chamberAlignment[usedSusibos][usedSusibos]; 
00298     TH1I *beamProfile[usedSusibos]; 
00299     TH1I *wbeamProfile[usedSusibos]; 
00300 
00301     TProfile *posCorr[usedSusibos][usedSusibos];
00302     TProfile *wposCorr[usedSusibos][usedSusibos];
00303     TH2I *posCorr2D[usedSusibos][usedSusibos];
00304     TH2I *wposCorr2D[usedSusibos][usedSusibos];
00305 
00306     TProfile *trackPosCorr[usedSusibos][nIteration];
00307     TProfile *trackwPosCorr[usedSusibos][nIteration];
00308     TH2I *trackPosCorr2D[usedSusibos][nIteration];
00309     TH2I *trackwPosCorr2D[usedSusibos][nIteration];
00310 
00311     Double_t m[usedSusibos][usedSusibos] = {{0.0}};
00312     Double_t b[usedSusibos][usedSusibos] = {{0.0}};
00313     Double_t sigma[usedSusibos] = {0.0};
00314 
00315     TH1I *residuals[usedSusibos][usedSusibos];
00316     TH1I *wresiduals[usedSusibos][usedSusibos];
00317 
00318     TH1I *trackResiduals[usedSusibos][nIteration];
00319     TH1I *wtrackResiduals[usedSusibos][nIteration];
00320 
00321     TH1I *hmTrack[nIteration];
00322     TH1I *hbTrack[nIteration];
00323 
00324     TH1D *rPulseShape[NUM_SPADIC_CHA];
00325     TH1D *cPulseShape[NUM_SPADIC_CHA];
00326     Int_t lineColor[NUM_SPADIC_CHA] = {kBlue, kRed, kGreen, kYellow, kBlack, kCyan, kMagenta, kRed+2};
00327     for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {
00328       rPulseShape[ch] = new TH1D (title + appending,title + appending,SPADIC_TRACE_SIZE,0,1.8/*SPADIC_TRACE_SIZE*/);
00329       rPulseShape[ch]->SetXTitle("time [#mus]"/*"time-bin"*/);
00330       rPulseShape[ch]->SetYTitle("ADC value");
00331       rPulseShape[ch]->SetLineColor(lineColor[ch]);
00332       cPulseShape[ch] = new TH1D (title + appending,title + appending,SPADIC_TRACE_SIZE,0,1.8/*SPADIC_TRACE_SIZE*/);
00333       cPulseShape[ch]->SetXTitle("time [#mus]"/*"time-bin"*/);
00334       cPulseShape[ch]->SetYTitle("ADC value");
00335       cPulseShape[ch]->SetLineColor(lineColor[ch]);
00336     }
00337     for (Int_t iteration = 0; iteration < nIteration; iteration++) {
00338       title.Form("hmTrack_%i",iteration);
00339       name.Form("m track (%i)",iteration);
00340       hmTrack[iteration] = new TH1I(title, name,500,-5,5);
00341       title.Form("hbTrack_%i",iteration);
00342       name.Form("b track (%i)",iteration);
00343       hbTrack[iteration] = new TH1I(title, name,3000,-30,30);
00344     }
00345 
00346     for (Int_t suid = 0; suid < usedSusibos; suid++) {
00347       title.Form("PRF_Profile%i",usedSuId[suid]);
00348       newtitle.Form("1DPRF_Profile%i",usedSuId[suid]);
00349       prf[suid] = (TProfile*)inputFileFit.Get(title)->Clone(newtitle);
00350       prf[suid]->Fit("gausFit","RQ0");
00351       sigma[suid] = gausFit->GetParameter(2);
00352       if (debug) 
00353         printf("  Susibo %i   A %.2f  mean %.2f sigma %.2f\n",suid,gausFit->GetParameter(0),gausFit->GetParameter(1),gausFit->GetParameter(2));
00354 
00355       title.Form("beamProfile_%i",usedSuId[suid]);
00356       name.Form("beam profile (Sus%i)",usedSuId[suid]);
00357       beamProfile[suid] = new TH1I(title, name, NUM_SPADIC_CHA*padWidth[suid]*100,-0.5*padWidth[suid],(NUM_SPADIC_CHA-0.5)*padWidth[suid]);
00358       title.Form("reonstructed hit position [mm] (Susibo %i)",usedSuId[suid]);
00359       beamProfile[suid]->SetXTitle(title);
00360       beamProfile[suid]->SetYTitle("#");
00361 
00362       title.Form("wbeamProfile_%i",usedSuId[suid]);
00363       name.Form("weighted beam profile (Sus%i)",usedSuId[suid]);
00364       wbeamProfile[suid] = new TH1I(title, name, NUM_SPADIC_CHA*padWidth[suid]*100,-0.5*padWidth[suid],(NUM_SPADIC_CHA-0.5)*padWidth[suid]);
00365       title.Form("reonstructed weighted centre of gravity hit position [mm] (Susibo %i)",usedSuId[suid]);
00366       wbeamProfile[suid]->SetXTitle(title);
00367       wbeamProfile[suid]->SetYTitle("#");
00368 
00369       for (Int_t iteration = 0; iteration < nIteration; iteration++) {
00370         title.Form("trackPosCorr_%i_%i",usedSuId[suid],iteration);
00371         name.Form("track  (Sus%i) iteration %i",usedSuId[suid],iteration);
00372         trackPosCorr[suid][iteration] = new TProfile(title, name, NUM_SPADIC_CHA*padWidth[suid],-0.5*padWidth[suid],(NUM_SPADIC_CHA-0.5)*padWidth[suid]);
00373         trackPosCorr[suid][iteration]->SetMarkerStyle(04);
00374         title.Form("reconstructed hit position [mm] (Susibo %i)",usedSuId[suid]);
00375         trackPosCorr[suid][iteration]->SetXTitle(title);
00376         title.Form("distance to track [mm] (Susibo %i)",usedSuId[suid]);
00377         trackPosCorr[suid][iteration]->SetYTitle(title);
00378 
00379         title.Form("trackwPosCorr_%i_%i",usedSuId[suid],iteration);
00380         name.Form("track  (Sus%i) iteration %i",usedSuId[suid],iteration);
00381         trackwPosCorr[suid][iteration] = new TProfile(title, name, NUM_SPADIC_CHA*padWidth[suid],-0.5*padWidth[suid],(NUM_SPADIC_CHA-0.5)*padWidth[suid]);
00382         trackwPosCorr[suid][iteration]->SetMarkerStyle(04);
00383         title.Form("reconstructed hit position [mm] (Susibo %i)",usedSuId[suid]);
00384         trackwPosCorr[suid][iteration]->SetXTitle(title);
00385         title.Form("distance to track [mm] (Susibo %i)",usedSuId[suid]);
00386         trackwPosCorr[suid][iteration]->SetYTitle(title);
00387 
00388         title.Form("trackPosCorr2D_%i_%i",usedSuId[suid],iteration);
00389         name.Form("track  (Sus%i) iteration %i",usedSuId[suid],iteration);
00390         trackPosCorr2D[suid][iteration] = new TH2I(title, name, 
00391                                                    NUM_SPADIC_CHA*padWidth[suid]*10,-0.5*padWidth[suid],(NUM_SPADIC_CHA-0.5)*padWidth[suid],
00392                                                    2*NUM_SPADIC_CHA*padWidth[suid]*10,-(NUM_SPADIC_CHA-0.5)*padWidth[suid],(NUM_SPADIC_CHA-0.5)*padWidth[suid]);
00393         title.Form("reconstructed hit position [mm] (Susibo %i)",usedSuId[suid]);
00394         trackPosCorr2D[suid][iteration]->SetXTitle(title);
00395         title.Form("distance to track [mm] (Susibo %i)",usedSuId[suid]);
00396         trackPosCorr2D[suid][iteration]->SetYTitle(title);
00397         trackPosCorr2D[suid][iteration]->SetContour(99);
00398 
00399 
00400         title.Form("trackwPosCorr2D_%i_%i",usedSuId[suid],iteration);
00401         name.Form("track  (Sus%i) iteration %i",usedSuId[suid],iteration);
00402         trackwPosCorr2D[suid][iteration] = new TH2I(title, name, 
00403                                                     NUM_SPADIC_CHA*padWidth[suid]*10,-0.5*padWidth[suid],(NUM_SPADIC_CHA-0.5)*padWidth[suid],
00404                                                     2*NUM_SPADIC_CHA*padWidth[suid]*10,-(NUM_SPADIC_CHA-0.5)*padWidth[suid],(NUM_SPADIC_CHA-0.5)*padWidth[suid]);
00405 
00406         title.Form("reconstructed hit position [mm] (Susibo %i)",usedSuId[suid]);
00407         trackwPosCorr2D[suid][iteration]->SetXTitle(title);
00408         title.Form("distance to track [mm] (Susibo %i)",usedSuId[suid]);
00409         trackwPosCorr2D[suid][iteration]->SetYTitle(title);
00410         trackwPosCorr2D[suid][iteration]->SetContour(99);
00411 
00412 
00413         title.Form("trackResiduals_%i_%i",usedSuId[suid],iteration);
00414         name.Form("track residuals (Sus%i) iteration %i",usedSuId[suid],iteration);
00415         trackResiduals[suid][iteration] = new TH1I(title, name, padWidth[suid]*100, -padWidth[suid], padWidth[suid]);
00416         trackResiduals[suid][iteration]->SetLineColor(1);
00417         trackResiduals[suid][iteration]->SetLineStyle(2);
00418         title.Form("distance to track [mm] (Susibo %i)",usedSuId[suid]);
00419         trackResiduals[suid][iteration]->SetXTitle(title);
00420         trackResiduals[suid][iteration]->SetYTitle("#");
00421 
00422         title.Form("wtrackRsiduals_%i_%i",usedSuId[suid],iteration);
00423         name.Form("weighted track residuals (Sus%i) iteration %i",usedSuId[suid],iteration);
00424         wtrackResiduals[suid][iteration] = new TH1I(title, name, padWidth[suid]*100, -padWidth[suid], padWidth[suid]);
00425         wtrackResiduals[suid][iteration]->SetLineColor(1);
00426         wtrackResiduals[suid][iteration]->SetLineStyle(1);
00427         title.Form("distance to track [mm] (Susibo %i)",usedSuId[suid]);
00428         wtrackResiduals[suid][iteration]->SetXTitle(title);
00429         wtrackResiduals[suid][iteration]->SetYTitle("#");
00430       }
00431 
00432       for (Int_t suid2 = 0; suid2 < usedSusibos; suid2++) {
00433         title.Form("chamberAlignment_%i_%i",usedSuId[suid],usedSuId[suid2]);
00434         newtitle.Form("ChamberAlignment_%i_%i",usedSuId[suid],usedSuId[suid2]);
00435         chamberAlignment[suid][suid2] = (TProfile*)inputFileFit.Get(title)->Clone(newtitle);
00436         chamberAlignment[suid][suid2]->Fit("alignmentFit","Q","0",0.25*NUM_SPADIC_CHA*padWidth[suid],0.75*NUM_SPADIC_CHA*padWidth[suid]);
00437         m[suid][suid2] = alignmentFit->GetParameter(0);
00438         b[suid][suid2] = alignmentFit->GetParameter(1);
00439         if (debug) 
00440           printf("%.2f * x + %.2f\n",m[suid][suid2],b[suid][suid2]);
00441 
00442         title.Form("posCorr_%i_%i",usedSuId[suid],usedSuId[suid2]);
00443         name.Form("position correlation (Sus%i,Sus%i)",usedSuId[suid],usedSuId[suid2]);
00444         posCorr[suid][suid2] = new TProfile(title, name, NUM_SPADIC_CHA*padWidth[suid],-0.5*padWidth[suid],(NUM_SPADIC_CHA-0.5)*padWidth[suid]);
00445         posCorr[suid][suid2]->SetMarkerStyle(04);
00446         title.Form("reconstructed hit position [mm] (Susibo %i)",usedSuId[suid]);
00447         posCorr[suid][suid2]->SetXTitle(title);
00448         title.Form("distance of reconstructed hit position (Susibo %i - Sisibo %i) [mm]",usedSuId[suid],usedSuId[suid2]);
00449         posCorr[suid][suid2]->SetYTitle(title);
00450 
00451         title.Form("wposCorr_%i_%i",usedSuId[suid],usedSuId[suid2]);
00452         name.Form("weighted position correlation (Sus%i,Sus%i)",usedSuId[suid],usedSuId[suid2]);
00453         wposCorr[suid][suid2] = new TProfile(title, name, NUM_SPADIC_CHA*padWidth[suid],-0.5*padWidth[suid],(NUM_SPADIC_CHA-0.5)*padWidth[suid]);
00454         wposCorr[suid][suid2]->SetMarkerStyle(04);
00455         title.Form("reconstructed hit position [mm] (Susibo %i)",usedSuId[suid]);
00456         wposCorr[suid][suid2]->SetXTitle(title);
00457         title.Form("distance of reconstructed hit position (Susibo %i - Sisibo %i) [mm]",usedSuId[suid],usedSuId[suid2]);
00458         wposCorr[suid][suid2]->SetYTitle(title);
00459 
00460         title.Form("posCorr2D_%i_%i",usedSuId[suid],usedSuId[suid2]);
00461         name.Form("position correlation 2D (Sus%i,Sus%i)",usedSuId[suid],usedSuId[suid2]);
00462         posCorr2D[suid][suid2] = new TH2I(title, name, 
00463                                           NUM_SPADIC_CHA*padWidth[suid]*10,-0.5*padWidth[suid],(NUM_SPADIC_CHA-0.5)*padWidth[suid],
00464                                           2*NUM_SPADIC_CHA*padWidth[suid2]*10,-(NUM_SPADIC_CHA-0.5)*padWidth[suid2],(NUM_SPADIC_CHA-0.5)*padWidth[suid2]);
00465         title.Form("reconstructed hit position [mm] (Susibo %i)",usedSuId[suid]);
00466         posCorr2D[suid][suid2]->SetXTitle(title);
00467         title.Form("distance of reconstructed hit position (Susibo %i - Sisibo %i) [mm]",usedSuId[suid],usedSuId[suid2]);
00468         posCorr2D[suid][suid2]->SetYTitle(title);
00469         posCorr2D[suid][suid2]->SetContour(99);
00470 
00471         title.Form("wposCorr2D_%i_%i",usedSuId[suid],usedSuId[suid2]);
00472         name.Form("weighted position correlation 2D (Sus%i,Sus%i)",usedSuId[suid],usedSuId[suid2]);
00473         wposCorr2D[suid][suid2] = new TH2I(title, name, 
00474                                            NUM_SPADIC_CHA*padWidth[suid]*10,-0.5*padWidth[suid],(NUM_SPADIC_CHA-0.5)*padWidth[suid],
00475                                            2*NUM_SPADIC_CHA*padWidth[suid2]*10,-(NUM_SPADIC_CHA-0.5)*padWidth[suid2],(NUM_SPADIC_CHA-0.5)*padWidth[suid2]);
00476         title.Form("reconstructed hit position [mm] (Susibo %i)",usedSuId[suid]);
00477         wposCorr2D[suid][suid2]->SetXTitle(title);
00478         title.Form("distance of reconstructed hit position (Susibo %i - Sisibo %i) [mm]",usedSuId[suid],usedSuId[suid2]);
00479         wposCorr2D[suid][suid2]->SetYTitle(title);
00480         wposCorr2D[suid][suid2]->SetContour(99);
00481 
00482         title.Form("residuals_%i_%i",usedSuId[suid],usedSuId[suid2]);
00483         name.Form("residuals (Sus%i,Sus%i)",usedSuId[suid],usedSuId[suid2]);
00484         residuals[suid][suid2] = new TH1I(title, name, NUM_SPADIC_CHA*padWidth[suid]*5, -NUM_SPADIC_CHA*padWidth[suid], NUM_SPADIC_CHA*padWidth[suid]);
00485         title.Form("distance of reconstructed hit position (Susibo %i - Sisibo %i) [mm]",usedSuId[suid],usedSuId[suid2]);
00486         residuals[suid][suid2]->SetXTitle(title);
00487         residuals[suid][suid2]->SetYTitle("#");
00488 
00489         title.Form("wresiduals_%i_%i",usedSuId[suid],usedSuId[suid2]);
00490         name.Form("weighted residuals (Sus%i,Sus%i)",usedSuId[suid],usedSuId[suid2]);
00491         wresiduals[suid][suid2] = new TH1I(title, name, NUM_SPADIC_CHA*padWidth[suid]*5, -NUM_SPADIC_CHA*padWidth[suid], NUM_SPADIC_CHA*padWidth[suid]);
00492         title.Form("distance of reconstructed hit position (Susibo %i - Sisibo %i) [mm]",usedSuId[suid],usedSuId[suid2]);
00493         wresiduals[suid][suid2]->SetXTitle(title);
00494         wresiduals[suid][suid2]->SetYTitle("#");
00495         wresiduals[suid][suid2]->SetLineColor(3);
00496       }
00497     }
00498     while (fileSet.good()) {
00499       TString filename;
00500       fileSet >> filename;
00501       printf("  %s\n",filename.Data());
00502       if (filename.Sizeof() > 5) {
00503         TString ifname = filename;
00504         TString filefolder = filename; // use file name from file set to make subfolder in output file
00505         if(!ifname.EndsWith(".root")) ifname.Append(".root");   
00506         ifname = folder + filename;
00507         if(!ifname.EndsWith(".root")) ifname.Append(".root");
00508         TFile inputFile(ifname,"READ");
00509 
00510    
00511         TTree* theTree=0;
00512         TKey* kee=0;
00513         TIter iter(inputFile.GetListOfKeys());
00514         while ( ( kee=dynamic_cast<TKey*>(iter()) ) !=0 ) {
00515           theTree = dynamic_cast<TTree*>(kee->ReadObj());
00516           if (theTree)
00517             break; // we take first Tree in file, disregarding its name...
00518         }
00519         if(theTree==0) {
00520           cout <<"Error: no Tree in file "<< ifname.Data() << endl;
00521           return;
00522         }
00523   
00524         outputFile.mkdir(filefolder);
00525 
00526         // Thresholds
00527         Double_t minAmplitude[8] =      {        15,        15,        15,        15,        15,        15,        15,       15   };
00528         Int_t HitTimeWindow[8][2] =     {   {13,20},   {13,20},   {13,20},   {13,20},   {13,20},   {13,20},   {13,20},   {13,20}  };
00529         Int_t IntegrationWindow[8][2] = {   {14,22},   {14,22},   {14,22},   {14,22},   {14,22},   {14,22},   {14,22},   {14,22}  };
00530 
00531         Double_t dxPosmm[usedSusibos] = {0.0};
00532         Double_t wdxPosmm[usedSusibos] = {0.0};
00533 
00534         TCernOct11UnpackEvent* evnt = new TCernOct11UnpackEvent;
00535         TGo4EventElement* theBase=evnt;
00536         evnt->synchronizeWithTree(theTree, &theBase);
00537         //evnt->synchronizeWithTree(theTree, &evnt);
00538 
00539         // with go4 4.4.x:
00540 
00541         // evnt->synchronizeWithTree(theTree);
00542         // theTree->SetBranchAddress(theTree->GetListOfBranches()->First()->GetName(), &evnt);
00543 
00544         Int_t entries = theTree->GetEntries();
00545         //Int_t allEvents =0;
00546         //if (entries > 100000)
00547         //entries = 100000;
00548         cout << "   " << entries << " Entries" << endl;
00549         allEvents += entries;
00550         if (debug) 
00551           entries *= 0.1;
00552 
00553         Int_t FFMEvents = 0;
00554         Int_t MSEvents = 0;
00555         Int_t pulseProbes = 0;
00556         /*
00557           Int_t nPion[usedSusibos][nMerge] = {{0}};
00558           Int_t nElectron[usedSusibos][nMerge] = {{0}};
00559           Int_t nGoodEvents[usedSusibos][nMerge] = {{0}};
00560         */
00561         TMbsCrateEvent* fCrateInputEvent;
00562         TSpadicEvent* SpadicInputEvent;
00563         TSpadicData* theSpadic;
00564         TSpadicData* theSpadicA;
00565         TSpadicData* theSpadicB;
00566         //MFSpadicEvent* theMFSpadic = new MFSpadicEvent();
00567         //Float_t Ch1(0), Ch2(0), Pb(0);
00568 
00569 
00570         Int_t c_x(1680*4), c_y(1050*4);
00571         TCanvas* cResiduals = new TCanvas("cResiduals","cResiduals",c_x,c_y);
00572         cResiduals->Divide(usedSusibos,usedSusibos);
00573         TCanvas* cPosCorr = new TCanvas("cPosCorr","cPosCorr",c_x,c_y);
00574         cPosCorr->Divide(usedSusibos,usedSusibos);
00575         TCanvas* cwPosCorr = new TCanvas("cwPosCorr","cwPosCorr",c_x,c_y);
00576         cwPosCorr->Divide(usedSusibos,usedSusibos);
00577         TCanvas* cPos = new TCanvas("cPos","cPos",c_x,c_y/8*nIteration);
00578         cPos->Divide(usedSusibos,nIteration*2);
00579         TCanvas* cwPos = new TCanvas("cwPos","cwPos",c_x,c_y/8*nIteration);
00580         cwPos->Divide(usedSusibos,nIteration*2);
00581         TLegend* leg[usedSusibos][usedSusibos];
00582         TLegend* trackLeg[usedSusibos][nIteration];
00583         TString pr;
00584         TCanvas* cShape = new TCanvas("cShape","cShape",1600,600);
00585         cShape->Divide(2,1);
00586 
00587         //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00588         //Int_t rs(0), imerge(0), ifilter(0);
00589 
00590 
00591         for (Int_t rs = 0; rs < recoSets; rs++) {
00592           for (Int_t imerge = 0; imerge < nMerge; imerge++) {
00593             for (Int_t ifilter = 0; ifilter < nFilter; ifilter++) {
00594               printf("\n\n                ########   %s   %s   %s   ########\n\n",recoSetName[rs].Data(),mergeName[imerge].Data(),filterName[ifilter].Data());
00595               /*
00596                 appending.Form("_%s_%s_%s",recoSetName[rs].Data(),mergeName[imerge].Data(),filterName[ifilter].Data());
00597                 prResult[rs][imerge][ifilter]->Reset();
00598                 prResult[rs][imerge][ifilter]->SetTitle("prResult" + appending);
00599                 prwResult[rs][imerge][ifilter]->Reset();
00600                 prwResult[rs][imerge][ifilter]->SetTitle("prwResult" + appending);
00601               */
00602               for (Int_t iteration = 0; iteration < nIteration; iteration++) {
00603                 appending.Form("_%s_%s_%s_%i",recoSetName[rs].Data(),mergeName[imerge].Data(),filterName[ifilter].Data(),iteration);
00604                 hmTrack[iteration]->Reset();
00605                 hmTrack[iteration]->SetNameTitle("hmTrack" + appending,"hmTrack" + appending);
00606                 hbTrack[iteration]->Reset();
00607                 hbTrack[iteration]->SetNameTitle("hbTrack" + appending,"hbTrack" + appending);
00608               }
00609               for (Int_t suid = 0; suid < usedSusibos; suid++){
00610                 appending.Form("_%s_%i_%s_%s",recoSetName[rs].Data(),usedSuId[suid],mergeName[imerge].Data(),filterName[ifilter].Data());
00611                 prf[suid]->Reset();
00612                 prf[suid]->SetNameTitle("prf" + appending,"prf" + appending);
00613                 beamProfile[suid]->Reset();
00614                 beamProfile[suid]->SetNameTitle("beamProfile" + appending,"beamProfile" + appending);
00615                 wbeamProfile[suid]->Reset();
00616                 wbeamProfile[suid]->SetNameTitle("wbeamProfile" + appending,"wbeamProfile" + appending);
00617         
00618 
00619 
00620                 for (Int_t iteration = 0; iteration < nIteration; iteration++) {
00621                   appending.Form("_%s_%i_%s_%s%i",recoSetName[rs].Data(),usedSuId[suid],mergeName[imerge].Data(),filterName[ifilter].Data(),iteration);
00622                   trackResiduals[suid][iteration]->Reset();
00623                   trackResiduals[suid][iteration]->SetNameTitle("trackResiduals" + appending,"trackResiduals" + appending);
00624                   wtrackResiduals[suid][iteration]->Reset();
00625                   wtrackResiduals[suid][iteration]->SetNameTitle("wtrackResiduals" + appending,"wtrackResiduals" + appending);
00626                   trackPosCorr[suid][iteration]->Reset();
00627                   trackPosCorr[suid][iteration]->SetNameTitle("trackPosCorr" + appending,"trackPosCorr" + appending);
00628                   trackwPosCorr[suid][iteration]->Reset();
00629                   trackwPosCorr[suid][iteration]->SetNameTitle("trackwPosCorr" + appending,"trackwPosCorr" + appending);
00630                   trackPosCorr2D[suid][iteration]->Reset();
00631                   trackPosCorr2D[suid][iteration]->SetNameTitle("trackPosCorr2D" + appending,"trackPosCorr2D" + appending);
00632                   trackwPosCorr2D[suid][iteration]->Reset();
00633                   trackwPosCorr2D[suid][iteration]->SetNameTitle("trackwPosCorr2D" + appending,"trackwPosCorr2D" + appending);
00634                   //cout <<  trackwPosCorr2D[suid][iteration]->GetTitle() << endl;
00635                 }
00636                 for (Int_t suid2 = 0; suid2 < usedSusibos; suid2++){
00637                   appending.Form("_%s_%i_%i_%s_%s",recoSetName[rs].Data(),usedSuId[suid],usedSuId[suid2],mergeName[imerge].Data(),filterName[ifilter].Data());
00638                   chamberAlignment[suid][suid2]->Reset();
00639                   chamberAlignment[suid][suid2]->SetNameTitle("chamberAlignment" + appending,"chamberAlignment" + appending);
00640                   posCorr[suid][suid2]->Reset();
00641                   posCorr[suid][suid2]->SetNameTitle("posCorr" + appending,"posCorr" + appending);
00642                   wposCorr[suid][suid2]->Reset();
00643                   wposCorr[suid][suid2]->SetNameTitle("wposCorr" + appending,"wposCorr" + appending);
00644                   posCorr2D[suid][suid2]->Reset();
00645                   posCorr2D[suid][suid2]->SetNameTitle("posCorr2D" + appending,"posCorr2D" + appending);
00646                   wposCorr2D[suid][suid2]->Reset();
00647                   wposCorr2D[suid][suid2]->SetNameTitle("wposCorr2D" + appending,"wposCorr2D" + appending);
00648                   residuals[suid][suid2]->Reset();
00649                   residuals[suid][suid2]->SetNameTitle("residuals" + appending,"residuals" + appending);
00650                   wresiduals[suid][suid2]->Reset();
00651                   wresiduals[suid][suid2]->SetNameTitle("wresiduals" + appending,"wresiduals" + appending);
00652                 }
00653               }
00654               appending.Form("_%s_%s_%s",recoSetName[rs].Data(),mergeName[imerge].Data(),filterName[ifilter].Data());
00655               //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00656               //evnt->synchronizeWithTree(theTree, &theBase);
00657               for (Int_t iteration = 0; iteration < nIteration; iteration++) {
00658                 FFMEvents = 0;
00659                 MSEvents = 0;
00660                 printf("\n   Interation %i\n",iteration);
00661                 for(Int_t i = 0; i < entries; ++i) {
00662                   Statusbar(i,entries);
00663                   theTree->GetEntry(i);
00664                   fCrateInputEvent=dynamic_cast<TMbsCrateEvent*>(evnt->GetSubEvent("MBSCRATE"));    
00665                   SpadicInputEvent=dynamic_cast<TSpadicEvent*>(evnt->GetSubEvent("SPADIC"));
00666                   //int sid=2; // example: get spadic of id 1
00667                   for (Int_t suid = 0; suid < usedSusibos; suid++){
00668    
00669                     theSpadic=dynamic_cast<TSpadicData*>(SpadicInputEvent->getEventElement(usedSuId[suid]));
00670                     Double_t pulse[NUM_SPADIC_CHA][SPADIC_TRACE_SIZE] = {{0.0}};
00671             
00672                     Bool_t overFlow = theSpadic->fSpadicOverFlowEvent;//false;
00673                     Bool_t underShoot = theSpadic->fSpadicUndershootsEvent;
00674                     dxPosmm[suid] = 0.0;
00675                     wdxPosmm[suid] = 0.0;
00676 
00677                     //for (Int_t imerge = 0; imerge < nMerge; imerge++) {
00678                     msMerge = Bool_t(imerge);
00679                     if((suid < 4) && (msMerge)) {
00680                       theSpadicA = dynamic_cast<TSpadicData*>(SpadicInputEvent->getEventElement(usedSuId[suid]));
00681                       if(theSpadicA==0) continue; // skip not configured ones
00682                       if(!theSpadicA->IsValid()) continue; // required to suppress empty spadic data!
00683 
00684                       theSpadicB = dynamic_cast<TSpadicData*>(SpadicInputEvent->getEventElement(topSuId[suid]));
00685                       if(theSpadicB==0) continue; // skip not configured ones
00686                       if(!theSpadicB->IsValid()) continue; // required to suppress empty spadic data!
00687 
00688                       if (theSpadicA->fSpadicOverFlowEvent || theSpadicB->fSpadicOverFlowEvent) overFlow = true;
00689                       if (theSpadicA->fSpadicUndershootsEvent || theSpadicB->fSpadicUndershootsEvent) underShoot = true;
00690                       //MergeData(pulse, theSpadicA, theSpadicB);
00691                       Double_t pulseA[NUM_SPADIC_CHA][SPADIC_TRACE_SIZE] = {{0.0}};
00692                       Double_t pulseB[NUM_SPADIC_CHA][SPADIC_TRACE_SIZE] = {{0.0}};
00693                       for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {   
00694                         for (Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++) {   
00695                           pulseA[ch][bin] = theSpadicA->fSpadicCompensated[ch][bin];
00696                           pulseB[ch][bin] = theSpadicB->fSpadicCompensated[ch][bin];
00697                         }
00698                       }
00699                       BaselineCorrection(pulseA);
00700                       BaselineCorrection(pulseB);
00701                       for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {   
00702                         for (Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++) {
00703                           if (pulseA[ch][bin] > 0)
00704                             pulse[ch][bin] = pulseA[ch][bin];
00705                           if (pulseB[7-ch][bin] > 0) // channel remaping for MS top spadics
00706                             pulse[ch][bin] += pulseB[7-ch][bin];
00707                         }
00708                       }      
00709                     }
00710                     else{
00711                       //CopyData(pulse, theSpadic);
00712                       for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {   
00713                         for (Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++) {   
00714                           pulse[ch][bin] = theSpadic->fSpadicCompensated[ch][bin];
00715                      
00716                         }
00717                       }
00718                       BaselineCorrection(pulse);
00719                     }
00720                     if(theSpadic==0) continue; // skip not configured ones
00721                     if(!theSpadic->IsValid()) continue; // required to suppress empty spadic data!
00722 
00723                     //if (ifilter == 0)
00724 
00725                     if (ifilter == 1)
00726                       if (overFlow)
00727                         continue;
00728                     if (ifilter == 2)
00729                       if (overFlow || underShoot)
00730                         continue;
00731                     if (ifilter == 3)
00732                       if (underShoot)
00733                         continue;
00734 
00735                     // Find pad with maximum signal
00736                     Double_t maxAmplitude = 0.0;
00737                     Double_t maxAmplitudeCh[NUM_SPADIC_CHA] = {0.0};
00738                     Int_t maxCh = -1;
00739                     Int_t maxTB[NUM_SPADIC_CHA] = {-1};
00740                     for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {     
00741                       for (Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++) {
00742                         if (theSpadic->fSpadicOverflows[ch][bin])
00743                           overFlow = true;
00744                         if (pulse[ch][bin] > maxAmplitudeCh[ch]){
00745                           maxAmplitudeCh[ch] = pulse[ch][bin];
00746                           maxTB[ch] = bin;
00747                         }
00748                       }
00749                       if(maxAmplitudeCh[ch] > maxAmplitude){
00750                         maxAmplitude = maxAmplitudeCh[ch];
00751                         maxCh = ch;
00752                       }
00753                     }
00754                     if (maxCh < 1 || maxCh > 6)  continue; // skip fragmented signals
00755                     Double_t intensCh[NUM_SPADIC_CHA] = {0.0};
00756                     Bool_t drawProbe = false;
00757                     for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {
00758                       for (Int_t bin = 0; bin < SPADIC_TRACE_SIZE; bin++) {
00759                         if (debug)
00760                           if (pulseProbes < 20 && suid < 4){
00761                             if (maxTB[maxCh] > HitTimeWindow[suid][0] && maxTB[maxCh] < HitTimeWindow[suid][1]){
00762                               if (maxAmplitude > 150 && maxAmplitude < 175 ) {
00763                                 cPulseShape[ch]->SetBinContent(bin+1,pulse[ch][bin]);
00764                                 rPulseShape[ch]->SetBinContent(bin+1,theSpadic->fSpadicPulse[ch][bin]);
00765                                 if (bin == 0 && ch == 0) pulseProbes++;
00766                                 drawProbe = true;
00767                               }
00768                             }
00769                           }
00770 
00771                         if (rs == 0) {
00772                           if (pulse[ch][bin] > 0.0) {
00773                             intensCh[ch] +=  pulse[ch][bin];              
00774                           }
00775                         }
00776                         if (rs == 1) {
00777                           if (bin > IntegrationWindow[suid][0] && bin < IntegrationWindow[suid][1]) {
00778                             if (pulse[ch][bin] > 0.0) {
00779                               intensCh[ch] +=  pulse[ch][bin];            
00780                             }
00781                           }
00782                         }
00783                         if (rs == 2) {
00784                           intensCh[ch] = pulse[ch][maxTB[maxCh]];
00785                         }
00786                       }
00787                     }           
00788                     if (drawProbe) {    
00789                       for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {
00790                         cShape->cd(1);
00791                         rPulseShape[ch]->GetYaxis()->SetRangeUser(-10,255);
00792                         if (ch == 0)                           
00793                           rPulseShape[ch]->DrawCopy();
00794                         else
00795                           rPulseShape[ch]->DrawCopy("same");
00796                         cShape->cd(2);
00797                         cPulseShape[ch]->GetYaxis()->SetRangeUser(-10,255);
00798                         if (ch == 0)                           
00799                           cPulseShape[ch]->DrawCopy();
00800                         else
00801                           cPulseShape[ch]->DrawCopy("same");
00802                       }
00803                       TString nshape;
00804                       nshape.Form("%i",pulseProbes);
00805                       //cShape->SaveAs("data2011/Analysis/Pics/PR/Shape" + nshape + "_" + TString(nameFileSet(0,nameFileSet.Sizeof()-5)) + appending + ".pdf");
00806                       for (Int_t ch = 0; ch < NUM_SPADIC_CHA; ch++) {
00807                         cPulseShape[ch]->Reset();
00808                         rPulseShape[ch]->Reset();
00809                       }
00810                     }
00811                     if (maxTB[maxCh] > HitTimeWindow[suid][0] && maxTB[maxCh] < HitTimeWindow[suid][1]){
00812                       if (maxAmplitude > minAmplitude[suid]) {
00813                         if (maxCh >= 1 && maxCh <= NUM_SPADIC_CHA - 1){
00814                           if (intensCh[maxCh-1] > 0.0 && intensCh[maxCh] > 0.0 && intensCh[maxCh+1] > 0.0) {                
00815                             Double_t dxPos = 0.5 * log((intensCh[maxCh+1] / intensCh[maxCh-1])) / 
00816                               log((pow(intensCh[maxCh],2) / (intensCh[maxCh+1] * intensCh[maxCh-1])));
00817                             Double_t wdxPos = 1. / (pow(intensCh[maxCh-1],2) + pow(intensCh[maxCh+1],2)) * 
00818                               (
00819                                pow(intensCh[maxCh-1],2) * 
00820                                (pow(sigma[suid],2) / padWidth[suid] * 
00821                                 log(intensCh[maxCh] / intensCh[maxCh-1]) - 
00822                                 0.5 * padWidth[suid]) + 
00823                                pow(intensCh[maxCh+1],2) * 
00824                                (pow(sigma[suid],2) / padWidth[suid] * 
00825                                 log(intensCh[maxCh+1] / intensCh[maxCh]) + 
00826                                 0.5 * padWidth[suid])
00827                                );
00828                             dxPosmm[suid] = dxPos * padWidth[suid] + maxCh * padWidth[suid];
00829                             beamProfile[suid]->Fill(dxPosmm[suid]);
00830                             wdxPosmm[suid] = wdxPos  + maxCh * padWidth[suid];
00831                             wbeamProfile[suid]->Fill(wdxPosmm[suid]);
00832                           } //(intensCh[maxCh-1] > 0.0 && intensCh[maxCh] > 0.0 && intensCh[maxCh+1] > 0.0)
00833                           else continue;
00834                         } //(maxCh >= 1 && maxCh <= NUM_SPADIC_CHA - 1)
00835                         else continue;
00836                       } //(maxAmplitude > minAmplitude[suid])
00837                       else continue;
00838                     } //(maxTB[maxCh] > HitTimeWindow[suid][0] && maxTB[maxCh] < HitTimeWindow[suid][1])
00839                     else continue;
00840                     //} // for merge
00841                   } //for suid
00842 
00843                   //Alignment Iteration Compensation
00844                   Double_t  dxPosmmComp[usedSusibos] = {0.0};
00845                   Double_t wdxPosmmComp[usedSusibos] = {0.0};
00846 
00847                   for (Int_t suid = 0; suid < usedSusibos; suid++) {
00848                     dxPosmmComp[suid]  = dxPosmm[suid];
00849                     wdxPosmmComp[suid] = wdxPosmm[suid];
00850                     for (Int_t iIt = 0; iIt < iteration;iIt++) {
00851                       if (dxPosmm[suid] != 0.0) 
00852                         //dxPosmmComp[suid]  -= (mTrack[suid][iIt] *  dxPosmmComp[suid] + bTrack[suid][iIt]);
00853                         dxPosmmComp[suid]  -= (mTrack[suid][iIt] *  dxPosmm[suid] + bTrack[suid][iIt]);
00854                       if (wdxPosmm[suid] != 0.0) 
00855                         //wdxPosmmComp[suid] -= (mTrack[suid][iIt] * wdxPosmmComp[suid] + bTrack[suid][iIt]);
00856                         wdxPosmmComp[suid] -= (mTrack[suid][iIt] * wdxPosmm[suid] + bTrack[suid][iIt]);
00857                     }
00858                   }
00859 
00860                   Double_t distance[usedSusibos] = {0.0};
00861                   Double_t wdistance[usedSusibos] = {0.0};
00862                   Bool_t trackfound = Distance2TrackMS(dxPosmmComp, distance/*, hmTrack[iteration], hbTrack[iteration]*/);
00863                   Bool_t wtrackfound = Distance2TrackMS(wdxPosmmComp, wdistance/*, hmTrack[iteration], hbTrack[iteration]*/);
00864                   Bool_t trackfoundF = Distance2TrackFFM(dxPosmmComp, distance/*, m, b*/);
00865                   Bool_t wtrackfoundF = Distance2TrackFFM(wdxPosmmComp, wdistance/*, m, b*/);
00866 
00867 
00868                   if (!trackfound && !wtrackfound && !trackfoundF && !wtrackfoundF) continue;
00869 
00870                   for (Int_t suid = 0; suid < usedSusibos; suid++) {
00871                     if (dxPosmm[suid] != 0.0) {
00872                       if (suid < 4){
00873                         if (trackfound) {
00874                           MSEvents++;
00875                           trackPosCorr[suid][iteration]->Fill(dxPosmm[suid], distance[suid]);
00876                           trackPosCorr2D[suid][iteration]->Fill(dxPosmm[suid], distance[suid]);
00877                           trackResiduals[suid][iteration]->Fill(distance[suid]);
00878                         }
00879                         if (wtrackfound) {
00880                           trackwPosCorr[suid][iteration]->Fill(wdxPosmm[suid], wdistance[suid]);
00881                           trackwPosCorr2D[suid][iteration]->Fill(wdxPosmm[suid], wdistance[suid]);
00882                           wtrackResiduals[suid][iteration]->Fill(wdistance[suid]);
00883                         }
00884                       }
00885                       else {
00886                         if (trackfoundF) {
00887                           FFMEvents++;
00888                           trackPosCorr[suid][iteration]->Fill(dxPosmm[suid], distance[suid]);
00889                           trackPosCorr2D[suid][iteration]->Fill(dxPosmm[suid], distance[suid]);
00890                           trackResiduals[suid][iteration]->Fill(distance[suid]);
00891                         }
00892                         if (wtrackfoundF) {
00893                           trackwPosCorr[suid][iteration]->Fill(wdxPosmm[suid], wdistance[suid]);
00894                           trackwPosCorr2D[suid][iteration]->Fill(wdxPosmm[suid], wdistance[suid]);
00895                           wtrackResiduals[suid][iteration]->Fill(wdistance[suid]);
00896                         }
00897                       }
00898                     }
00899                     if (iteration == 0) {
00900                       for (Int_t suid2 = 0; suid2 < usedSusibos; suid2++) {
00901                         if (dxPosmm[suid] != 0.0 && dxPosmm[suid2] != 0.0) {
00902                           if (suid == suid2) {
00903 
00904                           }
00905                           else {
00906                             //trackResiduals[suid]->Fill(distance[suid] - (m[suid][suid2] * dxPosmm[suid] + b[suid][suid2]));
00907                             posCorr[suid][suid2]->Fill(dxPosmm[suid], dxPosmm[suid] - dxPosmm[suid2]);
00908                             posCorr2D[suid][suid2]->Fill(dxPosmm[suid], dxPosmm[suid] - dxPosmm[suid2]);              
00909                             residuals[suid][suid2]->Fill((dxPosmm[suid] - dxPosmm[suid2]) - (m[suid][suid2] * dxPosmm[suid] + b[suid][suid2]));
00910                           }             
00911                         }
00912                         if (wdxPosmm[suid] != 0.0 && wdxPosmm[suid2] != 0.0) {
00913                           if (suid == suid2) {
00914         
00915                           }
00916                           else {
00917                             //wtrackResiduals[suid]->Fill(wdistance[suid] - (m[suid][suid2] * dxPosmm[suid] + b[suid][suid2]));
00918                             wposCorr[suid][suid2]->Fill(wdxPosmm[suid], wdxPosmm[suid] - wdxPosmm[suid2]);
00919                             wposCorr2D[suid][suid2]->Fill(wdxPosmm[suid], wdxPosmm[suid] - wdxPosmm[suid2]);        
00920                             wresiduals[suid][suid2]->Fill((wdxPosmm[suid] - wdxPosmm[suid2]) - (m[suid][suid2] * wdxPosmm[suid] + b[suid][suid2]));
00921                           }
00922                         }
00923                       }
00924                     }
00925                   }
00926                 } // for entries
00927                 //Fit new slop and offset
00928                 printf("    Fraction of tracks with one event per chamber:\n     MS: %.2f%% FFM: %.2f%%\n",MSEvents*100/Float_t(entries),FFMEvents*100/Float_t(entries));
00929                 for (Int_t suid = 0; suid < usedSusibos; suid++) {
00930 
00931                   trackLeg[suid][iteration] = new TLegend(0.15,0.6,0.5,0.85);
00932                   trackLeg[suid][iteration]->SetLineColor(0);
00933                   trackLeg[suid][iteration]->SetLineStyle(0);
00934                   trackLeg[suid][iteration]->SetFillColor(0);
00935                   trackLeg[suid][iteration]->SetFillStyle(0);
00936                   trackLeg[suid][iteration]->SetTextSize(0.035);
00937 
00938                   // .........................................................................
00939                   cwPos->cd(1+suid+iteration*usedSusibos*2);
00940                   alignmentFit->SetParameter(0,0);
00941                   alignmentFit->SetParameter(1,0);
00942                   printf("Susibo %2i:",usedSuId[suid]);
00943                   trackwPosCorr[suid][iteration]->Fit("alignmentFit","Q","0",0.25*NUM_SPADIC_CHA*padWidth[suid],0.75*NUM_SPADIC_CHA*padWidth[suid]);
00944                   printf("    m = %8.4f    b = %8.4f   |   ", alignmentFit->GetParameter(0), alignmentFit->GetParameter(1));
00945                   if (iteration > 0) {
00946                     wtrackResiduals[suid][iteration]->Fit("prFit","Q0");
00947                     printf("PR = %6.3f+-%6.3fmm    ||    ",prFit->GetParameter(2) / sqrt(2.),prFit->GetParError(2) / sqrt(2.));
00948                   }
00949                   trackwPosCorr2D[suid][iteration]->DrawCopy("colz");
00950                   trackwPosCorr[suid][iteration]->GetYaxis()->SetRangeUser(-(NUM_SPADIC_CHA-0.5)*padWidth[suid],(NUM_SPADIC_CHA-0.5)*padWidth[suid]);
00951                   trackwPosCorr[suid][iteration]->DrawCopy("same");
00952                   //trackPosCorr[suid][iteration]->Fit("alignmentFit","Q","0",0.25*NUM_SPADIC_CHA*padWidth[suid],0.75*NUM_SPADIC_CHA*padWidth[suid]);
00953 
00954 
00955                   // .........................................................................
00956                   cwPos->cd(1+suid+(iteration*usedSusibos*2+usedSusibos));
00957                   wtrackResiduals[suid][iteration]->DrawCopy();
00958 
00959                   prFit->SetParameter(0,0);
00960                   prFit->SetParameter(1,0);
00961                   prFit->SetParameter(2,0);
00962 
00963                   if (iteration == 0)
00964                     wtrackResiduals[suid][iteration]->Fit("prFit","Q0");
00965                   else
00966                     wtrackResiduals[suid][iteration]->Fit("prFit","RQ0");
00967                   pr.Form("weighted PR = %.3fmm",prFit->GetParameter(2) / sqrt(2.));
00968                   if (iteration == nIteration-1)
00969                     prwResult[rs][imerge][ifilter]->Fill(suid,prFit->GetParameter(2) / sqrt(2.)/*, 1./pow(prFit->GetParError(2) / sqrt(2.),2)*/);
00970                   prFit->DrawCopy("same");
00971                   trackLeg[suid][iteration]->AddEntry((TObject*)0,pr,"");
00972                   trackLeg[suid][iteration]->Draw("same");
00973             
00974         
00975                   // .........................................................................
00976                   cPos->cd(1+suid+iteration*usedSusibos*2);
00977                   alignmentFit->SetParameter(0,0);
00978                   alignmentFit->SetParameter(1,0);
00979                   trackPosCorr[suid][iteration]->Fit("alignmentFit","Q","0",0.25*NUM_SPADIC_CHA*padWidth[suid],0.75*NUM_SPADIC_CHA*padWidth[suid]);
00980                   mTrack[suid][iteration] = alignmentFit->GetParameter(0);
00981                   bTrack[suid][iteration] = alignmentFit->GetParameter(1);
00982                   printf("m = %8.4f    b = %8.4f   |    ", mTrack[suid][iteration], bTrack[suid][iteration]);
00983                   if (iteration > 0) {
00984                     trackResiduals[suid][iteration]->Fit("prFit","Q0");
00985                     printf("PR = %6.3f+-%6.3fmm",prFit->GetParameter(2) / sqrt(2.),prFit->GetParError(2) / sqrt(2.));
00986                   }
00987                   printf("\n");
00988                   trackPosCorr2D[suid][iteration]->DrawCopy("colz");
00989                   trackPosCorr[suid][iteration]->GetYaxis()->SetRangeUser(-(NUM_SPADIC_CHA-0.5)*padWidth[suid],(NUM_SPADIC_CHA-0.5)*padWidth[suid]);
00990                   trackPosCorr[suid][iteration]->DrawCopy("same");
00991                   //trackPosCorr[suid][iteration]->Fit("alignmentFit","Q","0",0.25*NUM_SPADIC_CHA*padWidth[suid],0.75*NUM_SPADIC_CHA*padWidth[suid]);
00992  
00993                    
00994 
00995                   // .........................................................................
00996                   cPos->cd(1+suid+(iteration*usedSusibos*2+usedSusibos));
00997                   trackResiduals[suid][iteration]->DrawCopy();
00998             
00999                   prFit->SetParameter(0,0);
01000                   prFit->SetParameter(1,0);
01001                   prFit->SetParameter(2,0);
01002             
01003                   if (iteration == 0)
01004                     trackResiduals[suid][iteration]->Fit("prFit","Q0");
01005                   else
01006                     trackResiduals[suid][iteration]->Fit("prFit","RQ0");
01007                   pr.Form("unweighted PR = %.3fmm",prFit->GetParameter(2) / sqrt(2.));
01008                   if (iteration == nIteration-1)
01009                     prResult[rs][imerge][ifilter]->Fill(suid,prFit->GetParameter(2) / sqrt(2.)/*, 1./pow(prFit->GetParError(2) / sqrt(2.),2)*/);
01010                   prFit->DrawCopy("same");
01011                   trackLeg[suid][iteration]->AddEntry((TObject*)0,pr,"");
01012                   trackLeg[suid][iteration]->Draw("same");
01013 
01014                   // .........................................................................
01015                   outputFile.cd(filefolder);
01016                   wtrackResiduals[suid][iteration]->Write("", TObject::kOverwrite);
01017                   trackwPosCorr2D[suid][iteration]->Write("", TObject::kOverwrite);
01018                   trackwPosCorr[suid][iteration]->Write("", TObject::kOverwrite);
01019                   trackResiduals[suid][iteration]->Write("", TObject::kOverwrite);
01020                   trackPosCorr2D[suid][iteration]->Write("", TObject::kOverwrite);
01021                   trackPosCorr[suid][iteration]->Write("", TObject::kOverwrite);
01022                 } // suid
01023 
01024                 if (iteration == nIteration-1) {
01025                   //cPos->SaveAs("data2011/Analysis/Pics/PR/Pos_" + TString(nameFileSet(0,nameFileSet.Sizeof()-5)) + "_" + filefolder + appending + picEx);
01026                   //cwPos->SaveAs("data2011/Analysis/Pics/PR/wPos_" + TString(nameFileSet(0,nameFileSet.Sizeof()-5)) + "_" + filefolder + appending + picEx);
01027                   outputFile.cd(filefolder);
01028                   cPos->Write("", TObject::kOverwrite);
01029                   cwPos->Write("", TObject::kOverwrite);
01030                 }
01031          
01032               } // interation
01033               //inputFile.Close();
01034               for (Int_t iteration = 0; iteration < nIteration; iteration++){
01035                 for (Int_t suid = 0; suid < usedSusibos; suid++) {
01036                   wtrackResiduals[suid][iteration]->Reset();
01037                   trackwPosCorr2D[suid][iteration]->Reset();
01038                   trackwPosCorr[suid][iteration]->Reset();
01039                   trackResiduals[suid][iteration]->Reset();
01040                   trackPosCorr2D[suid][iteration]->Reset();
01041                   trackPosCorr[suid][iteration]->Reset();
01042                 }
01043               }
01044 
01045               for (Int_t suid = 0; suid < usedSusibos; suid++) {
01046                 for (Int_t suid2 = 0; suid2 < usedSusibos; suid2++) {
01047                   cPosCorr->cd(suid2*usedSusibos+suid+1);
01048                   posCorr2D[suid][suid2]->DrawCopy("colz");
01049                   //posCorr[suid][suid2]->Fit("alignmentFit","Q","0",0.25*NUM_SPADIC_CHA*padWidth[suid],0.75*NUM_SPADIC_CHA*padWidth[suid]);
01050                   //alignmentFit->DrawCopy("same");
01051                   posCorr[suid][suid2]->GetYaxis()->SetRangeUser(-(NUM_SPADIC_CHA-0.5)*padWidth[suid2],(NUM_SPADIC_CHA-0.5)*padWidth[suid2]);
01052                   posCorr[suid][suid2]->DrawCopy("same");
01053 
01054                   outputFile.cd(filefolder);
01055                   posCorr[suid][suid2]->Write("", TObject::kOverwrite);
01056                   posCorr2D[suid][suid2]->Write("", TObject::kOverwrite);
01057                   posCorr[suid][suid2]->Reset();
01058                   posCorr2D[suid][suid2]->Reset();
01059 
01060                   cwPosCorr->cd(suid2*usedSusibos+suid+1);
01061                   wposCorr2D[suid][suid2]->DrawCopy("colz");
01062                   //wposCorr[suid][suid2]->Fit("alignmentFit","Q","0",0.25*NUM_SPADIC_CHA*padWidth[suid],0.75*NUM_SPADIC_CHA*padWidth[suid]);
01063                   //alignmentFit->DrawCopy("same");
01064                   wposCorr[suid][suid2]->GetYaxis()->SetRangeUser(-(NUM_SPADIC_CHA-0.5)*padWidth[suid2],(NUM_SPADIC_CHA-0.5)*padWidth[suid2]);
01065                   wposCorr[suid][suid2]->DrawCopy("same");
01066 
01067                   outputFile.cd(filefolder);
01068                   wposCorr[suid][suid2]->Write("", TObject::kOverwrite);
01069                   wposCorr2D[suid][suid2]->Write("", TObject::kOverwrite);
01070                   wposCorr[suid][suid2]->Reset();
01071                   wposCorr2D[suid][suid2]->Reset();
01072 
01073                   leg[suid][suid2] = new TLegend(0.15,0.6,0.5,0.85);
01074                   leg[suid][suid2]->SetLineColor(0);
01075                   leg[suid][suid2]->SetLineStyle(0);
01076                   leg[suid][suid2]->SetFillColor(0);
01077                   leg[suid][suid2]->SetFillStyle(0);
01078                   leg[suid][suid2]->SetTextSize(0.035);
01079                   cResiduals->cd(suid2*usedSusibos+suid+1);
01080                   if (suid == suid2) {
01081                     /*
01082                       trackResiduals[suid]->DrawCopy();
01083                       trackResiduals[suid]->Fit("prFit","RQ0");
01084                       pr.Form("PR = %.3fmm",prFit->GetParameter(2) / sqrt(2.));
01085                       leg[suid][suid2]->AddEntry((TObject*)0,pr,"");
01086                       prFit->DrawCopy("same");
01087                       trackResiduals[suid]->DrawCopy("same");
01088                       wtrackResiduals[suid]->DrawCopy("same");
01089                       leg[suid][suid2]->Draw("same");
01090                     */
01091                   }
01092                   else {
01093                     residuals[suid][suid2]->DrawCopy();
01094                     residuals[suid][suid2]->Fit("prFit","RQ0");
01095                     pr.Form("PR = %.3fmm",prFit->GetParameter(2) / sqrt(2.));
01096                     leg[suid][suid2]->AddEntry((TObject*)0,pr,"");
01097                     prFit->DrawCopy("same");
01098                     residuals[suid][suid2]->DrawCopy("same");
01099                     wresiduals[suid][suid2]->DrawCopy("same");
01100                     leg[suid][suid2]->Draw("same");
01101                   }
01102                   outputFile.cd(filefolder);
01103                   residuals[suid][suid2]->Write("", TObject::kOverwrite);
01104                   wresiduals[suid][suid2]->Write("", TObject::kOverwrite);
01105                   residuals[suid][suid2]->Reset();
01106                   wresiduals[suid][suid2]->Reset();
01107                 }
01108                 outputFile.cd(filefolder);
01109                 beamProfile[suid]->Write("", TObject::kOverwrite);
01110                 wbeamProfile[suid]->Write("", TObject::kOverwrite);
01111                 beamProfile[suid]->Reset();
01112                 wbeamProfile[suid]->Reset();
01113                 //trackResiduals[suid]->Write("", TObject::kOverwrite);
01114                 //wtrackResiduals[suid]->Write("", TObject::kOverwrite);
01115               }
01116 
01117               for (Int_t iteration = 0; iteration < nIteration; iteration++){ 
01118                 outputFile.cd(filefolder);
01119                 hmTrack[iteration]->Write("", TObject::kOverwrite);
01120                 hbTrack[iteration]->Write("", TObject::kOverwrite);
01121                 hmTrack[iteration]->Reset();
01122                 hbTrack[iteration]->Reset();
01123               }
01124 
01125               outputFile.cd(filefolder);
01126               cwPosCorr->Write("", TObject::kOverwrite);
01127               cPosCorr->Write("", TObject::kOverwrite);
01128               cResiduals->Write("", TObject::kOverwrite);
01129 
01130               //cwPosCorr->SaveAs("data2011/Analysis/Pics/PR/wPosCorr_" + TString(nameFileSet(0,nameFileSet.Sizeof()-5)) + "_" + filefolder + appending + picEx);
01131               //cPosCorr->SaveAs("data2011/Analysis/Pics/PR/PosCorr_" + TString(nameFileSet(0,nameFileSet.Sizeof()-5)) + "_" + filefolder + appending + picEx);
01132               //cResiduals->SaveAs("data2011/Analysis/Pics/PR/PR_" + TString(nameFileSet(0,nameFileSet.Sizeof()-5)) + "_" + filefolder + appending + picEx);
01133 
01134               outputFile.cd();
01135 
01136 
01137             }
01138           }
01139         }
01140         cwPosCorr->Clear("D");
01141         cPosCorr->Clear("D");
01142         cResiduals->Clear("D");
01143         cPos->Clear("D");
01144         cwPos->Clear("D");
01145         /*
01146           delete cwPosCorr;
01147           delete cPosCorr;
01148           delete cResiduals;
01149           delete cPos;
01150           delete cwPos;
01151         */
01152         inputFile.Close();
01153       } // if filename is ok
01154     } // if fileset.open()
01155     TCanvas* cPrResults = new TCanvas("cPrResults","cPrResults",1600,600);
01156     for (Int_t rs = 0; rs < recoSets; rs++) {
01157       for (Int_t imerge = 0; imerge < nMerge; imerge++) {
01158         for (Int_t ifilter = 0; ifilter < nFilter; ifilter++) {   
01159           cPrResults->Divide(2,1);
01160           cPrResults->cd(1);
01161           prResult[rs][imerge][ifilter]->GetYaxis()->SetRangeUser(0,0.5);
01162           prResult[rs][imerge][ifilter]->DrawCopy("P");
01163           cPrResults->cd(2);
01164           prwResult[rs][imerge][ifilter]->GetYaxis()->SetRangeUser(0,0.5);
01165           prwResult[rs][imerge][ifilter]->DrawCopy("P");
01166           outputFile.cd();
01167           prResult[rs][imerge][ifilter]->Write("", TObject::kOverwrite);
01168           prwResult[rs][imerge][ifilter]->Write("", TObject::kOverwrite);
01169           cPrResults->Write("", TObject::kOverwrite);
01170           cPrResults->SaveAs("data2011/Analysis/Pics/PR/PR_Results_" + TString(nameFileSet(0,nameFileSet.Sizeof()-5)) + appending + ".pdf");
01171           cPrResults->Clear("D");
01172         }
01173       }
01174     }
01175     outputFile.Close();
01176     inputFileFit.Close();
01177   } // if fileSet.open()
01178 }