#include "TCTScan.h"

ClassImp(TCTScan)

TCTScan::TCTScan(Int_t num)
{
  // Constructor 
  //    Int_t num ; number of MeasureWF to be stored in the array
Num=num;
//filenames = new TClonesArray("TString",Num);
vf =new TClonesArray("MeasureWF",Num);
//vf->BypassStreamer(kFALSE);
}


TCTScan::~TCTScan()
{
  // Destructor
  // if(pt!=NULL) delete pt;
 delete vf;
//delete filenames;
}

void TCTScan::SetNum(Int_t num)
{
  // Constructor 
  //    Int_t num ; number of MeasureWF to be stored in the array
Num=num;
//filenames = new TClonesArray("TString",Num);
 if(vf!=NULL) delete vf;
vf =new TClonesArray("MeasureWF",Num);
//vf->BypassStreamer(kFALSE);
}

void TCTScan::AddMeasureWF(Int_t index,Char_t *name,Float_t offset,Int_t bins)
{
  //	Add new MeasureWF
  //		Int_t index    ;  index 
  //		Char_t *name   ;  file name of the measurement
  //		Float_t offset ;  start time shift
  TClonesArray &entryp = *vf;
  new(entryp[index]) MeasureWF(name,bins+1,offset);
}

void TCTScan::AddMeasureWF(Int_t index, MeasureWF *mwf)
{
  //	Add new MeasureWF
  //		Int_t index    ;  index 
  //		MeasuredWF *mwf;  measured vaweform
  TClonesArray &entryp = *vf;
  new(entryp[index]) MeasureWF(mwf->GetEntries());
   entryp[index]=mwf;			       
}

void TCTScan::Merge(TCTScan *scan)
{
  //	Merges two TCTScans!
Int_t newNum=scan->GetNum()+Num;
vf->ExpandCreate(newNum);
TClonesArray &entryp = *vf;
printf("Expanding finishedn");
for(Int_t i=Num;i<newNum;i++)
  entryp[i]=scan->Get(i-Num);
printf("Looping finishedn");
Num=newNum;
  
}

void TCTScan::Info(Int_t start,Int_t end)
{
  //Shows Info
if(start==-1111) start=0;
if(end==-1111) end=Num-1;
for(Int_t i=start;i<=end;i++)
  {
    printf(">>>>>>>> ENTRY NUMBER = %d !n",i);
    Get(i)->Info();
  printf("n");
  }
}
void TCTScan::NormArray(Int_t num,Float_t *array)
{
  //array normalization
Float_t max=0;
Int_t maxi=0,i=0;
for(i=0;i<num;i++) if(array[i]>max) {max=array[i]; maxi=i;}
for(i=0;i<num;i++) array[i]/=max;
}

TGraphErrors *TCTScan::FDV(Float_t mint, Float_t maxt,Float_t x0,Float_t x1,Float_t x2, Float_t x3,Int_t model,Int_t opt,Float_t StartTime, Float_t fluence,Float_t Errorx, Float_t Errory)
{
  //Draws FDV vs. time or vs. tempeature plot! CCE is drawn and fitted!
  //	     Float mint ; Time window 
  //	     Float maxt ; 
  //	     Float_t x0 ;  first point of the first line
  //	     Float_t x1 ;  second point of the first line
  //	     Float_t x2 ;  first point of the second line
  //	     Float_t x3 ;  second point of the second line
  //         Int_t model;  lin-lin=0  scale  sqrt-lin=1
  //	     Int_t opt  ;  opt=0; time  ,  opt=1; temperature	opt=2; Neff time; opt=3; Neff temperature
  //         Int_t StartTime; Start Time;
  //	     Float_t Fluence; fluence (10^12);
  //	     Float_t Errorx ; 
  //	     Float_t Errory ;  
Float_t *fdv=new Float_t[Num];
Float_t *x=new Float_t[Num];
Float_t *xe=new Float_t[Num];
Float_t *fdve=new Float_t[Num];
Int_t k=0;
if(StartTime==-1111) StartTime=Get(0)->GetTime();
for(Int_t i=0;i<Num;i++)
  {
    //    if(TMath::Abs(Get(i)->GetCurrent(Get(i)->GetEntries()-1))!=TMath::Abs(Get(i)->GetCurrent(Get(i)->GetEntries()-2)))
    //      {
    fdv[i]=Get(i)->GetFDV(mint,maxt,x0,x1,x2,x3,model,0); 
    if(opt>1) {fdv[i]=fdv[i]/(69.3*fluence)*100; fdve[i]=Errory;}    
    if(opt==0 || opt==2)  x[i]=Get(i)->GetTime()-StartTime; else x[i]=Get(i)->GetT();
    xe[i]=Errorx;
    //printf("%d %e %e n",i,fdv[k],Get(i)->GetCurrent(Get(i)->GetEntries()-1));
    //    printf("time %f ,%f mint=%f,maxt=%f,x0=%f,x1=%f,x2=%f,x3=%f n",x[i],fdv[i],mint, maxt, x0 ,x1,x2,x3);
    //    k++;
    //     }
  }
TGraphErrors *it=new TGraphErrors(Num,x,fdv,xe,fdve);
 it->SetLineWidth(4);
 it->SetMarkerStyle(21);
if(opt==0) it->SetTitle("FDV vs. Time"); else  if(opt==1) it->SetTitle("FDV vs. Temperature");
else if(opt==2) it->SetTitle("Neff vs. Time");  if(opt==3) it->SetTitle("Neff vs. Temperature");

if(opt==0 || opt==2) it->Draw("APL"); else it->Draw("AP");
if(opt==0 || opt==2) it->GetHistogram()->SetXTitle("time[h]"); else it->GetHistogram()->SetXTitle("T[C]");
if(opt==0 || opt==1) it->GetHistogram()->SetYTitle("FDV [V]"); else it->GetHistogram()->SetYTitle("N_{eff}/#Phi_{eq} [10^{-2} cm^{-1}]");
it->GetHistogram()->Draw();
if(opt==0 || opt==2) it->Draw("APL"); else it->Draw("AP");
return(it);
}

TGraph *TCTScan::FDV(Elec *el,Float_t x0,Float_t x1,Float_t x2, Float_t x3,Int_t model,Int_t opt)
{
  //Draws FDV vs. time or vs. tempeature plot using electronics. CCE is drawn and fitted with shaped signal max!
  //	     Elec *el   ;  electronics
  //	     Float_t x0 ;  first point of the first line
  //	     Float_t x1 ;  second point of the first line
  //	     Float_t x2 ;  first point of the second line
  //	     Float_t x3 ;  second point of the second line
  //         Int_t model;  lin-lin=0  scale  sqrt-lin=1
  //	     Int_t opt  ;  opt=0; temperature  ,  opt=1; time	
Float_t *fdv=new Float_t[Num];
Float_t *x=new Float_t[Num];
Float_t StartTime=Get(0)->GetTime();

for(Int_t i=0;i<Num;i++)
  {
    fdv[i]=Get(i)->GetFDV(el,x0,x1,x2,x3,model,0);
    if(opt==0)  x[i]=Get(i)->GetTime()-StartTime; else x[i]=Get(i)->GetT();
    //    printf("time %f ,%f mint=%f,maxt=%f,x0=%f,x1=%f,x2=%f,x3=%f n",x[i],fdv[i],mint, maxt, x0 ,x1,x2,x3);
  }
TGraph *it=new TGraph(Num,x,fdv);
 it->SetLineWidth(4);
 it->SetMarkerStyle(21);
if(opt==0) it->SetTitle("FDV vs. Time"); else  it->SetTitle("FDV vs. Temperature");
if(opt==0) it->Draw("APL"); else it->Draw("AP");
if(opt==0) it->GetHistogram()->SetXTitle("time[h]"); else it->GetHistogram()->SetXTitle("T[C]");
 it->GetHistogram()->SetYTitle("FDv[V]");
it->GetHistogram()->Draw();
if(opt==0 ) it->Draw("APL"); else it->Draw("AP");
return(it);
}

TGraph *TCTScan::CCEE(Elec *el,Float_t Volt,Int_t opt)
{
  // Draws CCE vs. temperature or time at given voltage using electronics
  //		Elec *el     ; electronics
  //		Float_t Volt ; voltage
  //	        Int_t opt    ;  opt=0; temperature  ,  opt=1; time     
TH1F *his=new TH1F();
Float_t *I=new Float_t[Num];
Float_t *x=new Float_t[Num];
Double_t xmax;
Double_t xmin;
Float_t StartTime=Get(0)->GetTime();
for(Int_t i=0;i<Num;i++)
  {
    Get(i)->GetHistogram(Volt,his);
    xmax=his->GetXaxis()->GetXmax();
    xmin=his->GetXaxis()->GetXmin();
    his->GetXaxis()->SetLimits(xmin*1e-9,xmax*1e-9);
    el->preamp(his);
    el->CRshape(his);
    el->RCshape(his);
    el->RCshape(his);
    I[i]=his->GetMaximum();
    if(opt==0)  x[i]=Get(i)->GetTime()-StartTime; else x[i]=Get(i)->GetT(Volt);
  }

NormArray(Num,I);
TGraph *it=new TGraph(Num,x,I);
 it->SetLineWidth(4);
 it->SetMarkerStyle(21);
if(opt==0) it->SetTitle("CCE vs. Time"); else  it->SetTitle("CCE (electronics) vs. Temperature");
if(opt==0) it->Draw("APL"); else it->Draw("AP"); 
if(opt==0) it->GetHistogram()->SetXTitle("time[h]"); else it->GetHistogram()->SetXTitle("T[C]");
 it->GetHistogram()->SetYTitle("Charge [arb.]");
it->GetHistogram()->Draw();
if(opt==0) it->Draw("APL"); else it->Draw("AP"); 
return(it);
}


TGraph *TCTScan::CCE(Float_t Volt,Int_t opt,Float_t mint, Float_t maxt)
{
  // Draws CCE vs. temperature or time at given voltage using charge
  //		Float_t Volt ; voltage
  //	        Int_t opt    ;  opt=0; temperature  ,  opt=1; time 
  //		Float_t maxt ; integration interval 
  //		Float_t mint ; 
Float_t *I=new Float_t[Num];
Float_t *x=new Float_t[Num];
Float_t StartTime=Get(0)->GetTime();
for(Int_t i=0;i<Num;i++)
  {
    I[i]=Get(i)->Integral(Volt,mint,maxt);
    if(opt==0)  x[i]=Get(i)->GetTime()-StartTime; else x[i]=Get(i)->GetT(Volt);
  }


TGraph *it=new TGraph(Num,x,I);
 it->SetLineWidth(4);
 it->SetMarkerStyle(21);
if(opt==0) it->SetTitle("CCE vs. Time"); else  it->SetTitle("CCE vs. Temperature");
if(opt==0) it->Draw("APL"); else it->Draw("AP");
if(opt==0) it->GetHistogram()->SetXTitle("time[h]"); else it->GetHistogram()->SetXTitle("T[C]");
 it->GetHistogram()->SetYTitle("CCE");
it->GetHistogram()->Draw();
if(opt==0) it->Draw("APL"); else it->Draw("AP");

return(it);
}

void TCTScan::ReadFile(Char_t *name, Int_t *step,Int_t count,Float_t shift,Int_t bins)
{
  // Reads multiple files of type .tct
  //	Char_t *name ; file name -> the part of the string that changes should be replaced with @
  //	Float_t step ; step size
  //	Int_t count  ; index in the vector (TCTScan) to put the first measurement in
  //    Float_t shift; (same shift as in the constructor of the MeasureWF)
Int_t index,k=0,offset=0;
Int_t i;
Char_t filename[300];
Char_t extension[100];

Char_t v[6];
while(name[k]!='@' && name[k]!='0') k++; index=k; 
strcpy(extension,&name[index+1]); name[index]='0';


for(i=0;i<count;i++)
  {

    sprintf(filename,"%s%d%s",name,step[i],extension);
    printf("%d %sn",count,filename);
    AddMeasureWF(i,filename,shift,bins);
  }
}

Int_t TCTScan::i2a(Char_t *v,Float_t vol,Int_t sign)
{
  // int to ascii with signs (same as in MeasureWF)
Int_t k=0;
       if(sign) if(vol>0) v[k++]='+'; else {v[k++]='-'; vol=-vol;}
       if(!sign) vol=vol>0?vol:-vol; 
       if((Int_t) vol/100!=0) v[k++]=(Char_t)(((Int_t) vol)/100)+48;
       v[k++]=(Char_t)(((Int_t) vol%100)/10)+48;
       v[k++]=(Char_t)((Int_t) vol%10)+48; 
       v[k++]='0';
       return k-1;
}

void TCTScan::DrawMulti(Float_t volt,Float_t low,Float_t high,Int_t Start,Int_t End, Int_t Step,Int_t Deconv)
{
  // Draws Multiple vaweforms at given voltage vs. temperature or time
  //		Float_t volt ; voltage
  //		Float_t low  ; time interval to display
  //		Float_t high ;
  //		Int_t Start  ; start index in the vector
  // 		Int_t End    ; end index in the vector
  // 		Int_t Step   ; step (default=1) 
  //		Int_t Deconv ; Deconvolute=1, No Deconvolution=0
Char_t v[5];
Int_t cii=0;
Int_t colori[]={1,2,3,4,5,6,7,13,28,30,34,38,40,31,46,49,1,2,3,4,5,6,7};
Float_t max=0;
Int_t i,maxi,color; 
Elec *elTCT=new Elec(10e-12,50,100e-9);
TH1F *ch;
TH1F *his,*hisd,*hism;
if(Start==-1) Start=0;
if(End==-1) End=Num-1;

 for(i=Start;i<=End;i+=Step) 
   {
     ch=Get(i)->GetHA(volt); 
     if(ch->GetMaximum()>max) {max=ch->GetMaximum(); maxi=i;}
   }
 if(max==0) maxi=End; //printf("maxi=%d",maxi);

ch=Get(maxi)->GetHA(volt);   
elTCT->SetCp(Get(max)->CAP);
 

 if(Deconv)
   {
     hism=new TH1F("Deconvoluted","Deconvoluted",ch->GetNbinsX(),ch->GetXaxis()->GetXmin(),
     ch->GetXaxis()->GetXmax()); 
     hisd=new TH1F("Deconvolution","Deconvoution",ch->GetNbinsX(),ch->GetXaxis()->GetXmin(),
     ch->GetXaxis()->GetXmax()); 
     ch->Copy(*hism); elTCT->Revpreamp((TH1F *)hism,1e9);
   } 
   else hism=ch;
 
 hism->SetXTitle("t[ns]");
 hism->SetYTitle("I[V/50#Omega]");
 hism->SetLabelSize(0.045,"X");
 hism->SetLabelSize(0.045,"Y");
 i2a(v,volt,0); //title without signs
 TString title="TCT Measurement @ U=";title=title+v; title=title+" V";
 hism->SetTitle((const char *)title);

if(low!=-1111. || high!=-1111.) hism->GetXaxis()->SetRange( hism->GetXaxis()->FindBin(low),hism->GetXaxis()->FindBin(high));
 hism->DrawCopy();
 
 for(i=Start;i<=End;i+=Step)
    {
      color=colori[cii];
    //      color=i/7*40+i%7+1; 
      his=Get(i)->GetHA(volt);   
      if(Deconv) 
	 {
	  his->Copy(*hisd); 
	  elTCT->SetCp(Get(i)->CAP);
	  elTCT->Revpreamp((TH1F *)hisd,1e9);
	 } else hisd=his;
      //      color=colori[i-Start];
      hisd->SetLineColor((Color_t)color);
      hisd->DrawCopy("SAME"); 
      cii++;    
    }
    Legend(hism,Start,End,Step);
}

void TCTScan::DrawMultiCS(Float_t volt,Float_t low,Float_t high,Int_t Start,Int_t End)
{
  // Draws multi current shapes (see MeasureWF) at given voltage vs. temperature or time
  //		Float_t volt ; voltage
  //		Float_t low  ; time interval to display
  //		Float_t high ;
  //		Int_t Start  ; start index in the vector
  // 		Int_t End    ; end index in the vector
Char_t v[5];
Int_t colori[]={1,2,3,4,5,6,7,13,28,30,34,38,40,31,46,49,1,2,3,4,5,6,7};
Float_t max=0;
Int_t i,maxi,color; 
TH1F *ch;
TH1F *his;
if(Start==-1) Start=0;
if(End==-1) End=Num-1;

 for(i=Start;i<=End;i++) 
   {
     ch=Get(i)->ChargeShape(volt); 
     if(ch->GetMaximum()>max) {max=ch->GetMaximum(); maxi=i;}
   }
 if(max==0) maxi=End; //printf("maxi=%d",maxi);

ch=Get(maxi)->ChargeShape(volt);   

if(low!=-1111. || high!=-1111.) ch->GetXaxis()->SetRange( ch->GetXaxis()->FindBin(low),ch->GetXaxis()->FindBin(high));
  
 ch->SetXTitle("t[ns]");
 ch->SetYTitle("Charge [arb.]");
 ch->SetLabelSize(0.045,"X");
 ch->SetLabelSize(0.045,"Y");
 i2a(v,volt,0);
 TString title="TCT Measurement @ U=";title=title+v; title=title+" V";
 ch->SetTitle((const char *)title);
 ch->DrawCopy();

 for(i=Start;i<=End;i++)
    {
      //      color=i/7*40+i%7+1; 
      his=Get(i)->ChargeShape(volt);   
      color=colori[i-Start];
      his->SetLineColor((Color_t)color);
      his->DrawCopy("SAME");
      }
    Legend(ch,Start,End);
}

void TCTScan::DrawMultiCCE(Float_t minit,Float_t maxit,Int_t Start,Int_t End,Int_t model)
{
  // Draws multiple CCE (see MeasureWF) vs. temperature or time
  //		Float_t minit ; integration time interval
  //		Float_t maxit ;
  //		Int_t Start   ; start index in the vector
  // 		Int_t End     ; end index in the vector
  //		Int_t model   ; lin-lin=0  scale  sqrt-lin=1
Char_t v[5];
Int_t colori[]={1,2,3,4,5,6,7,13,28,30,34,38,40,31,46,49,1,2,3,4,5,6,7};
Float_t max=0;
Int_t i,maxi,color; 
TGraph *ch;
TH1F *his;
if(Start==-1) Start=0;
if(End==-1) End=Num-1;

 for(i=Start;i<=End;i++) 
   {
     ch=Get(i)->CCE(minit,maxit,model,0);
     color=colori[i-Start];
     ch->SetLineColor((Color_t)color);  
     ch->SetLineWidth(2);
     if(i==Start) 
       {
	 ch->Draw("AL"); 
	 ch->SetTitle("Charge vs. Voltage (diff. T)");
	 if(model!=1) ch->GetHistogram()->SetXTitle("U[V]"); else ch->GetHistogram()->SetXTitle("Sqrt U[ Sqrt V]");//	 ch->GetHistogram()->SetXTitle("U[V]");
	 ch->GetHistogram()->SetYTitle("Charge [arb.]");
	 ch->GetHistogram()->Draw(); 
	 ch->Draw("AL"); 
	 his=ch->GetHistogram();
       }
     else ch->Draw("L");	
   }

   Legend(his,Start,End);
}


TGraph *TCTScan::DrawIT(Float_t Volt,Int_t opt)
{
  // Draws current at given voltage vs. temperature (time)
Int_t i;
Float_t *x=new Float_t [Num];
Float_t *y=new Float_t [Num];

Float_t StartTime=Get(0)->GetTime();
for(i=0;i<Num;i++)
  {
    y[i]=Get(i)->GetCurrent(Volt)*1e6;    
    if(opt==0)  x[i]=Get(i)->GetTime()-StartTime; else x[i]=Get(i)->GetT();
  }
 TGraph *it=new TGraph(Num,x,y);

 it->SetLineWidth(4);
 it->SetMarkerStyle(21);
 if(opt==0) it->SetTitle("Current vs. Time"); else it->SetTitle("Current vs. Temperature");
if(opt==0) it->Draw("APL"); else it->Draw("AP"); 
if(opt==0) it->GetHistogram()->SetXTitle("time[h]"); else it->GetHistogram()->SetXTitle("T[C]");
 it->GetHistogram()->SetYTitle("I[uA]");
 it->GetHistogram()->Draw();
if(opt==0) it->Draw("APL"); else it->Draw("AP"); 
return(it);

}


void TCTScan::Legend(TH1F *ch,Int_t start, Int_t end,Int_t Step)
{
  //Draws legend
Float_t minx,miny,maxy,maxx,x1,x2,y1,y2;
TString title,utit="T=";
Char_t v[5]; 
Int_t color,cii=0;;
Int_t colori[]={1,2,3,4,5,6,7,13,28,30,34,38,40,31,46,49,1,2,3,4,5,6,7};
TText *text;

minx=ch->GetXaxis()->GetBinCenter(ch->GetXaxis()->GetFirst());
maxx=ch->GetXaxis()->GetBinCenter(ch->GetXaxis()->GetLast());
miny=ch->GetMinimum();
maxy=ch->GetMaximum();

x1=(maxx-minx)*0.6+minx; 
x2=(maxx-minx)*0.9+minx;
y2=(maxy-miny)*0.35+miny;
y1=(maxy-miny)*0.95+miny;

//printf("coords: x1=%f y1=%f x2=%f y2=%fn",x1,y1,x2,y2);

 if(pt!=NULL) delete pt;
 pt=new TPaveText(x1,y1,x2,y2);

for(Int_t i=start;i<=end;i+=Step)
  { 
   i2a(v,(Int_t) Get(i)->GetT()); 
   //   color=i/7*40+i%7+1; 
   color=colori[cii];
   title=utit+v; title=title+" C";
   text=pt->AddText((const char *)title);
   text->SetTextColor((Color_t)color);
   text->SetTextSize(0.05);
   cii++;
  }
pt->Draw();
}

/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
ClassImp(PSTCT)

 PSTCT::PSTCT(char *FileName, Float_t time0, Int_t Bin)
{
  //////////////////////////////////////////////////////////////////////////// 
  // Class for manipulation of the 3D/Position Sensitive TCT Measurements
  // 
  // char *FileName;  name of the file with data
  // Float_t time0 ;  time shift of the points. Used mainly to set the arrival of the laser pulse at t=0 ns;
  // Int Bin; Selects the timeformat of the measurements; 
  //          0 - ascii (default) older format has a type 11 while newer has type 22
  //          1 - binary (prefered in new measurements)
  // 
  // Example of use :
  // // Convert into MeasureWF along projection
  // PSTCT aa("../Meritve/scanz-grobo-1.tct", 92.2,1); // The second parameter is to set the scale such that signal start at t=0;
  // aa.CorrectBaseLine();   // Baseline correction 
  // aa.PrintInfo();         // Information about the read data 
  // // example of projection the data along Y, 41 point at indexes of Z=23,X=0, U1=0, U2=0;
  // MeasureWF *wf=aa.Projection(0,1,0,0,23,0,0,41); 
  // strcpy(wf.suffix," #mum");
  // strcpy(wf.prefix,"Y=");
  // //Draw waveforms and substract the 0 waveform from the rest (cancel oscilations
  // wf.DrawMulti(-1,40,1,40,4,0,0);



Int_t i,j;
Char_t filef[5];
 float header[100];
for(i=0;i<3;i++) WFOnOff[i]=1;
Date=TArrayI(6);

if(!Bin) sprintf(filef,"r+"); else sprintf(filef,"rb+");
if((in=fopen((const Char_t *)FileName,filef))==NULL) {printf("n Error opening file for readingn"); return;}

 if(!Bin)
   {
fscanf(in,"%d",&type); 
 if(!(type==11 || type==22)) {
   printf("Can not read other formats than waveform: %d!n",type);}
   else
     {
fscanf(in,"%d %d %d %d %d %dn",&Date[0],&Date[1],&Date[2],&Date[3],&Date[4],&Date[5]);
fscanf(in,"%dn",&abstime);
fscanf(in,"%f %f %dn",&x0,&dx,&Nx);
fscanf(in,"%f %f %dn",&y0,&dy,&Ny);
fscanf(in,"%f %f %dn",&z0,&dz,&Nz);
if(type==22) fscanf(in,"%d %d %dn",&WFOnOff[0],&WFOnOff[1],&WFOnOff[2]); //Read in the wafeform off on
fscanf(in,"%d",&NU1); U1=TArrayF(NU1); 
 for(i=0;i<NU1;i++) fscanf(in,"%f",&U1[i]); // printf("%d %d %fn",NU1,i,U1[i]);}

 if(type==22) fscanf(in,"%d",&NU2); else NU2=1;
 U2=TArrayF(NU2); 
 if(type==22) for(i=0;i<NU2;i++) fscanf(in,"%f",&U2[i]);

 I2=TArrayF(NU2*NU1); 
 I1=TArrayF(NU2*NU1); 

 fscanf(in,"%f %f %dn",&t0,&dt,&NP);

 numxyz=Nx*Ny*Nz;
 
 for(i=0;i<8;i++) xyz[i]=new Float_t [numxyz*NU1*NU2];

 if(WFOnOff[0]) {histo1 =new TClonesArray("TH1F",numxyz*NU1*NU2); histo1->BypassStreamer(kFALSE);}
 if(WFOnOff[1]) {histo2 =new TClonesArray("TH1F",numxyz*NU1*NU2); histo2->BypassStreamer(kFALSE);}
 if(WFOnOff[2]) {histo3 =new TClonesArray("TH1F",numxyz*NU1*NU2); histo3->BypassStreamer(kFALSE);}
     }
 ReadWFs(time0);
   }
 else
   {
     int read=fread((void *)header,sizeof(float),100,in);
     swooip(header,read);
     printf("read char=%dn",read); 
	 type=(Int_t) header[0]; //for simplicity reasons only - not used
	 for(i=0;i<6;i++) Date[i]=(Int_t) header[i+1];
     abstime=(int) header[7];
	 x0=header[8];  dx=header[9]; Nx=(int)header[10]; 
     y0=header[11]; dy=header[12]; Ny=(int)header[13]; 
     z0=header[14]; dz=header[15]; Nz=(int)header[16]; 
     for(i=0;i<3;i++) WFOnOff[i]=(int)header[17+i];
     NU1=(int)header[20];U1=TArrayF(NU1);
     for(i=0;i<NU1;i++) U1[i]=header[21+i];
     NU2=(int)header[21+NU1];U2=TArrayF(NU2);
     for(i=0;i<NU2;i++) U2[i]=header[22+NU1+i];
     t0=header[22+NU1+NU2];
     dt=header[23+NU1+NU2];
     NP=(int)header[24+NU1+NU2];

    I2=TArrayF(NU2*NU1); 
    I1=TArrayF(NU2*NU1); 
    numxyz=Nx*Ny*Nz;
 
    for(i=0;i<8;i++) xyz[i]=new Float_t [numxyz*NU1*NU2];
 if(WFOnOff[0]) {histo1 =new TClonesArray("TH1F",numxyz*NU1*NU2); histo1->BypassStreamer(kFALSE);}
 if(WFOnOff[1]) {histo2 =new TClonesArray("TH1F",numxyz*NU1*NU2); histo2->BypassStreamer(kFALSE);}
 if(WFOnOff[2]) {histo3 =new TClonesArray("TH1F",numxyz*NU1*NU2); histo3->BypassStreamer(kFALSE);}
 //for(i=0;i<50;i++) printf("%d %fn",i,header[i]);
 
 fseek(in,(25+NU1+NU2)*sizeof(Float_t),SEEK_SET);
 ReadWFsBin(time0);

   }

 RefInd=-1;
 //Setting the color map
}


 void  PSTCT::ReadWFsBin(Float_t time0)
{
  Int_t i,ii,j,k,q,r,numread;
  Float_t data,tU1,tU2,tI1,tI2;
  Char_t hisname1[100];
  Char_t hisname2[100];
  Char_t hisname3[100];
  TClonesArray &entryp1 = *histo1;  
  TClonesArray &entryp2 = *histo2;  
  TClonesArray &entryp3 = *histo3;  
  Float_t buf[10000];

  for(q=0;q<NU1;q++)
       {
     for(r=0;r<NU2;r++)
         {

	       fread((void *)buf,sizeof(Float_t),4,in); swooip(buf,4);
               tU1=buf[0]; tU2=buf[1]; tI1=buf[2]; tI2=buf[3];
	       U1[q]=tU1; I1[r+q*NU2]=tI1;
	       U2[r]=tU2; I2[r+q*NU2]=tI2;
	     
	       //	       printf("%d %d :: %f %f %f %fn",r,q,tU1,tU2,tI1,tI2);
    
	   for(i=0;i<numxyz;i++)
	     {

        ii=i+numxyz*r+(NU2*numxyz)*q;
	
               fread((void *)buf,sizeof(Float_t),4,in);  swooip(buf,4); 
	       //	        printf("%d :: %f %f %f %fn",ii,buf[0],buf[1],buf[2],buf[3]);
	        for(j=0;j<4;j++) xyz[j][ii]=buf[j]; 
		xyz[7][ii]=xyz[3][ii]; 
	        xyz[3][ii]=tU1; xyz[4][ii]=tU2; 
		xyz[5][ii]=tI1; xyz[6][ii]=tI2;
           

	if(WFOnOff[0])
	  {
	    sprintf(hisname1,"Ch. 1:x=%5.3e,y=%5.3e,z=%5.3e,U1=%4.2f, U2=%4.2f ",xyz[0][ii],xyz[1][ii],xyz[2][ii],xyz[3][ii],xyz[4][ii]);
	    new(entryp1[ii]) TH1F((const Char_t *)(hisname1),(const Char_t *)(hisname1),NP,t0*1e9-time0,(NP*dt+t0)*1e9-time0);
	  }
	if(WFOnOff[1])
	  {
	    sprintf(hisname2,"Ch. 2:x=%5.3e,y=%5.3e,z=%5.3e,U1=%4.2f, U2=%4.2f ",xyz[0][ii],xyz[1][ii],xyz[2][ii],xyz[3][ii],xyz[4][ii]);
            new(entryp2[ii]) TH1F((const Char_t *)(hisname2),(const Char_t *)(hisname2),NP,t0*1e9-time0,(NP*dt+t0)*1e9-time0);
	  }
	if(WFOnOff[2])
	  {
	    sprintf(hisname3,"Ch. 3:x=%5.3e,y=%5.3e,z=%5.3e,U1=%4.2f, U2=%4.2f ",xyz[0][ii],xyz[1][ii],xyz[2][ii],xyz[3][ii],xyz[4][ii]);
            new(entryp3[ii]) TH1F((const Char_t *)(hisname3),(const Char_t *)(hisname3),NP,t0*1e9-time0,(NP*dt+t0)*1e9-time0);
	  }


	for(k=0;k<3;k++)
	  {
	   if(WFOnOff[k]==1)
	      {
		//    printf("reading ...... %d %d ... ",k,ftell(in)/4);
		numread=fread(buf,sizeof(Float_t),NP,in);  swooip(buf, NP); // printf("%f %f %fn",buf[NP-3],buf[NP-2],buf[NP-1]);
		//		printf("read ...... %d(%d)n",numread,ii);  if(ii==17) return;
		//fread(buf,sizeof(Float_t),NP,in);  swooip(buf, NP);
		for(j=0;j<NP;j++)
		  {
		         //if(ii<2) if(j%20!=0) printf("%5.2f", buf[j]); else printf("%5.2fn", buf[j]);
		    switch(k)
		      {
			case 0: 
			  ((TH1F*)entryp1[ii])->SetBinContent(j+1,buf[j]);
			  break;
			case 1: 
			  ((TH1F*)entryp2[ii])->SetBinContent(j+1,buf[j]);
			  break;
			case 2: 
			  ((TH1F*)entryp3[ii])->SetBinContent(j+1,buf[j]);
			  break;

		      }
		  }
		
		//printf("k=%d n",k);
	
	      }
	
	  }
	
       }
	
	 }

       }

    //fclose(in); 
}



 void  PSTCT::ReadWFs(Float_t time0)
{
  Int_t i,ii,j,k,q,r;
  Float_t data,tU1,tU2,tI1,tI2;
  Char_t hisname1[100];
  Char_t hisname2[100];
  Char_t hisname3[100];
  TClonesArray &entryp1 = *histo1;  
  TClonesArray &entryp2 = *histo2;  
  TClonesArray &entryp3 = *histo3;  
   

  for(q=0;q<NU1;q++)
       {
     for(r=0;r<NU2;r++)
         {


           if(type==22) 
	     {
	       fscanf(in,"%f %f %f %f",&tU1,&tU2,&tI1,&tI2);	     
	       //	       printf("%d %d :: %f %f %f %fn",r,q,tU1,tU2,tI1,tI2);	     
		 U1[q]=tU1; I1[r+q*NU2]=tI1;
	         U2[r]=tU2; I2[r+q*NU2]=tI2;
	     }
	     
    
	   for(i=0;i<numxyz;i++)
	     {

        ii=i+numxyz*r+(NU2*numxyz)*q;
	
			if(type==22) 
			  {
					for(j=0;j<4;j++) fscanf(in,"%f",&xyz[j][ii]); xyz[7][ii]=xyz[3][ii]; 
							  xyz[3][ii]=tU1; xyz[4][ii]=tU2; xyz[5][ii]=tI1; xyz[6][ii]=tI2;
			  } 
			 if(type==11) for(j=0;j<5;j++) fscanf(in,"%f",&xyz[j][ii]);

			if(WFOnOff[0])
			  {
				sprintf(hisname1,"Ch. 1:x=%5.3e,y=%5.3e,z=%5.3e,U1=%4.2f, U2=%4.2f ",xyz[0][ii],xyz[1][ii],xyz[2][ii],xyz[3][ii],xyz[4][ii]);
				new(entryp1[ii]) TH1F((const Char_t *)(hisname1),(const Char_t *)(hisname1),NP,t0*1e9-time0,(NP*dt+t0)*1e9-time0);
			  }
			if(WFOnOff[1])
			  {
				sprintf(hisname2,"Ch. 2:x=%5.3e,y=%5.3e,z=%5.3e,U1=%4.2f, U2=%4.2f ",xyz[0][ii],xyz[1][ii],xyz[2][ii],xyz[3][ii],xyz[4][ii]);
					new(entryp2[ii]) TH1F((const Char_t *)(hisname2),(const Char_t *)(hisname2),NP,t0*1e9-time0,(NP*dt+t0)*1e9-time0);
			  }
			if(WFOnOff[2])
			  {
				sprintf(hisname3,"Ch. 3:x=%5.3e,y=%5.3e,z=%5.3e,U1=%4.2f, U2=%4.2f ",xyz[0][ii],xyz[1][ii],xyz[2][ii],xyz[3][ii],xyz[4][ii]);
					new(entryp3[ii]) TH1F((const Char_t *)(hisname3),(const Char_t *)(hisname3),NP,t0*1e9-time0,(NP*dt+t0)*1e9-time0);
			  }


			for(k=0;k<3;k++)
			  {
			   if(WFOnOff[k]==1)
				  {
				for(j=0;j<NP;j++)
				  {
					fscanf(in,"%e",&data); //printf("%e ",data);
					switch(k)
					  {
					case 0: 
					  ((TH1F*)entryp1[ii])->SetBinContent(j+1,data);
					  break;
					case 1: 
					  ((TH1F*)entryp2[ii])->SetBinContent(j+1,data);
					  break;
					case 2: 
					  ((TH1F*)entryp3[ii])->SetBinContent(j+1,data);
					  break;

					  }
				  }
			  //	  printf("k=%d n",k);
				  }
			  }
       }
	 }
  }
    //fclose(in); 
}

 TH1F *PSTCT::GetHA(int ch , int index)  
{
  TH1F *his;
  switch(ch)
    {
    case 0: 
      his=(TH1F *)histo1->At(index); 
      his->SetLineColor(1);
      break;
    case 1: his=(TH1F *)histo2->At(index); 
      his->SetLineColor(2);
      break;
    case 2: his=(TH1F *)histo3->At(index); 
      his->SetLineColor(4);
      break;
    default: his=NULL; break;
    }
return(his);
}

 TH1F *PSTCT::GetHA(Int_t ch , Int_t x, Int_t y, Int_t z, Int_t nu1, Int_t nu2)  
{
  return(GetHA(ch,indx(x,y,z,nu1,nu2)));
}



//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////



 TH1F *PSTCT::ModGetHA(Int_t ch, Int_t index1, Float_t lrange, Float_t rrange)  
{
        //
	// his = his1 - his 2
	// his2 default graf (index2 = 0 default)

  TH1F *his1, *his2, *his;
  
  Float_t B,A,C;
  Int_t i,p,lbin,rbin,q, index2 = RefInd;  
   
  if(RefInd<0 || RefInd>numxyz*NU1*NU2) return GetHA(ch,index1);                             
   
  switch(ch)
    {
    case 0: 
      his1=(TH1F *)histo1->At(index1);
	  his2=(TH1F *)histo1->At(index2);
      his=new TH1F("txt","txt",his1->GetNbinsX(),his1->GetBinCenter(1),his1->GetBinCenter(his1->GetNbinsX()));
	  his->SetLineColor(1);
      break;
    case 1: 
	  his1=(TH1F *)histo2->At(index1);
	  his2=(TH1F *)histo2->At(index2);
      his=new TH1F("txt","txt",his1->GetNbinsX(),his1->GetBinCenter(1),his1->GetBinCenter(his1->GetNbinsX()));		
      his->SetLineColor(2);
      break;
    case 2: 
	  his1=(TH1F *)histo3->At(index1);
	  his2=(TH1F *)histo3->At(index2);
      his=new TH1F("txt","txt",his1->GetNbinsX(),his1->GetBinCenter(1),his1->GetBinCenter(his1->GetNbinsX()));		
      his->SetLineColor(4);
      break;
    default: his=NULL; break;
    }

  if (lrange==-1111 && rrange==-1111)
  {
	  lbin = 1;
	  rbin = his1->GetNbinsX();	  
  }
  else
  {
	for(p=0;p<=his1->GetNbinsX();p++)
		{			
			if (lrange >= his1->GetBinCenter(p) && lrange <= his1->GetBinCenter(p+1)) lbin = p;
			if (rrange >= his1->GetBinCenter(p) && rrange <= his1->GetBinCenter(p+1)) rbin = p;
		}	
  }

		for(i=0;i<his1->GetNbinsX();i++) 
		{	
			if(lbin<=i && i<=rbin)
			{
				B=his1->GetBinContent(i);
				A=his2->GetBinContent(i);		
				C = B-A;
				his->SetBinContent(i,C);			
			}else his->SetBinContent(i,his1->GetBinContent(i));
		}	

	return his;
}

 TH1F *PSTCT::ModGetHA(Int_t ch , Int_t x, Int_t y, Int_t z, Int_t nu1, Int_t nu2, Float_t lrange, Float_t rrange)  
{
  return(ModGetHA(ch,indx(x,y,z,nu1,nu2),lrange,rrange));
}



//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////



 TH1F *PSTCT::AverageGetHA(Int_t ch , Int_t start, Int_t stop, Int_t avrg, Int_t v1, Int_t v2, Int_t v3, Int_t v4)  
{
        // Function is used to get the average wafeform:
        // Int_t ch; channel number - amplifier number
        // Int_t start; start index of the average
        // Int_t stop;  stop index of the average
        // Int_t avrg; decides over which direction it calcultes the average
	//avrg = 0  -> x  
	//avrg = 1  -> y
	//avrg = 2  -> z
	//avrg = 3  -> U1
	//avrg = 4  -> U2
        // Int_t v1,v2,v3,v4; are the indexes of the other paramters in the same order as for average

Int_t i,j,t,n,u1s,u2s,u1e,u2e,q1,q2;
Int_t *s[5];
TH1F *his1,*his2;
Float_t zbir;

switch(avrg)
{
	case 0:			
		s[0] = &i,s[1] = &v1,s[2] = &v2,s[3] = &v3,s[4] = &v4;
		if (stop>Nx) stop = Nx;
		if (start<0) start = 0;
	break;
	case 1:		
		s[0] = &v1,s[1] = &i,s[2] = &v2,s[3] = &v3,s[4] = &v4;
		if (stop>Ny) stop = Ny;
		if (start<0) start = 0;
	break;
	case 2:			
		s[0] = &v1,s[1] = &v2,s[2] = &i,s[3] = &v3,s[4] = &v4;
		if (stop>Nz) stop = Nz;
		if (start<0) start = 0;
	break;
	case 3:	
		s[0] = &v1,s[1] = &v2,s[2] = &v3,s[3] = &i,s[4] = &v4;
		if (stop>NU1) stop = NU1;
		if (start<0) start = 0;
	break;
	case 4:			
		s[0] = &v1,s[1] = &v2,s[2] = &v3,s[3] = &v4,s[4] = &i;
		if (stop>NU2) stop = NU2;
		if (start<0) start = 0;
	break;
}

if ((start==-1111 && stop==-1111) || (start==0 && stop==0)) n=1;
else n=stop-start;
	

	his1=GetHA(ch,*s[0],*s[1],*s[2],*s[3],*s[4]);
	his2 = new TH1F();
	his1->Copy(*his2);
    his2->Reset();

		 for(i=start;i<stop;i++)
		 {			   				
		   his1=GetHA(ch,*s[0],*s[1],*s[2],*s[3],*s[4]);
		   his2->Add(his1);			  
		 }	 
	
	his2->Scale(1./n);

  return his2;
}



//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


 Int_t PSTCT::indx(int x, int y, int z, int nu1, int nu2) 
{
  if(x>Nx-1 || x<0) {printf("index x: out of rangen"); return 0;}
  if(y>Ny-1 || y<0) {printf("index y: out of rangen"); return 0;}
  if(z>Nz-1 || z<0) {printf("index z: out of rangen"); return 0;}
  if(nu1>NU1-1 || nu1<0) {printf("index nu1: out of rangen"); return 0;}
  if(nu2>NU2-1 || nu2<0) {printf("index nu2: out of rangen"); return 0;}

  return( (x+Nx*y+(Nx*Ny)*z)+numxyz*nu2+(NU2*numxyz)*nu1 );
};

 void PSTCT::cords(int x, int y, int z, int nu1, int nu2) 
{
  Int_t ix=indx(x,y,z,nu1,nu2);
  printf("X(%d)=%f, Y(%d)=%f, Z(%d)=%f :: ",ix,xyz[0][ix],ix,xyz[1][ix],ix,xyz[2][ix]);
  printf("%d-> U1=%f, U2=%f, I1=%f, I2=%f :: %4.0fn",ix,xyz[3][ix],xyz[4][ix],xyz[5][ix],xyz[6][ix],xyz[7][ix]);
};


/*TH2F *PSTCT::Draw(Int_t ch, Int_t mode, Int_t nu1, Int_t nu2, Float_t tlow, Float_t thi)
{
	// case: 0 - select x,y
	// case: 1 - select x,z
	// case: 2 - select x,U1
	// case: 3 - select x,U2
	// case: 4 - select y,z
	// case: 5 - select y,U1
	// case: 6 - select y,U2
	// case: 7 - select z,U1
	// case: 8 - select z,U2
	// case: 9 - select U1,U2

	Char_t txt[100];
Float_t integral,max,min;
Int_t i,j,left,right;

 sprintf(txt,"2d his (%f %f)",U1[nu1],U2[nu1]);

TH2F *his2=new TH2F(txt,txt,Nx,x0,Nx*dx+x0,Ny,y0,Ny*dy+y0);
TH1F *his;

    for(j=0;j<Ny;j++)
     for(i=0;i<Nx;i++)
     {
       switch(mode)
	 {
     case 0:		// integral
	   his=GetHA(ch,i,j,0,nu1,nu2);
	   right=his->GetXaxis()->FindBin(thi);
	   left=his->GetXaxis()->FindBin(tlow);
	   integral=his->Integral(left,right);
	   break;
	 case 1:		// max
	    integral=GetHA(ch,i,j,0,nu1,nu2)->GetMaximum();
	    break;
	 case 2:		// min
	    integral=GetHA(ch,i,j,0,nu1,nu2)->GetMinimum();
	    break; 
	 }
       his2->SetBinContent(i,j,integral);
     }
    
	   max=his2->GetMaximum();
	   min=his2->GetMinimum();

	   if(TMath::Abs(max)>TMath::Abs(min)) 
	     his2->SetMinimum(-TMath::Abs(max)); 
	        else 
	     his2->SetMaximum(TMath::Abs(min));

        return his2;
	}*/

 Float_t PSTCT::GetWidth(TH1F *his, Float_t tleft, Float_t tright, Float_t minwidth, Float_t maxwidth, Int_t percentage)
{
	//his -> histogram
	//tleft -> left margin, default 0
	//tright -> right margin, default 25
	//minwidth -> min signal width, default -1111
	//maxwidth -> max signal width, default -1111
	//percentage -> height of min and max, default 50 (50%)

	//return Dt -> ok	
	//return -No of points -> if minwidth & maxwidth are default values
	//return -53 -> error

Float_t q [100];
Int_t j = 0,i,i1,i2,index1,index2,pp;
Float_t val,valn;
Float_t Dt=0,pDt;
Float_t max=0,pmax,max2;
	

	for(i=his->GetXaxis()->FindBin(tleft);i<his->GetXaxis()->FindBin(tright);i++) 
	{
		pmax=his->GetBinContent(i);

		if(max<pmax)
			{
				max=pmax;
				max2=max*percentage/100;			
			}
	}	

	for(i=his->GetXaxis()->FindBin(tleft);i<his->GetXaxis()->FindBin(tright);i++) 
		
	{
		valn=his->GetBinContent(i+1);
		val=his->GetBinContent(i);
		
		if((val<=max2 && valn>max2) || (val>max2 && valn<=max2))	
			 {
				q[j]=his->GetBinCenter(i);				
			    j++;
			 }
	}	

	//printf("Maximum :: %fn",max);

	if(j==2) return Dt = q[1]-q[0];	
	else if(j>2)
	{
		if((minwidth==-1111 || maxwidth==-1111)||(minwidth==0 || maxwidth==0)) return -j;		
		else
		{	
			index1=j-1;
			while(index1>0)			
			{	
				index2 = index1-1;
				while(index2>=0)			
				{							
					pDt=q[index1]-q[index2];
					//printf("q[%d] :: %f, q[%d] :: %f, pDt :: %fn",index1,q[index1],index2,q[index2],pDt);
					if(maxwidth>=pDt && pDt>=minwidth)
					{
						if(pDt>Dt)
						{
							Dt=pDt;
						}
					}
					index2--;
				}
				index1--;
			}			
			return Dt;
		}		
	}else return -53;
}


 TH2F *PSTCT::Draw(Int_t ch, Int_t stype, Int_t mode, Int_t v1, Int_t v2, Int_t v3, Float_t tlow, Float_t thi)
{
        //  Int_t ch; channel number
        //  Int_t stype
	// case: 0 - select x,y
	// case: 1 - select x,z
	// case: 2 - select x,U1
	// case: 3 - select x,U2
	// case: 4 - select y,z
	// case: 5 - select y,U1
	// case: 6 - select y,U2
	// case: 7 - select z,U1
	// case: 8 - select z,U2
	// case: 9 - select U1,U2
        // Int_t mode; 
        //          mode =0 -> Integral
        //          mode =1 -> maximum signal
        //          mode =2 -> minimum signal
        //          mode =3 -> width of the signal
        // Int_t v1,v2,v3; indexes of the remaining three parameters (order always follows x,y,z,U1,U2) 
        // Float_t tlow, thi; Time windos of the plots
        // 
        // 
        // Example :
        // TH2F *plot=aa.Draw(0,4,0,0,0,0,0,20);
        // plot->Draw("COLZ");

Char_t txt[100];
Float_t integral,max,min,width;
Int_t i,j,left,right,t,pX,pY,u1s,u1e,u2s,u2e;
Int_t *s[5];

TH2F *his2;
TH1F *his;

switch(stype)
{
case 0://x,y
	t = indx(0,0,v1,v2,v3);
	sprintf(txt,"2d his (Z[%d]=%f U1[%d]=%f U2[%d]=%f)",v1,xyz[2][t],v2,xyz[3][t],v3,xyz[4][t]);
	his2 = new TH2F(txt,txt,Nx,x0,Nx*dx+x0,Ny,y0,Ny*dy+y0);
	s[0] = &i,s[1] = &j,s[2] = &v1,s[3] = &v2,s[4] = &v3;
	pX=Nx; pY=Ny;
break;
case 1://x,z
	t = indx(0,v1,0,v2,v3);
	sprintf(txt,"2d his (Y[%d]=%f U1[%d]=%f U2[%d]=%f)",v1,xyz[1][t],v2,xyz[3][t],v3,xyz[4][t]);
	his2 = new TH2F(txt,txt,Nx,x0,Nx*dx+x0,Nz,z0,Nz*dz+z0);	
	s[0] = &i,s[1] = &v1,s[2] = &j,s[3] = &v2,s[4] = &v3;
	pX=Nx; pY=Nz;		
break;
case 2://x,U1
	t = indx(0,v1,v2,0,v3);
	sprintf(txt,"2d his (Y[%d]=%f Z[%d]=%f U2[%d]=%f)",v1,xyz[1][t],v2,xyz[2][t],v3,xyz[4][t]);
	if(U1[NU1-1]<U1[0]) {u1s=U1[NU1-1]; u1e=U1[0];} else { u1s=U1[0]; u1e=U1[NU1-1];}
	his2 = new TH2F(txt,txt,(Int_t) Nx,(Float_t) x0,Nx*dx+x0,NU1,u1s,u1e);
	s[0] = &i,s[1] = &v1,s[2] = &v2,s[3] = &j,s[4] = &v3;
	pX=Nx; pY=NU1;		
break;
case 3://x,U2
	t = indx(0,v1,v2,v3,0);
	sprintf(txt,"2d his (Y[%d]=%f Z[%d]=%f U1[%d]=%f)",v1,xyz[1][t],v2,xyz[2][t],v3,xyz[3][t]);
	if(U2[NU2-1]<U2[0]) {u2s=U2[NU2-1]; u2e=U2[0];} else { u2s=U2[0]; u2e=U2[NU2-1];}
	his2 = new TH2F(txt,txt,Nx,x0,Nx*dx+x0,NU2,u2s,u2e);
	s[0] = &i,s[1] = &v1,s[2] = &v2,s[3] = &v3,s[4] = &j;
	pX=Nx; pY=NU2;
break;
case 4://y,z	
	t = indx(v1,0,0,v2,v3);
	sprintf(txt,"2d his (X[%d]=%f U1[%d]=%f U2[%d]=%f)",v1,xyz[2][t],v2,xyz[3][t],v3,xyz[4][t]);
	his2 = new TH2F(txt,txt,Ny,y0,Ny*dy+y0,Nz,z0,Nz*dz+z0);
	s[0] = &v1,s[1] = &i,s[2] = &j,s[3] = &v2,s[4] = &v3;
	pX=Ny; pY=Nz;			
break;
case 5://y,U1	
	t = indx(v1,0,v2,0,v3);
	sprintf(txt,"2d his (X[%d]=%f Z[%d]=%f U2[%d]=%f)",v1,xyz[0][t],v2,xyz[2][t],v3,xyz[4][t]);
	if(U1[NU1-1]<U1[0]) {u1s=U1[NU1-1]; u1e=U1[0];} else { u1s=U1[0]; u1e=U1[NU1-1];}
	his2 = new TH2F(txt,txt,Ny,y0,Ny*dy+y0,NU1,u1s,u1e);
	s[0] = &v1,s[1] = &i,s[2] = &v2,s[3] = &j,s[4] = &v3;
	pX=Ny; pY=NU1;
	
break;
case 6://y,U2	
	t = indx(v1,0,v2,v3,0);
	sprintf(txt,"2d his (X[%d]=%f Z[%d]=%f U1[%d]=%f)",v1,xyz[0][t],v2,xyz[2][t],v3,xyz[3][t]);
	if(U2[NU2-1]<U2[0]) {u2s=U2[NU2-1]; u2e=U2[0];} else { u2s=U2[0]; u2e=U2[NU2-1];}
	his2 = new TH2F(txt,txt,Ny,y0,Ny*dy+y0,NU2,u2s,u2e);
	s[0] = &v1,s[1] = &i,s[2] = &v2,s[3] = &v3,s[4] = &j;
	pX=Ny; pY=NU2;			
break;
case 7://z,U1	
	t = indx(v1,v2,0,0,v3);
	sprintf(txt,"2d his (X[%d]=%f Y[%d]=%f U2[%d]=%f)",v1,xyz[0][t],v2,xyz[1][t],v3,xyz[4][t]);
	if(U1[NU1-1]<U1[0]) {u1s=U1[NU1-1]; u1e=U1[0];} else { u1s=U1[0]; u1e=U1[NU1-1];}
	his2 = new TH2F(txt,txt,Nz,z0,Nz*dz+z0,NU1,u1s,u1e);
	s[0] = &v1,s[1] = &v2,s[2] = &i,s[3] = &j,s[4] = &v3;
	pX=Nz; pY=NU1;			
break;
case 8://z,U2	
	t = indx(v1,v2,0,v3,0);
	sprintf(txt,"2d his (X[%d]=%f Y[%d]=%f U1[%d]=%f)",v1,xyz[0][t],v2,xyz[1][t],v3,xyz[3][t]);
	if(U2[NU2-1]<U2[0]) {u2s=U2[NU2-1]; u2e=U2[0];} else { u2s=U2[0]; u2e=U2[NU2-1];}
	his2 = new TH2F(txt,txt,Nz,z0,Nz*dz+z0,NU2,u2s,u2e);
	s[0] = &v1,s[1] = &v2,s[2] = &i,s[3] = &v3,s[4] = &j;
	pX=Nz; pY=NU2;
break;
case 9://U1,U2	
	t = indx(v1,v2,v3,0,0);
	sprintf(txt,"2d his (X[%d]=%f Y[%d]=%f Z[%d]=%f)",v1,xyz[0][t],v2,xyz[1][t],v3,xyz[2][t]);
	if(U2[NU2-1]<U2[0]) {u2s=U2[NU2-1]; u2e=U2[0];} else { u2s=U2[0]; u2e=U2[NU2-1];}
	if(U1[NU1-1]<U1[0]) {u1s=U1[NU1-1]; u1e=U1[0];} else { u1s=U1[0]; u1e=U1[NU1-1];}
    his2 = new TH2F(txt,txt,NU1,u1s,u1e,NU2,u2s,u2e);
	s[0] = &v1,s[1] = &v2,s[2] = &v3,s[3] = &i,s[4] = &j;
	pX=NU1; pY=NU2;		
break;
}

			for(j=0;j<pY;j++)
			 for(i=0;i<pX;i++)
			 {
				
			   switch(mode)
			 {			
			 case 0:		// integral					
			   his=GetHA(ch,*s[0],*s[1],*s[2],*s[3],*s[4]);			   
			   right=his->GetXaxis()->FindBin(thi);
			   left=his->GetXaxis()->FindBin(tlow);	
			   integral=his->Integral(left,right);
			   //printf("Vrednost integrala :: %fn",integral);
			   break;
			 case 1:		// max	
				his=GetHA(ch,*s[0],*s[1],*s[2],*s[3],*s[4]);
				integral=GetHA(ch,*s[0],*s[1],*s[2],*s[3],*s[4])->GetMaximum();
				width = GetWidth(his,0,40);
				//if(width > 0)
				//{
				//	printf("[i,j,width] = [%d  %d  %f]n",i,j,width);
				//}				
				break;
			 case 2:		// min		  	   
				integral=GetHA(ch,*s[0],*s[1],*s[2],*s[3],*s[4])->GetMinimum();
				break; 
			 case 3:		// Dt
				his=GetHA(ch,*s[0],*s[1],*s[2],*s[3],*s[4]);	
				integral = GetWidth(his,30,60,3,5);
				break; 
			 }
			   //			   printf("i=%d,j=%d int=%f:: %d %d %d %d %dn",i,j,integral,*s[0],*s[1],*s[2],*s[3],*s[4]);
			   //printf("i=%d,j=%d int=%f:: %d %d %d %d %dn",i,j,integral,*s[0],*s[1],*s[2],*s[3],*s[4]);
			   his2->SetBinContent(i+1,j+1,integral);
			 }	   

			   max=his2->GetMaximum();
			   min=his2->GetMinimum();

			   if(TMath::Abs(max)>TMath::Abs(min)) 
				 his2->SetMinimum(-TMath::Abs(max)); 
					else 
				 his2->SetMaximum(TMath::Abs(min));
			   his2->Draw("SURF2");
				return his2;
}


// void PSTCT::DrawMulti(Int_t ch,Int_t dir, Int_t x1,Int_t x2);
// {
//   Int_t i=0,ix;
//   for(i=x1;i<=x2;i++)
//     {
//       switch(dir)
// 	{
// 	case 0: ix=indx(i,
// 	}
      
//     }
// }


 MeasureWF *PSTCT::Projection(Int_t num, Int_t *List)
{

  Float_t Delta;
  Int_t i=0,ix;
  MeasureWF *MWF=new MeasureWF(num);
  TH1F *his;
    for(i=0;i<num;i++)
        {
	ix=indx(List[i*6+1],List[i*6+2],List[i*6+3],List[i*6+4],List[i*6+5]);     
	his=GetHA(List[i*6],ix);
	MWF->AddHisto(i,i,his);
        }
    strcpy(MWF->suffix,"");
    MWF->DrawMode=false;
    strcpy(MWF->prefix,"wf="); 
    return MWF;
}





 MeasureWF *PSTCT::Projection(int ch, int dir,int x,int y,int z, int nu1, int nu2, int num)
{
  // Projection parameters
  // int ch  -> channel number 
  // int dir -> direction of the projection
  // int x   -> x0 of the projection
  // int y   -> y0 of the projection
  // int z   -> z0 of the projection
  // int nu1 -> voltage 1
  // int nu1 -> voltage 2
  // int num -> number of wfs
  
  Float_t Delta;
  Int_t i=0,ix;
  MeasureWF *MWF=new MeasureWF(num);
  TH1F *his;
    for(i=0;i<num;i++)
        {
         switch(dir)  
	   {
	   case 0: ix=indx(i+x,y,z,nu1,nu2); Delta=i*dx; strcpy(MWF->prefix,"x=");   strcpy(MWF->suffix," #mum"); break;
	   case 1: ix=indx(x,y+i,z,nu1,nu2); Delta=i*dy; strcpy(MWF->prefix,"y=");   strcpy(MWF->suffix," #mum"); break;
	   case 2: ix=indx(x,y,z+i,nu1,nu2); Delta=i*dz; strcpy(MWF->prefix,"z=");   strcpy(MWF->suffix," #mum"); break;
	   case 3: ix=indx(x,y,z,nu1+i,nu2); Delta=U1[nu1+i]; strcpy(MWF->prefix,"U1="); strcpy(MWF->suffix," V");    break;
	   case 4: ix=indx(x,y,z,nu1,nu2+i); Delta=U2[nu2+i]; strcpy(MWF->prefix,"U2="); strcpy(MWF->suffix," V");    break;
	   default: ix=indx(i+x,y,z,nu1,nu2); Delta=i*dx; break;
       	   }
	his=GetHA(ch,ix);
	MWF->AddHisto(i,Delta,his);
        }
    MWF->DrawMode=false;
    //    strcpy(MWF->suffix," #mum");
    return MWF;
}



 void PSTCT::DrawList(Int_t num, Int_t *List)
{
  Float_t Delta;
  Int_t i=0,ix;
  TH1F *his;
    for(i=0;i<num;i++)
        {
	  his=GetHA(List[i*6],List[i*6+1],List[i*6+2],List[i*6+3],List[i*6+4],List[i*6+5]);
	  if(i==0) his->DrawCopy(); else his->DrawCopy("SAME");
        }

}

 void PSTCT::DrawList(Int_t num, Int_t *ListC,Int_t *ListP)
{
  // Function draws the list 
  Float_t Delta;
  Int_t i=0,ix;
  TH1F *his;
    for(i=0;i<num;i++)
        {
	  his=GetHA(ListC[i],ListP[i]);
	  if(i==0) his->Draw(); else his->Draw("SAME");
        }

}

 void PSTCT::CorrectBaseLine(Float_t xc)
{
  // Function corrects the baseline (DC offset) of all wafeforms
  // Float_t xc ; time denoting the start of the pulse 
  //              correction factor is calculated from all the bins before xc 
  Int_t right[3],left[3];
  Int_t i,j,k;
  Double_t sum=0,corr[3];
  TH1F *his[3];
  Int_t Num=numxyz*NU1*NU2; //number of all waveforms

  
  for(j=0;j<Num;j++)
    {
    if(j==0)  printf("Baseline correction (%d waveforms) :: ",Num);

      if(WFOnOff[0]==1) his[0]=((TH1F *)histo1->At(j));
      if(WFOnOff[1]==1) his[1]=((TH1F *)histo2->At(j));
      if(WFOnOff[2]==1) his[2]=((TH1F *)histo3->At(j));

    for(i=0;i<3;i++)
      {
	if(WFOnOff[i]==1)
	  {
	    right[i]=his[i]->GetXaxis()->FindBin(xc);
	    left[i]=1;
	    his[i]->Integral(left[i],right[i]);
	    corr[i]=his[i]->Integral(left[i],right[i])/(right[i]-left[i]);
	  }
      }

    if(j%100==0) printf("."); 
      //printf("%d :: Baseline correction = %e , Integral before trigger=%e , Nbins=%d!n",j,corr,his->Integral(left,right),right-left);

    for(k=0;k<3;k++)
      if(WFOnOff[k]==1)
      for(i=1;i<his[k]->GetNbinsX();i++)
	his[k]->SetBinContent(i,his[k]->GetBinContent(i)-corr[k]);
      
    }
  
printf(" finishedn"); 

}


 void PSTCT::PrintInfo()
{
  // Function prints the information about the class and its members
  
  Int_t i,j;
  printf("Format of the file %dn",type); 
  printf("Date and time of the meaurement: %d.%d.%d %d:%d:%dn",Date[0],Date[1],Date[2],Date[3],Date[4],Date[5]);
  printf("Active osciloscope ch: Ch1=%d Ch2=%d Ch3=%dn",WFOnOff[0],WFOnOff[1],WFOnOff[2]);
  printf("Number of points %d (X=%d, Y=%d, Z=%d)n",Nx*Ny*Nz,Nx,Ny,Nz);
  printf("Positions: r0=(%f,%f,%f) dr=(%f,%f,%f) n",x0,y0,z0,dx,dy,dz);
  printf("Time scale: points=%d, t0=%e, dt=%en",NP,t0,dt);
  printf("Voltages: NU1=%d , NU2=%d::n",NU1,NU2);
  for(i=0;i<NU1;i++) 
     for(j=0;j<NU2;j++) 
      printf("U1,U2(%f,%f)::I1,I2(%f,%f)n",U1[i],U2[j],I1[j+NU2*i],I2[j+NU2*i]); printf("n");
    

  
}


 void  PSTCT::swoo(char *a, char *b) {
  // byte swaping (LABVIEW,HPUX g++)<->(LINUX g++, WINNT cl) 
  char c = *a;
  *a = *b; *b = c;
}

 void PSTCT::swooip(float *in, int s) {
 // byte swaping (LABVIEW,HPUX g++)<->(LINUX g++, WINNT cl) 
  char *sr, b;
  while(s--) {
    sr=(char *)in;
    swoo(&sr[0], &sr[3]);
    swoo(&sr[1], &sr[2]);
    in++;
  }
}


ROOT page - Class index - Top of the page

This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.