HLT/comp/AliL3Modeller.cxx

   1 //$Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   4 //*-- Copyright &copy ASV
   5
   6 #include "AliL3StandardIncludes.h"
   7
   8 #include "AliL3Modeller.h"
   9 #include "AliL3MemHandler.h"
  10 #ifdef use_aliroot
  11 #include "AliL3FileHandler.h"
  12 #endif
  13 #include "AliL3TrackArray.h"
  14 #include "AliL3ModelTrack.h"
  15 #include "AliL3DigitData.h"
  16 #include "AliL3Transform.h"
  17
  18 #if GCCVERSION == 3
  19 using namespace std;
  20 #endif
  21
  22 //_____________________________________________________________
  23 // AliL3Modeller
  24 //
  25 // Class for modeling TPC data.
  26 //
  27 // This performs the cluster finding, based on track parameters.
  28 // Basically it propagates the tracks to all padrows, and looks
  29 // for a corresponding cluster. For the moment only cog is calculated,
  30 // and no deconvolution is done.
  31
  32 ClassImp(AliL3Modeller)
  33
  34 AliL3Modeller::AliL3Modeller()
  35 {
  36   fMemHandler=0;
  37   fTracks=0;
  38   fTrackThreshold=0;
  39   fPadOverlap=0;
  40   SetOverlap();
  41   SetTrackThreshold();
  42 }
  43
  44
  45 AliL3Modeller::~AliL3Modeller()
  46 {
  47   if(fMemHandler)
  48     delete fMemHandler;
  49   if(fTracks)
  50     delete fTracks;
  51 }
  52
  53 void AliL3Modeller::Init(Int_t slice,Int_t patch,Char_t *trackdata,Char_t *path,Bool_t houghtracks,Bool_t binary)
  54 {
  55   fSlice = slice;
  56   fPatch = patch;
  57
  58   sprintf(fPath,"%s",path);
  59
  60   fTracks = new AliL3TrackArray("AliL3ModelTrack");
  61
  62   Char_t fname[100];
  63   AliL3MemHandler *file = new AliL3MemHandler();
  64   //sprintf(fname,"%s/tracks_tr_%d_0.raw",trackdata,fSlice); //output tracks from the tracker (no merging)
  65   sprintf(fname,"%s/tracks_ho_%d.raw",trackdata,fSlice);
  66   if(!file->SetBinaryInput(fname))
  67     {
  68       cerr<<"AliL3Modeller::Init : Error opening trackfile: "<<fname<<endl;
  69       return;
  70     }
  71   file->Binary2TrackArray(fTracks);
  72   file->CloseBinaryInput();
  73   delete file;
  74
  75   if(houghtracks)
  76     cout<<"AliL3Modeller is assuming local hough tracksegments!"<<endl;
  77   else
  78     cout<<"AliL3Modeller is assuming global tracks!"<<endl;
  79
  80   if(!houghtracks)
  81     fTracks->QSort();
  82
  83   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
  84     {
  85       AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
  86       if(!track) continue;
  87       track->Init(fSlice,fPatch);
  88
  89       //Only if the tracks has been merged across sector boundaries:
  90       //if(!houghtracks)
  91       //track->Rotate(fSlice,kTRUE); //!!!!!!!!!!!!!!!!!!!
  92
  93       track->CalculateHelix();
  94     }
  95
  96   CalculateCrossingPoints();
  97
  98   CheckForOverlaps();
  99
 100   UInt_t ndigits=0;
 101   AliL3DigitRowData *digits=0;
 102 #ifdef use_aliroot
 103   fMemHandler = new AliL3FileHandler();
 104   fMemHandler->Init(slice,patch);
 105   if(binary == kFALSE)
 106     {
 107       sprintf(fname,"%s/digitfile.root",fPath);
 108       fMemHandler->SetAliInput(fname);
 109       digits = fMemHandler->AliDigits2Memory(ndigits);
 110     }
 111   else
 112     {
 113       sprintf(fname,"%sdigits_%d_%d.raw",fPath,fSlice,fPatch);
 114       if(!fMemHandler->SetBinaryInput(fname))
 115         {
 116           cerr<<"AliL3Modeller::Init : Error opening file "<<fname<<endl;
 117           return;
 118         }
 119       digits=(AliL3DigitRowData*)fMemHandler->CompBinary2Memory(ndigits);
 120     }
 121 #else
 122   fMemHandler = new AliL3MemHandler();
 123   fMemHandler->Init(slice,patch);
 124   if(binary == kFALSE)
 125     {
 126       cerr<<"AliL3Modeller::Init : Compile with AliROOT if you want rootfile as input"<<endl;
 127       return;
 128     }
 129   else
 130     {
 131       sprintf(fname,"%sdigits_%d_%d.raw",fPath,fSlice,fPatch);
 132       if(!fMemHandler->SetBinaryInput(fname))
 133         {
 134           cerr<<"AliL3Modeller::Init : Error opening file "<<fname<<endl;
 135           return;
 136         }
 137     }
 138   digits=(AliL3DigitRowData*)fMemHandler->CompBinary2Memory(ndigits);
 139 #endif
 140
 141   SetInputData(digits);
 142 }
 143
 144 void AliL3Modeller::FindClusters()
 145 {
 146
 147   if(!fTracks)
 148     {
 149       cerr<<"AliL3Modeller::Process : No tracks"<<endl;
 150       return;
 151     }
 152   if(!fRowData)
 153     {
 154       cerr<<"AliL3Modeller::Process : No data "<<endl;
 155       return;
 156     }
 157
 158   AliL3DigitRowData *rowPt = fRowData;
 159   AliL3DigitData *digPt=0;
 160
 161   Int_t ntimes = AliL3Transform::GetNTimeBins()+1;
 162   Int_t npads = AliL3Transform::GetNPads(AliL3Transform::GetLastRow(fPatch))+1;//Max num of pads.
 163   Int_t bounds = ntimes*npads;
 164   Digit *row = new Digit[bounds];
 165
 166   Int_t seq_charge;
 167   Int_t pad,time,index;
 168   Short_t charge;
 169   Cluster cluster;
 170
 171   for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
 172     {
 173       fCurrentPadRow = i;
 174       memset((void*)row,0,ntimes*npads*sizeof(Digit));
 175       digPt = (AliL3DigitData*)rowPt->fDigitData;
 176       //cout<<"Loading row "<<i<<" with "<<(Int_t)rowPt->fNDigit<<" digits"<<endl;
 177       for(UInt_t j=0; j<rowPt->fNDigit; j++)
 178         {
 179           pad = digPt[j].fPad;
 180           time = digPt[j].fTime;
 181           charge = digPt[j].fCharge;
 182           row[ntimes*pad+time].fCharge = charge;
 183           row[ntimes*pad+time].fUsed = kFALSE;
 184           //cout<<"Row "<<i<<" pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
 185         }
 186
 187       for(Int_t k=0; k<fTracks->GetNTracks(); k++)
 188         {
 189           AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(k);
 190           if(!track) continue;
 191
 192           if(track->GetPadHit(i)<0 || track->GetTimeHit(i)<0 || track->GetOverlap(i)>=0)
 193             {
 194               track->SetCluster(i,0,0,0,0,0,0); //The track has left the patch.
 195               continue;
 196             }
 197
 198           Int_t hitpad = (Int_t)rint(track->GetPadHit(i));
 199           Int_t hittime = (Int_t)rint(track->GetTimeHit(i));
 200           //cout<<"Checking track on row "<<i<<" with pad "<<hitpad<<" time "<<hittime<<endl;
 201           pad = hitpad;
 202           time = hittime;
 203           Int_t padsign=-1;
 204           Int_t timesign=-1;
 205
 206           memset(&cluster,0,sizeof(Cluster));
 207
 208           Int_t npads=0;
 209           while(1)//Process this padrow
 210             {
 211               if(pad < 0 || pad >= AliL3Transform::GetNPads(i))
 212                 {
 213                   //cout<<"Pad = "<<pad<<" on row "<<i<<endl;
 214                   FillCluster(track,&cluster,i,npads);
 215                   break;
 216                 }
 217               seq_charge=0;
 218               timesign=-1;
 219               time = hittime;
 220
 221               while(1) //Process sequence on this pad:
 222                 {
 223                   if(time < 0) break;
 224                   index = ntimes*pad + time;
 225                   if(index < 0 || index >= bounds)
 226                     {
 227                       cerr<<"AliL3Modeller::FindClusters : Index out of range : "<<index
 228                         <<" on row "<<i<<" pad "<<pad<<" time "<<time<<endl;
 229                       break;
 230                     }
 231
 232                   charge = row[index].fCharge;
 233                   if(charge==0 && timesign==-1) //zero charge on this timebin, perform checks:
 234                     {
 235                       if(seq_charge==0 && abs(time-hittime) <= fTimeOverlap) //No charge found on this pad, look further.
 236                         {
 237                           time--;
 238                           continue;
 239                         }
 240                       else //Boundary reached, or we have found one end of the sequence,->start looking in the other time direction
 241                         {
 242                           time = hittime+1;
 243                           timesign=1;
 244                           continue;
 245                         }
 246                     }
 247                   else if(charge==0 && timesign==1)//zero charge on this timebin, perform checks:
 248                     {
 249                       if(seq_charge==0 && abs(time-hittime) <= fTimeOverlap)//No charge found on this pad, look further
 250                         {
 251                           time++;
 252                           continue;
 253                         }
 254                       else //Boundary reached, or we have found the other end of the sequence, stop looking on this pad.
 255                         {
 256                           //if(fCurrentPadRow==31)
 257                           //   cerr<<"Breaking off at pad "<<pad<<" and time "<<time<<endl;
 258                           break;
 259                         }
 260                     }
 261
 262                   if(row[ntimes*pad+time].fUsed==kTRUE) //Don't use digits several times. This leads to mult. rec.tracks.
 263                     {
 264                       time += timesign;
 265                       continue;
 266                     }
 267
 268                   seq_charge += charge;
 269
 270                   //Update the cluster parameters with this timebin
 271                   cluster.fTime += time*charge;
 272                   cluster.fPad += pad*charge;
 273                   cluster.fCharge += charge;
 274                   cluster.fSigmaY2 += pad*pad*charge;
 275                   cluster.fSigmaZ2 += time*time*charge;
 276
 277                   row[ntimes*pad+time].fUsed = kTRUE;
 278                   time += timesign;
 279                 }
 280
 281
 282               if(seq_charge)//There was something on this pad, so keep looking on the neighbouring pad
 283                 {
 284                   pad += padsign;
 285                   npads++;
 286                 }
 287               else //Nothing more on this pad, goto next pad
 288                 {
 289                   if(padsign==-1)
 290                     {
 291                       if(cluster.fCharge==0 && abs(pad-hitpad) <= fPadOverlap && pad > 0)
 292                         {
 293                           pad--; //In this case, we haven't found anything yet,
 294                         }        //so we will try to expand our search within the natural boundaries.
 295                       else
 296                         {
 297                           pad=hitpad+1;
 298                           padsign=1;
 299                         }
 300                       continue;
 301                     }
 302
 303                   else if(padsign==1)
 304                     {
 305                       if(cluster.fCharge==0 && abs(pad-hitpad) <= fPadOverlap && pad < AliL3Transform::GetNPads(i)-2)
 306                         {
 307                           pad++;     //In this case, we haven't found anything yet,
 308                           continue;  //so we will try to expand our search within the natural boundaries.
 309                         }
 310                       else //We are out of range, or cluster if finished.
 311                         {
 312                           //if(fCurrentPadRow==31)
 313                           //cout<<"Out of range; charge "<<cluster.fCharge<<" paddiff "<<abs(pad-hitpad)<<endl;
 314                           FillCluster(track,&cluster,i,npads);
 315                           break;
 316                         }
 317                     }
 318                   else //Nothing more in this cluster
 319                     {
 320                       //if(fCurrentPadRow==31)
 321                       //cout<<"Filling final cluster"<<endl;
 322                       FillCluster(track,&cluster,i,npads);
 323                       break;
 324                     }
 325                 }
 326             }
 327           //cout<<"done"<<endl;
 328         }
 329       FillZeros(rowPt,row);
 330       fMemHandler->UpdateRowPointer(rowPt);
 331     }
 332   delete [] row;
 333   cout<<"done processing"<<endl;
 334
 335
 336   //Debug:
 337   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
 338     {
 339       AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
 340       if(!track) continue;
 341       if(track->GetNClusters() != AliL3Transform::GetNRows(fPatch))
 342         cerr<<endl<<"Mismatching hitcounts; nclusters: "<<track->GetNClusters()<<" nrows "<<AliL3Transform::GetNRows(fPatch)<<endl<<endl;
 343     }
 344
 345 }
 346
 347 void AliL3Modeller::FillCluster(AliL3ModelTrack *track,Cluster *cluster,Int_t row,Int_t npads)
 348 {
 349   if(cluster->fCharge==0)
 350     {
 351       track->SetCluster(row,0,0,0,0,0,0);
 352       return;
 353     }
 354   Float_t fcharge = (Float_t)cluster->fCharge;
 355   Float_t fpad = ((Float_t)cluster->fPad/fcharge);
 356   Float_t ftime = ((Float_t)cluster->fTime/fcharge);
 357   Float_t sigmaY2,sigmaZ2;
 358   CalcClusterWidth(cluster,sigmaY2,sigmaZ2);
 359   track->SetCluster(row,fpad,ftime,fcharge,sigmaY2,sigmaZ2,npads);
 360 }
 361
 362 void AliL3Modeller::FillZeros(AliL3DigitRowData *rowPt,Digit *row)
 363 {
 364   //Fill zero where data has been used.
 365
 366   Int_t ntimes = AliL3Transform::GetNTimeBins()+1;
 367   AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
 368   for(UInt_t j=0; j<rowPt->fNDigit; j++)
 369     {
 370       Int_t pad = digPt[j].fPad;
 371       Int_t time = digPt[j].fTime;
 372       if(row[ntimes*pad+time].fUsed==kTRUE)
 373         digPt[j].fCharge = 0;
 374     }
 375 }
 376
 377 void AliL3Modeller::WriteRemaining()
 378 {
 379   //Write remaining (nonzero) digits to file.
 380
 381   AliL3DigitRowData *rowPt;
 382   rowPt = (AliL3DigitRowData*)fRowData;
 383   Int_t digitcount=0;
 384   Int_t ndigits[(AliL3Transform::GetNRows(fPatch))];
 385   for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
 386     {
 387       AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
 388       ndigits[(i-AliL3Transform::GetFirstRow(fPatch))]=0;
 389       for(UInt_t j=0; j<rowPt->fNDigit; j++)
 390         {
 391           if(digPt[j].fCharge==0) continue;
 392           digitcount++;
 393           ndigits[(i-AliL3Transform::GetFirstRow(fPatch))]++;
 394         }
 395       //cout<<"Difference "<<(int)ndigits[(i-AliL3Transform::GetFirstRow(fPatch))]<<" "<<(int)rowPt->fNDigit<<endl;
 396       fMemHandler->UpdateRowPointer(rowPt);
 397     }
 398
 399   Int_t size = digitcount*sizeof(AliL3DigitData) + AliL3Transform::GetNRows(fPatch)*sizeof(AliL3DigitRowData);
 400   Byte_t *data = new Byte_t[size];
 401   memset(data,0,size);
 402   AliL3DigitRowData *tempPt = (AliL3DigitRowData*)data;
 403   rowPt = (AliL3DigitRowData*)fRowData;
 404
 405   for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
 406     {
 407       Int_t localcount=0;
 408       tempPt->fRow = i;
 409       tempPt->fNDigit = ndigits[(i-AliL3Transform::GetFirstRow(fPatch))];
 410       AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
 411       for(UInt_t j=0; j<rowPt->fNDigit; j++)
 412         {
 413           if(digPt[j].fCharge==0) continue;
 414           if(localcount >= ndigits[(i-AliL3Transform::GetFirstRow(fPatch))])
 415             {
 416               cerr<<"AliL3Modeller::WriteRemaining : Digitarray out of range!!"<<endl;
 417               return;
 418             }
 419           tempPt->fDigitData[localcount].fCharge = digPt[j].fCharge;
 420           tempPt->fDigitData[localcount].fPad = digPt[j].fPad;
 421           tempPt->fDigitData[localcount].fTime = digPt[j].fTime;
 422
 423           localcount++;
 424         }
 425       if(ndigits[(i-AliL3Transform::GetFirstRow(fPatch))] != localcount)
 426         {
 427           cerr<<"AliL3Modeller::WriteRemaining : Mismatch in digitcount"<<endl;
 428           return;
 429         }
 430       fMemHandler->UpdateRowPointer(rowPt);
 431       Byte_t *tmp = (Byte_t*)tempPt;
 432       Int_t size = sizeof(AliL3DigitRowData) + ndigits[(i-AliL3Transform::GetFirstRow(fPatch))]*sizeof(AliL3DigitData);
 433       tmp += size;
 434       tempPt = (AliL3DigitRowData*)tmp;
 435     }
 436
 437   Char_t fname[100];
 438   AliL3MemHandler *mem = new AliL3MemHandler();
 439   sprintf(fname,"%s/comp/remains_%d_%d.raw",fPath,fSlice,fPatch);
 440   mem->SetBinaryOutput(fname);
 441   mem->Memory2CompBinary((UInt_t)AliL3Transform::GetNRows(fPatch),(AliL3DigitRowData*)data);
 442   mem->CloseBinaryOutput();
 443   delete mem;
 444   delete [] data;
 445 }
 446
 447
 448 void AliL3Modeller::CalculateCrossingPoints()
 449 {
 450   cout<<"Calculating crossing points on "<<fTracks->GetNTracks()<<" tracks"<<endl;
 451   if(!fTracks)
 452     {
 453       cerr<<"AliL3Modeller::CalculateCrossingPoints(): No tracks"<<endl;
 454       return;
 455     }
 456   Float_t hit[3];
 457
 458   Int_t sector,row;
 459   for(Int_t i=AliL3Transform::GetLastRow(fPatch); i>=AliL3Transform::GetFirstRow(fPatch); i--)
 460     {
 461       for(Int_t j=0; j<fTracks->GetNTracks(); j++)
 462         {
 463           AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(j);
 464           if(!track) continue;
 465
 466           if(!track->GetCrossingPoint(i,hit))
 467             {
 468               cerr<<"AliL3Modeller::CalculateCrossingPoints : Track "<<j<<" does not intersect row "<<i<<" :"<<endl<<
 469                 "First point "<<track->GetFirstPointX()<<
 470                 " nhits "<<track->GetNHits()<<endl;//" tgl "<<track->GetTgl()<<" psi "<<track->GetPsi()<<" charge "<<track->GetCharge()<<endl;
 471                 //"Center "<<track->GetCenterX()<<" "<<track->GetCenterY()<<endl<<endl<<
 472                 //"--------"<<endl;
 473               fTracks->Remove(j);
 474               continue;
 475             }
 476           //cout<<"X "<<hit[0]<<" Y "<<hit[1]<<" Z "<<hit[2]<<" tgl "<<track->GetTgl()<<endl;
 477
 478           AliL3Transform::Slice2Sector(fSlice,i,sector,row);
 479           AliL3Transform::Local2Raw(hit,sector,row);
 480           //cout<<"Pad "<<hit[1]<<" time "<<hit[2]<<" in sector "<<sector<<" row "<<row<<endl;
 481           if(hit[1]<0 || hit[1]>AliL3Transform::GetNPads(i) ||
 482              hit[2]<0 || hit[2]>AliL3Transform::GetNTimeBins())
 483             {//Track is leaving the patch, so flag the track hits (<0)
 484               track->SetPadHit(i,-1);
 485               track->SetTimeHit(i,-1);
 486               continue;
 487             }
 488
 489
 490           track->SetPadHit(i,hit[1]);
 491           track->SetTimeHit(i,hit[2]);
 492
 493           //if(hit[1]<0 || hit[2]>445)
 494           //if(hit[2]<0 || hit[2]>445)
 495           //cout<<"pad "<<hit[1]<<" time "<<hit[2]<<" pt "<<track->GetPt()<<" psi "<<track->GetPsi()<<" tgl "<<track->GetTgl()<<" firstpoint "<<track->GetFirstPointX()<<" "<<track->GetFirstPointY()<<" "<<track->GetFirstPointZ()<<endl;
 496           //cout<<"Crossing pad "<<hit[1]<<" time "<<hit[2]<<endl;
 497         }
 498     }
 499   fTracks->Compress();
 500   cout<<"And there are "<<fTracks->GetNTracks()<<" tracks remaining"<<endl;
 501 }
 502
 503 void AliL3Modeller::CheckForOverlaps()
 504 {
 505   //Flag the tracks that overlap
 506
 507   cout<<"Checking for overlaps...";
 508   Int_t counter=0;
 509   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
 510     {
 511       AliL3ModelTrack *track1 = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
 512       if(!track1) continue;
 513       for(Int_t j=i+1; j<fTracks->GetNTracks(); j++)
 514         {
 515           AliL3ModelTrack *track2 = (AliL3ModelTrack*)fTracks->GetCheckedTrack(j);
 516           if(!track2) continue;
 517           for(Int_t k=AliL3Transform::GetFirstRow(fPatch); k<=AliL3Transform::GetLastRow(fPatch); k++)
 518             {
 519               if(track1->GetPadHit(k)<0 || track1->GetTimeHit(k)<0 ||
 520                  track2->GetPadHit(k)<0 || track2->GetTimeHit(k)<0)
 521                 continue;
 522
 523               if(track1->GetOverlap(k)>=0 || track2->GetOverlap(k)>=0) continue;
 524
 525               if(abs((Int_t)rint(track1->GetPadHit(k))-(Int_t)rint(track2->GetPadHit(k))) <= fPadOverlap &&
 526                  abs((Int_t)rint(track1->GetTimeHit(k))-(Int_t)rint(track2->GetTimeHit(k))) <= fTimeOverlap)
 527                 {
 528                   track2->SetOverlap(k,i);
 529                   //track1->SetOverlap(k,j);
 530                   counter++;
 531                 }
 532             }
 533         }
 534     }
 535   cout<<"found "<<counter<<" done"<<endl;
 536 }
 537
 538
 539 void AliL3Modeller::CalcClusterWidth(Cluster *cl,Float_t &sigmaY2,Float_t &sigmaZ2)
 540 {
 541
 542   Float_t padw,timew;
 543
 544   padw = AliL3Transform::GetPadPitchWidth(fPatch);
 545
 546   Float_t charge = (Float_t)cl->fCharge;
 547   Float_t pad = (Float_t)cl->fPad/charge;
 548   Float_t time = (Float_t)cl->fTime/charge;
 549   Float_t s2 = (Float_t)cl->fSigmaY2/charge - pad*pad;
 550
 551   //Save the sigmas in pad and time:
 552
 553   sigmaY2 = (s2);// + 1./12);//*padw*padw;
 554
 555   /*Constants added by offline
 556     if(s2 != 0)
 557     {
 558     sigmaY2 = sigmaY2*0.108;
 559     if(fPatch<3)
 560     sigmaY2 = sigmaY2*2.07;
 561     }
 562   */
 563
 564   s2 = (Float_t)cl->fSigmaZ2/charge - time*time;
 565   timew = AliL3Transform::GetZWidth();
 566   sigmaZ2 = (s2);// +1./12);//*timew*timew;
 567
 568
 569
 570   /*Constants added by offline
 571     if(s2 != 0)
 572     {
 573     sigmaZ2 = sigmaZ2*0.169;
 574     if(fPatch < 3)
 575     sigmaZ2 = sigmaZ2*1.77;
 576     }
 577   */
 578 }