HLT/comp/AliL3Modeller.cxx

   1 //$Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo$fi.uib.no>
   4 //*-- Copyright &copy ASV
   5
   6 #include <stream.h>
   7 #include <iostream.h>
   8 #include <math.h>
   9
  10 #include "AliL3Modeller.h"
  11 #include "AliL3MemHandler.h"
  12 #include "AliL3TrackArray.h"
  13 #include "AliL3ModelTrack.h"
  14 #include "AliL3DigitData.h"
  15 #include "AliL3Transform.h"
  16
  17 #include "AliL3Defs.h"
  18
  19 ClassImp(AliL3Modeller)
  20
  21 AliL3Modeller::AliL3Modeller()
  22 {
  23   fMemHandler=0;
  24   fTracks=0;
  25   fTransform=0;
  26 }
  27
  28
  29 AliL3Modeller::~AliL3Modeller()
  30 {
  31   if(fMemHandler)
  32     delete fMemHandler;
  33   if(fTracks)
  34     delete fTracks;
  35   if(fTransform)
  36     delete fTransform;
  37
  38 }
  39
  40 void AliL3Modeller::Init(Int_t slice,Int_t patch,Char_t *path)
  41 {
  42   fSlice = slice;
  43   fPatch = patch;
  44   fPadOverlap=4;
  45   fTimeOverlap=4;
  46   fTransform = new AliL3Transform();
  47   fTracks = new AliL3TrackArray("AliL3ModelTrack");
  48
  49   AliL3MemHandler *file = new AliL3MemHandler();
  50   if(!file->SetBinaryInput("tracks.raw"))
  51     {
  52       cerr<<"AliL3Modeller::Init : Error opening trackfile"<<endl;
  53       return;
  54     }
  55   file->Binary2TrackArray(fTracks);
  56   file->CloseBinaryInput();
  57   delete file;
  58
  59   fTracks->QSort();
  60   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
  61     {
  62       AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
  63       if(!track) continue;
  64       track->Init(fSlice,fPatch);
  65       track->Rotate(fSlice,kTRUE);
  66       track->CalculateHelix();
  67     }
  68
  69   CalculateCrossingPoints();
  70   CheckForOverlaps();
  71
  72   fMemHandler = new AliL3MemHandler();
  73   Char_t fname[100];
  74   sprintf(fname,"%sdigits_%d_%d.raw",path,fSlice,fPatch);
  75   if(!fMemHandler->SetBinaryInput(fname))
  76     {
  77       cerr<<"AliL3Modeller::Init : Error opening file "<<fname<<endl;
  78       return;
  79     }
  80   UInt_t ndigits;
  81   AliL3DigitRowData *digits=(AliL3DigitRowData*)fMemHandler->CompBinary2Memory(ndigits);
  82
  83   SetInputData(digits);
  84 }
  85
  86 void AliL3Modeller::Process()
  87 {
  88
  89   if(!fTracks)
  90     {
  91       cerr<<"AliL3Modeller::Process : No tracks"<<endl;
  92       return;
  93     }
  94   if(!fRowData)
  95     {
  96       cerr<<"AliL3Modeller::Process : No data "<<endl;
  97       return;
  98     }
  99
 100   AliL3DigitRowData *rowPt = fRowData;
 101   AliL3DigitData *digPt=0;
 102
 103   Int_t ntimes = fTransform->GetNTimeBins()+1;
 104   Int_t npads = fTransform->GetNPads(NRows[fPatch][1])+1;//Max num of pads.
 105   Digit *row = new Digit[(ntimes)*(npads)];
 106
 107   Int_t seq_charge;
 108   Int_t pad,time;
 109   Short_t charge;
 110   Cluster cluster;
 111
 112   for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
 113     {
 114       memset((void*)row,0,ntimes*npads*sizeof(Digit));
 115       digPt = (AliL3DigitData*)rowPt->fDigitData;
 116       for(UInt_t j=0; j<rowPt->fNDigit; j++)
 117         {
 118           pad = digPt[j].fPad;
 119           time = digPt[j].fTime;
 120           charge = digPt[j].fCharge;
 121           row[ntimes*pad+time].fCharge = charge;
 122           row[ntimes*pad+time].fUsed = kFALSE;
 123         }
 124
 125       for(Int_t k=0; k<fTracks->GetNTracks(); k++)
 126         {
 127           AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(k);
 128           if(!track) continue;
 129           if(track->GetOverlap()>=0) continue;//Track is overlapping
 130           if(track->GetPadHit(i)<0 || track->GetTimeHit(i)<0)
 131             {
 132               track->SetCluster(0,0,0,0,0); //The track has left the patch.
 133               continue;
 134             }
 135
 136           Int_t hitpad = (Int_t)rint(track->GetPadHit(i));
 137           Int_t hittime = (Int_t)rint(track->GetTimeHit(i));
 138           //cout<<"Checking track with pad "<<hitpad<<" time "<<hittime<<endl;
 139           pad = hitpad;
 140           time = hittime;
 141           Int_t padsign=-1;
 142           Int_t timesign=-1;
 143
 144           memset(&cluster,0,sizeof(Cluster));
 145
 146           while(1)//Process this padrow
 147             {
 148               seq_charge=0;
 149               timesign=-1;
 150               time = hittime;
 151               while(1) //Process sequence on this pad:
 152                 {
 153                   charge = row[ntimes*pad+time].fCharge;
 154                   if(charge==0 && timesign==-1)
 155                     {time=hittime+1; timesign=1; continue;}
 156                   else if(charge==0 && timesign==1)
 157                     break;
 158                   //cout<<"Doing pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
 159
 160                   seq_charge += charge;
 161
 162                   cluster.fTime += time*charge;
 163                   cluster.fPad += pad*charge;
 164                   cluster.fCharge += charge;
 165                   cluster.fSigmaY2 += pad*pad*charge;
 166                   cluster.fSigmaZ2 += time*time*charge;
 167
 168                   row[ntimes*pad+time].fUsed = kTRUE;
 169                   time += timesign;
 170                 }
 171               //cout<<"Finished on pad "<<pad<<" and time "<<time<<endl;
 172               if(seq_charge)
 173                 pad += padsign;
 174               else //Nothing more on this pad, goto next pad
 175                 {
 176                   if(padsign==-1)
 177                     {
 178                       if(cluster.fCharge==0 && abs(pad-hitpad) < fPadOverlap/2)
 179                         {
 180                           pad--; //In this case, we haven't found anything yet,
 181                         }        //so we will try to expand our search within the natural boundaries.
 182                       else
 183                         {
 184                           pad=hitpad+1;
 185                           padsign=1;
 186                         }
 187                       continue;
 188                     }
 189                   else if(padsign==1 && cluster.fCharge==0 && abs(pad-hitpad) < fPadOverlap/2)
 190                     {
 191                       pad++;
 192                       continue;
 193                     }
 194                   else //Nothing more in this cluster
 195                     {
 196                       Float_t fcharge = (Float_t)cluster.fCharge;
 197                       Float_t fpad = ((Float_t)cluster.fPad/fcharge);
 198                       Float_t ftime = ((Float_t)cluster.fTime/fcharge);
 199                       Float_t sigmaY2,sigmaZ2;
 200                       CalcClusterWidth(&cluster,sigmaY2,sigmaZ2);
 201                       //cout<<"row "<<i<<" charge "<<fcharge<<endl;
 202                       track->SetCluster(fpad,ftime,fcharge,sigmaY2,sigmaZ2);
 203                       break;
 204                     }
 205                 }
 206               // pad += padsign;
 207             }
 208         }
 209       FillZeros(rowPt,row);
 210       fMemHandler->UpdateRowPointer(rowPt);
 211     }
 212   delete [] row;
 213
 214 }
 215
 216 void AliL3Modeller::FillZeros(AliL3DigitRowData *rowPt,Digit *row)
 217 {
 218   //Fill zero where data has been used.
 219
 220   Int_t ntimes = fTransform->GetNTimeBins()+1;
 221   AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
 222   for(UInt_t j=0; j<rowPt->fNDigit; j++)
 223     {
 224       Int_t pad = digPt[j].fPad;
 225       Int_t time = digPt[j].fTime;
 226       if(row[ntimes*pad+time].fUsed==kTRUE)
 227         digPt[j].fCharge = 0;
 228     }
 229 }
 230
 231 void AliL3Modeller::WriteRemaining(Char_t *output)
 232 {
 233   //Write remaining (nonzero) digits to file.
 234
 235   cout<<"Writing remaining data to file "<<output<<endl;
 236   AliL3DigitRowData *rowPt;
 237   rowPt = (AliL3DigitRowData*)fRowData;
 238   Int_t digitcount=0;
 239   Int_t ndigits[(NumRows[fPatch])];
 240   for(Int_t i=NRows[fPatch][0]; i<NRows[fPatch][1]; i++)
 241     {
 242       AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
 243       ndigits[(i-NRows[fPatch][0])]=0;
 244       for(UInt_t j=0; j<rowPt->fNDigit; j++)
 245         {
 246           if(digPt[j].fCharge==0) continue;
 247           digitcount++;
 248           ndigits[(i-NRows[fPatch][0])]++;
 249         }
 250       //cout<<"Difference "<<(int)ndigits[(i-NRows[fPatch][0])]<<" "<<(int)rowPt->fNDigit<<endl;
 251       fMemHandler->UpdateRowPointer(rowPt);
 252     }
 253
 254   Int_t size = digitcount*sizeof(AliL3DigitData) + NumRows[fPatch]*sizeof(AliL3DigitRowData);
 255   Byte_t *data = new Byte_t[size];
 256   memset(data,0,size);
 257   AliL3DigitRowData *tempPt = (AliL3DigitRowData*)data;
 258   rowPt = (AliL3DigitRowData*)fRowData;
 259
 260   for(Int_t i=NRows[fPatch][0]; i<NRows[fPatch][1]; i++)
 261     {
 262       Int_t localcount=0;
 263       tempPt->fRow = i;
 264       tempPt->fNDigit = ndigits[(i-NRows[fPatch][0])];
 265       AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
 266       for(UInt_t j=0; j<rowPt->fNDigit; j++)
 267         {
 268           if(digPt[j].fCharge==0) continue;
 269           if(localcount >= ndigits[(i-NRows[fPatch][0])])
 270             {
 271               cerr<<"AliL3Modeller::WriteRemaining : Digitarray out of range!!"<<endl;
 272               return;
 273             }
 274           tempPt->fDigitData[localcount].fCharge = digPt[j].fCharge;
 275           tempPt->fDigitData[localcount].fPad = digPt[j].fPad;
 276           tempPt->fDigitData[localcount].fTime = digPt[j].fTime;
 277           localcount++;
 278         }
 279       if(ndigits[(i-NRows[fPatch][0])] != localcount)
 280         {
 281           cerr<<"AliL3Modeller::WriteRemaining : Mismatch in digitcount"<<endl;
 282           return;
 283         }
 284       fMemHandler->UpdateRowPointer(rowPt);
 285       Byte_t *tmp = (Byte_t*)tempPt;
 286       Int_t size = sizeof(AliL3DigitRowData) + ndigits[(i-NRows[fPatch][0])]*sizeof(AliL3DigitData);
 287       tmp += size;
 288       tempPt = (AliL3DigitRowData*)tmp;
 289     }
 290
 291   AliL3MemHandler *mem = new AliL3MemHandler();
 292   mem->SetBinaryOutput(output);
 293   mem->Memory2CompBinary((UInt_t)NumRows[fPatch],(AliL3DigitRowData*)data);
 294   mem->CloseBinaryOutput();
 295   delete mem;
 296 }
 297
 298
 299 void AliL3Modeller::CalculateCrossingPoints()
 300 {
 301   cout<<"Calculating crossing points on "<<fTracks->GetNTracks()<<" tracks"<<endl;
 302   if(!fTracks)
 303     {
 304       cerr<<"AliL3Modeller::CalculateCrossingPoints(): No tracks"<<endl;
 305       return;
 306     }
 307   Float_t hit[3];
 308   for(Int_t i=NRows[fPatch][1]; i>=NRows[fPatch][0]; i--)
 309     {
 310       for(Int_t j=0; j<fTracks->GetNTracks(); j++)
 311         {
 312           AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(j);
 313           if(!track) continue;
 314           if(track->GetNHits() < 100)
 315             fTracks->Remove(j);
 316           if(!track->GetCrossingPoint(i,hit))
 317             {
 318               cerr<<"AliL3Modeller::CalculateCrossingPoints : Track does not intersect line "<<endl;
 319               continue;
 320             }
 321           //cout<<" x "<<track->GetPointX()<<" y "<<track->GetPointY()<<" z "<<track->GetPointZ()<<endl;
 322
 323           fTransform->Local2Raw(hit,fSlice,i);
 324           if(hit[1]<0 || hit[1]>fTransform->GetNPads(i) ||
 325              hit[2]<0 || hit[2]>fTransform->GetNTimeBins())
 326             {//Track is leaving the patch, so flag the track hits (<0)
 327               track->SetPadHit(i,-1);
 328               track->SetTimeHit(i,-1);
 329               continue;
 330             }
 331
 332           track->SetPadHit(i,hit[1]);
 333           track->SetTimeHit(i,hit[2]);
 334
 335           //if(hit[1]<0 || hit[2]>445)
 336           //if(hit[2]<0 || hit[2]>445)
 337           //cout<<"pad "<<hit[1]<<" time "<<hit[2]<<" pt "<<track->GetPt()<<" psi "<<track->GetPsi()<<" tgl "<<track->GetTgl()<<" firstpoint "<<track->GetFirstPointX()<<" "<<track->GetFirstPointY()<<" "<<track->GetFirstPointZ()<<endl;
 338           //cout<<"Crossing pad "<<hit[1]<<" time "<<hit[2]<<endl;
 339         }
 340     }
 341   fTracks->Compress();
 342   //cout<<"And there are "<<fTracks->GetNTracks()<<" tracks remaining"<<endl;
 343 }
 344
 345 void AliL3Modeller::CheckForOverlaps()
 346 {
 347   //Flag the tracks that overlap
 348
 349   cout<<"Checking for overlaps"<<endl;
 350
 351   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
 352     {
 353       AliL3ModelTrack *track1 = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
 354       if(!track1) continue;
 355       for(Int_t j=i+1; j<fTracks->GetNTracks(); j++)
 356         {
 357           AliL3ModelTrack *track2 = (AliL3ModelTrack*)fTracks->GetCheckedTrack(j);
 358           if(!track2) continue;
 359           for(Int_t k=NRows[fPatch][0]; k<NRows[fPatch][1]; k++)
 360             {
 361               if(track1->GetPadHit(k)<0 || track1->GetTimeHit(k)<0 ||
 362                  track2->GetPadHit(k)<0 || track2->GetTimeHit(k)<0)
 363                 continue;
 364               if(fabs(track1->GetPadHit(k)-track2->GetPadHit(k)) < fPadOverlap &&
 365                  fabs(track1->GetTimeHit(k)-track2->GetTimeHit(k)) < fTimeOverlap)
 366                 {
 367                   track1->SetOverlap(j);
 368                   track2->SetOverlap(i);
 369                 }
 370             }
 371         }
 372     }
 373
 374 }
 375
 376
 377 void AliL3Modeller::CalcClusterWidth(Cluster *cl,Float_t &sigmaY2,Float_t &sigmaZ2)
 378 {
 379
 380   Float_t padw,timew;
 381   if(fPatch < 3)
 382     padw = fTransform->GetPadPitchWidthLow();
 383   else
 384     padw = fTransform->GetPadPitchWidthUp();
 385   Float_t charge = (Float_t)cl->fCharge;
 386   Float_t pad = (Float_t)cl->fPad/charge;
 387   Float_t time = (Float_t)cl->fTime/charge;
 388   Float_t s2 = (Float_t)cl->fSigmaY2/charge - pad*pad;
 389   sigmaY2 = (s2 + 1./12)*padw*padw;
 390
 391   if(s2 != 0)
 392     {
 393       sigmaY2 = sigmaY2*0.108;
 394       if(fPatch<3)
 395         sigmaY2 = sigmaY2*2.07;
 396     }
 397
 398   s2 = (Float_t)cl->fSigmaZ2/charge - time*time;
 399   timew = fTransform->GetZWidth();
 400   sigmaZ2 = (s2 +1./12)*timew*timew;
 401   if(s2 != 0)
 402     {
 403       sigmaZ2 = sigmaZ2*0.169;
 404       if(fPatch < 3)
 405         sigmaZ2 = sigmaZ2*1.77;
 406     }
 407
 408 }