]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPC.cxx
Changes towards speeding up the code. Thanks to Marian Ivanov.
[u/mrichter/AliRoot.git] / TPC / AliTPC.cxx
index 08ffd8c3732507e0dedba69b1c7cc06f68a09337..a840e735f7b1255352a86c175588d033fddd9fea 100644 (file)
 
 /*
 $Log$
+Revision 1.52  2002/02/18 09:26:09  kowal2
+Removed compiler warning
+
+Revision 1.51  2002/01/21 17:13:21  kowal2
+New track hits using root containers. Setting active sectors added.
+
+Revision 1.50  2001/12/06 14:16:19  kowal2
+meaningfull printouts
+
+Revision 1.49  2001/11/30 11:55:37  hristov
+Noise table created in Hits2SDigits (M.Ivanov)
+
+Revision 1.48  2001/11/24 16:10:21  kowal2
+Faster algorithms.
+
+Revision 1.47  2001/11/19 10:25:34  kowal2
+Nearest integer instead of integer when converting to ADC counts
+
+Revision 1.46  2001/11/07 06:47:12  kowal2
+Removed printouts
+
+Revision 1.45  2001/11/03 13:33:48  kowal2
+Updated algorithms in Hits2SDigits, SDigits2Digits,
+Hits2ExactClusters.
+Added method Merge
+
 Revision 1.44  2001/08/30 09:28:48  hristov
 TTree names are explicitly set via SetName(name) and then Write() is called
 
@@ -175,6 +201,7 @@ Introduction of the Copyright and cvs Log
 #include "AliDigits.h"
 #include "AliSimDigits.h"
 #include "AliTPCTrackHits.h"
+#include "AliTPCTrackHitsV2.h"
 #include "AliPoints.h"
 #include "AliArrayBranch.h"
 
@@ -196,6 +223,31 @@ Introduction of the Copyright and cvs Log
 
 ClassImp(AliTPC) 
 
+//_____________________________________________________________________________
+// helper class for fast matrix and vector manipulation - no range checking
+// origin - Marian Ivanov
+
+class AliTPCFastMatrix : public TMatrix {
+public :
+  AliTPCFastMatrix(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb);
+  inline Float_t & UncheckedAt(Int_t rown, Int_t coln) const  {return  (fIndex[coln])[rown];} //fast acces   
+  inline Float_t   UncheckedAtFast(Int_t rown, Int_t coln) const  {return  (fIndex[coln])[rown];} //fast acces   
+};
+
+AliTPCFastMatrix::AliTPCFastMatrix(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb):
+  TMatrix(row_lwb, row_upb,col_lwb,col_upb)
+   {
+   };
+//_____________________________________________________________________________
+class AliTPCFastVector : public TVector {
+public :
+  AliTPCFastVector(Int_t size);
+  inline Float_t & UncheckedAt(Int_t index) const  {return  fElements[index];} //fast acces  
+};
+
+AliTPCFastVector::AliTPCFastVector(Int_t size):TVector(size){
+};
+
 //_____________________________________________________________________________
 AliTPC::AliTPC()
 {
@@ -206,13 +258,16 @@ AliTPC::AliTPC()
   fHits     = 0;
   fDigits   = 0;
   fNsectors = 0;
-  //MI changes
   fDigitsArray = 0;
   fClustersArray = 0;
   fDefaults = 0;
-  fTrackHits = 0;  
-  fHitType = 2;  
-  fTPCParam = 0; 
+  fTrackHits = 0; 
+  fTrackHitsOld = 0;   
+  fHitType = 4; //default CONTAINERS - based on ROOT structure 
+  fTPCParam = 0;    
+  fNoiseTable = 0;
+  fActiveSectors =0;
+
 }
  
 //_____________________________________________________________________________
@@ -226,18 +281,25 @@ AliTPC::AliTPC(const char *name, const char *title)
   //
   // Initialise arrays of hits and digits 
   fHits     = new TClonesArray("AliTPChit",  176);
-  gAlice->AddHitList(fHits);
-  //MI change  
+  gAlice->AddHitList(fHits); 
   fDigitsArray = 0;
   fClustersArray= 0;
   fDefaults = 0;
   //
-  fTrackHits = new AliTPCTrackHits;  //MI - 13.09.2000
+  fTrackHits = new AliTPCTrackHitsV2;  
   fTrackHits->SetHitPrecision(0.002);
   fTrackHits->SetStepPrecision(0.003);  
-  fTrackHits->SetMaxDistance(100); 
+  fTrackHits->SetMaxDistance(100);
 
-  fHitType = 2;
+  fTrackHitsOld = new AliTPCTrackHits;  //MI - 13.09.2000
+  fTrackHitsOld->SetHitPrecision(0.002);
+  fTrackHitsOld->SetStepPrecision(0.003);  
+  fTrackHitsOld->SetMaxDistance(100); 
+
+  fNoiseTable =0;
+
+  fHitType = 4;
+  fActiveSectors = 0;
   //
   // Initialise counters
   fNsectors = 0;
@@ -274,6 +336,9 @@ AliTPC::~AliTPC()
   delete fDigits;
   delete fTPCParam;
   delete fTrackHits; //MI 15.09.2000
+  delete fTrackHitsOld; //MI 10.12.2001
+  if (fNoiseTable) delete [] fNoiseTable;
+
 }
 
 //_____________________________________________________________________________
@@ -288,7 +353,7 @@ void AliTPC::AddHit(Int_t track, Int_t *vol, Float_t *hits)
     TClonesArray &lhits = *fHits;
     new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits);
   }
-  if (fHitType&2)
+  if (fHitType>1)
    AddHit2(track,vol,hits);
 }
  
@@ -803,6 +868,97 @@ void AliTPC::CreateMaterials()
     
 }
 
+void AliTPC::GenerNoise(Int_t tablesize)
+{
+  //
+  //Generate table with noise
+  //
+  if (fTPCParam==0) {
+    // error message
+    fNoiseDepth=0;
+    return;
+  }
+  if (fNoiseTable)  delete[] fNoiseTable;
+  fNoiseTable = new Float_t[tablesize];
+  fNoiseDepth = tablesize; 
+  fCurrentNoise =0; //!index of the noise in  the noise table 
+  
+  Float_t norm = fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac();
+  for (Int_t i=0;i<tablesize;i++) fNoiseTable[i]= gRandom->Gaus(0,norm);      
+}
+
+Float_t AliTPC::GetNoise()
+{
+  // get noise from table
+  //  if ((fCurrentNoise%10)==0) 
+  //  fCurrentNoise= gRandom->Rndm()*fNoiseDepth;
+  if (fCurrentNoise>=fNoiseDepth) fCurrentNoise=0;
+  return fNoiseTable[fCurrentNoise++];
+  //gRandom->Gaus(0, fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac()); 
+}
+
+
+Bool_t  AliTPC::IsSectorActive(Int_t sec)
+{
+  //
+  // check if the sector is active
+  if (!fActiveSectors) return kTRUE;
+  else return fActiveSectors[sec]; 
+}
+
+void    AliTPC::SetActiveSectors(Int_t * sectors, Int_t n)
+{
+  // activate interesting sectors
+  if (fActiveSectors) delete [] fActiveSectors;
+  fActiveSectors = new Bool_t[fTPCParam->GetNSector()];
+  for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kFALSE;
+  for (Int_t i=0;i<n;i++) 
+    if ((sectors[i]>=0) && sectors[i]<fTPCParam->GetNSector())  fActiveSectors[sectors[i]]=kTRUE;
+    
+}
+
+void    AliTPC::SetActiveSectors()
+{
+  //
+  // activate sectors which were hitted by tracks 
+  //loop over tracks
+  if (fHitType==0) return;  // if Clones hit - not short volume ID information
+  if (fActiveSectors) delete [] fActiveSectors;
+  fActiveSectors = new Bool_t[fTPCParam->GetNSector()];
+  for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kFALSE;
+  TBranch * branch=0;
+  if (fHitType>1) branch = gAlice->TreeH()->GetBranch("TPC2");
+  else branch = gAlice->TreeH()->GetBranch("TPC");
+  Stat_t ntracks = gAlice->TreeH()->GetEntries();
+  // loop over all hits
+  for(Int_t track=0;track<ntracks;track++){
+    ResetHits();
+    //
+    if (fTrackHits && fHitType&4) {
+      TBranch * br1 = gAlice->TreeH()->GetBranch("fVolumes");
+      TBranch * br2 = gAlice->TreeH()->GetBranch("fNVolumes");    
+      br1->GetEvent(track);
+      br2->GetEvent(track);
+      Int_t *volumes = fTrackHits->GetVolumes();
+      for (Int_t j=0;j<fTrackHits->GetNVolumes(); j++)
+       fActiveSectors[volumes[j]]=kTRUE;
+    }
+    
+    //
+    if (fTrackHitsOld && fHitType&2) {
+      TBranch * br = gAlice->TreeH()->GetBranch("fTrackHitsInfo");
+      br->GetEvent(track);
+      AliObjectArray * ar = fTrackHitsOld->fTrackHitsInfo;
+      for (UInt_t j=0;j<ar->GetSize();j++){
+       fActiveSectors[((AliTrackHitsInfo*)ar->At(j))->fVolumeID] =kTRUE;
+      } 
+    }    
+  }
+  
+}  
+
+
+
 
 void AliTPC::Digits2Clusters(TFile *of, Int_t eventnumber)
 {
@@ -910,17 +1066,11 @@ void AliTPC::Hits2Clusters(TFile *of, Int_t eventn)
       tH->GetEvent(track);
       //
       //  Get number of the TPC hits
-      //
-      // nhits=fHits->GetEntriesFast();
-      //
-     
+      //     
        tpcHit = (AliTPChit*)FirstHit(-1);
 
       // Loop over hits
       //
-       //   for(Int_t hit=0;hit<nhits;hit++){
-       //tpcHit=(AliTPChit*)fHits->UncheckedAt(hit);
-
        while(tpcHit){
  
         if (tpcHit->fQ == 0.) {
@@ -929,9 +1079,6 @@ void AliTPC::Hits2Clusters(TFile *of, Int_t eventn)
         }
        sector=tpcHit->fSector; // sector number
 
-
-       //      if(sector != isec) continue; //terminate iteration
-
        if(sector != isec){
         tpcHit = (AliTPChit*) NextHit();
         continue; 
@@ -1008,7 +1155,6 @@ void AliTPC::Hits2Clusters(TFile *of, Int_t eventn)
   } // end of loop over sectors  
 
   cerr<<"Number of made clusters : "<<nclusters<<"                        \n";
-
   carray.GetTree()->SetName(cname);
   carray.GetTree()->Write();
 
@@ -1037,8 +1183,8 @@ void AliTPC::Hits2ExactClustersSector(Int_t isec)
   //  Int_t sector,nhits;
   Int_t ipart;
   const Int_t kcmaxhits=30000;
-  TVector * xxxx = new TVector(kcmaxhits*4);
-  TVector & xxx = *xxxx;
+  AliTPCFastVector * xxxx = new AliTPCFastVector(kcmaxhits*4);
+  AliTPCFastVector & xxx = *xxxx;
   Int_t maxhits = kcmaxhits;
   //construct array for each padrow
   for (Int_t i=0; i<fTPCParam->GetNRow(isec);i++) 
@@ -1051,39 +1197,25 @@ void AliTPC::Hits2ExactClustersSector(Int_t isec)
   TTree *tH = gAlice->TreeH();
   Stat_t ntracks = tH->GetEntries();
   Int_t npart = gAlice->GetNtrack();
-  TBranch * br = tH->GetBranch("fTrackHitsInfo");
   //MI change
   TBranch * branch=0;
-  if (fHitType&2) branch = tH->GetBranch("TPC2");
+  if (fHitType>1) branch = tH->GetBranch("TPC2");
   else branch = tH->GetBranch("TPC");
 
   //------------------------------------------------------------
   // Loop over tracks
   //------------------------------------------------------------
-  
+
   for(Int_t track=0;track<ntracks;track++){ 
     Bool_t isInSector=kTRUE;
     ResetHits();
-        
-    if (fHitType&2) {
-      isInSector=kFALSE;    
-      br->GetEvent(track);
-      AliObjectArray * ar = fTrackHits->fTrackHitsInfo;
-      for (UInt_t j=0;j<ar->GetSize();j++){
-       if (  ((AliTrackHitsInfo*)ar->At(j))->fVolumeID==isec) isInSector=kTRUE;
-      }
-    }
+     isInSector = TrackInVolume(isec,track);
     if (!isInSector) continue;
     //MI change
     branch->GetEntry(track); // get next track
-
-
     //
     //  Get number of the TPC hits and a pointer
     //  to the particles
-    //
-    //nhits=fHits->GetEntriesFast();
-    //
     // Loop over hits
     //
     Int_t currentIndex=0;
@@ -1100,11 +1232,6 @@ void AliTPC::Hits2ExactClustersSector(Int_t isec)
        continue; 
       }
 
-      //for(Int_t hit=0;hit<nhits;hit++){
-      //tpcHit=(AliTPChit*)fHits->UncheckedAt(hit);
-      //if (tpcHit==0) continue;
-      //sector=tpcHit->fSector; // sector number
-      //if(sector != isec) continue; 
       ipart=tpcHit->Track();
       if (ipart<npart) particle=gAlice->Particle(ipart);
       
@@ -1212,7 +1339,7 @@ void AliTPC::Hits2ExactClustersSector(Int_t isec)
 void AliTPC::SDigits2Digits2(Int_t eventnumber)  
 {
   //create digits from summable digits
-
+  GenerNoise(500000); //create teble with noise
   char  sname[100];
   char  dname[100];
   sprintf(sname,"TreeS_%s_%d",fTPCParam->GetTitle(),eventnumber);
@@ -1250,6 +1377,7 @@ void AliTPC::SDigits2Digits2(Int_t eventnumber)
       cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
       continue;
     }
+    if (!IsSectorActive(sec)) continue;
     AliSimDigits * digrow =(AliSimDigits*) arr->CreateRow(sec,row);
     Int_t nrows = digrow->GetNRows();
     Int_t ncols = digrow->GetNCols();
@@ -1260,26 +1388,41 @@ void AliTPC::SDigits2Digits2(Int_t eventnumber)
     digarr.ExpandTrackBuffer();
 
     
-    for (Int_t rows=0;rows<nrows; rows++){
-      for (Int_t col=0;col<ncols; col++){
-       Float_t  q  = digarr.GetDigitFast(rows,col);
-       //add noise
-       q = gRandom->Gaus(q,fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac()*16); 
-       q/=16;  //conversion faktor
-        q=(Short_t)q; // digits are short integers
-
-       if (q > zerosup){
-          if(q > fTPCParam->GetADCSat()) q = (Short_t)(fTPCParam->GetADCSat());
-         digrow->SetDigitFast(q,rows,col);  
-         for (Int_t tr=0;tr<3;tr++)    
-          
-((AliSimDigits*)digrow)->SetTrackIDFast(digarr.GetTrackIDFast(rows,col,tr)-2,rows,col,tr);
-       }
-       }   
-      }        
+    Short_t * pamp0 = digarr.GetDigits();
+    Int_t   * ptracks0 = digarr.GetTracks();
+    Short_t * pamp1 = digrow->GetDigits();
+    Int_t   * ptracks1 = digrow->GetTracks();
+    Int_t  nelems =nrows*ncols;
+    Int_t saturation = fTPCParam->GetADCSat();
+    //use internal structure of the AliDigits - for speed reason
+    //if you cahnge implementation
+    //of the Alidigits - it must be rewriten -
+    for (Int_t i= 0; i<nelems; i++){
+      //      Float_t q = *pamp0;
+      //q/=16.;  //conversion faktor
+      //Float_t noise= GetNoise(); 
+      //q+=noise;      
+      //q= TMath::Nint(q);
+      Float_t q = TMath::Nint(Float_t(*pamp0)/16.+GetNoise());
+      if (q>zerosup){
+       if (q>saturation) q=saturation;      
+       *pamp1=(Short_t)q;
+       //if (ptracks0[0]==0)
+       //  ptracks1[0]=1;
+       //else
+       ptracks1[0]=ptracks0[0];        
+       ptracks1[nelems]=ptracks0[nelems];
+       ptracks1[2*nelems]=ptracks0[2*nelems];
+      }
+      pamp0++;
+      pamp1++;
+      ptracks0++;
+      ptracks1++;       
+    }
+
     arr->StoreRow(sec,row);
     arr->ClearRow(sec,row);   
-    cerr<<sec<<"\t"<<row<<"\n";   
+    // cerr<<sec<<"\t"<<row<<"\n";   
   }  
 
     
@@ -1375,8 +1518,8 @@ void AliTPC::Merge(TTree * intree, Int_t *mask, Int_t nin, Int_t outid)
 
         if (q> zerosup){
 
-          if(q > fTPCParam->GetADCSat()) q = (Short_t)(fTPCParam->GetADCSat());
-          digrow->SetDigitFast(q,rows,col);  
+          if(q > fTPCParam->GetADCSat()) q = fTPCParam->GetADCSat();
+          digrow->SetDigitFast((Short_t)q,rows,col);  
           for (Int_t tr=0;tr<3;tr++){
             if (tr<labptr)
               ((AliSimDigits*)digrow)->SetTrackIDFast(label[tr],rows,col,tr);
@@ -1389,7 +1532,7 @@ void AliTPC::Merge(TTree * intree, Int_t *mask, Int_t nin, Int_t outid)
       arr->StoreRow(sec,row);
 
       arr->ClearRow(sec,row);   
-      cerr<<sec<<"\t"<<row<<"\n";  
   }  
   
   delete digarr;
@@ -1467,6 +1610,7 @@ void AliTPC::Hits2Digits(Int_t eventnumber)
 
 
   if(fDefaults == 0) SetDefaults();  // check if the parameters are set
+  GenerNoise(500000); //create teble with noise
 
   //setup TPCDigitsArray 
 
@@ -1480,9 +1624,9 @@ void AliTPC::Hits2Digits(Int_t eventnumber)
 
   fDigitsSwitch=0; // standard digits
 
-  cerr<<"Digitizing TPC...\n";
+  cerr<<"Digitizing TPC -- normal digits...\n";
 
- for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) Hits2DigitsSector(isec);
+ for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) if (IsSectorActive(isec)) Hits2DigitsSector(isec);
 
   // write results
 
@@ -1512,7 +1656,7 @@ void AliTPC::Hits2SDigits2(Int_t eventnumber)
 
 
   if(fDefaults == 0) SetDefaults();
-
+  GenerNoise(500000); //create table with noise
   //setup TPCDigitsArray 
 
   if(GetDigitsArray()) delete GetDigitsArray();
@@ -1523,7 +1667,7 @@ void AliTPC::Hits2SDigits2(Int_t eventnumber)
   arr->MakeTree(fDigitsFile);
   SetDigitsArray(arr);
 
-  cerr<<"Digitizing TPC...\n"; 
+  cerr<<"Digitizing TPC -- summable digits...\n"; 
 
   fDigitsSwitch=1; // summable digits
   
@@ -1531,7 +1675,7 @@ void AliTPC::Hits2SDigits2(Int_t eventnumber)
 
   fTPCParam->SetZeroSup(0);
 
- for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) Hits2DigitsSector(isec);
+ for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) if (IsSectorActive(isec)) Hits2DigitsSector(isec);
 
 
   // write results
@@ -1562,6 +1706,7 @@ void AliTPC::Hits2SDigits()
   Int_t eventnumber = gAlice->GetEvNumber();
 
   if(fDefaults == 0) SetDefaults();
+  GenerNoise(500000); //create table with noise
 
   //setup TPCDigitsArray 
 
@@ -1573,11 +1718,11 @@ void AliTPC::Hits2SDigits()
   arr->MakeTree(fDigitsFile);
   SetDigitsArray(arr);
 
-  cerr<<"Digitizing TPC...\n"; 
+  cerr<<"Digitizing TPC -- summable digits...\n"; 
 
   //  fDigitsSwitch=1; // summable digits  -for the moment direct
 
for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) Hits2DigitsSector(isec);
 for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) if (IsSectorActive(isec)) Hits2DigitsSector(isec);
 
 
   // write results
@@ -1632,7 +1777,7 @@ void AliTPC::Hits2DigitsSector(Int_t isec)
 
       //--------------------------------------------------------
       //   Digitize this sector, row by row
-      //   row[i] is the pointer to the TObjArray of TVectors,
+      //   row[i] is the pointer to the TObjArray of AliTPCFastVectors,
       //   each one containing electrons accepted on this
       //   row, assigned into tracks
       //--------------------------------------------------------
@@ -1649,11 +1794,9 @@ void AliTPC::Hits2DigitsSector(Int_t isec)
 
        fDigitsArray->StoreRow(isec,i);
 
-               Int_t ndig = dig->GetDigitSize(); 
-        
-       //printf("*** Sector, row, compressed digits %d %d %d ***\n",isec,i,ndig);
-       
+       Int_t ndig = dig->GetDigitSize(); 
+       if (gDebug > 10) printf("*** Sector, row, compressed digits %d %d %d ***\n",isec,i,ndig);        
+       
         fDigitsArray->ClearRow(isec,i);  
 
    
@@ -1697,10 +1840,11 @@ void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows)
   //  Integrated signal for this row
   //  and a single track signal
   //    
-  TMatrix *m1   = new TMatrix(0,nofPads,0,nofTbins); // integrated
-  TMatrix *m2   = new TMatrix(0,nofPads,0,nofTbins); // single
+
+  AliTPCFastMatrix *m1 = new AliTPCFastMatrix(0,nofPads,0,nofTbins); // integrated
+  AliTPCFastMatrix *m2 = new AliTPCFastMatrix(0,nofPads,0,nofTbins); // single
   //
-  TMatrix &total  = *m1;
+  AliTPCFastMatrix &total  = *m1;
 
   //  Array of pointers to the label-signal list
 
@@ -1732,28 +1876,27 @@ void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows)
   Int_t tracks[3];
 
   AliDigits *dig = fDigitsArray->GetRow(isec,irow);
-  for(Int_t ip=0;ip<nofPads;ip++){
-    for(Int_t it=0;it<nofTbins;it++){
-
-      Float_t q = total(ip,it);
-
-      Int_t gi =it*nofPads+ip; // global index
-
+  Int_t gi=-1;
+  Float_t fzerosup = zerosup+0.5;
+  for(Int_t it=0;it<nofTbins;it++){
+    Float_t *pq = &(total.UncheckedAt(0,it));
+    for(Int_t ip=0;ip<nofPads;ip++){
+      gi++;
+      Float_t q=*pq;      
+      pq++;
       if(fDigitsSwitch == 0){
-
-        q = gRandom->Gaus(q,fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac()); 
-
-        q = (Int_t)q;
-
-        if(q <=zerosup) continue; // do not fill zeros
+       q+=GetNoise();
+        if(q <=fzerosup) continue; // do not fill zeros
+        q = TMath::Nint(q);
         if(q > fTPCParam->GetADCSat()) q = fTPCParam->GetADCSat();  // saturation
 
       }
 
       else {
        if(q <= 0.) continue; // do not fill zeros
+       if(q>2000.) q=2000.;
        q *= 16.;
-       q = (Int_t)q;
+       q = TMath::Nint(q);
       }
 
       //
@@ -1761,7 +1904,7 @@ void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows)
       //    
 
       for(Int_t j1=0;j1<3;j1++){
-        tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -1;
+       tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -2;
       }
 
 //Begin_Html
@@ -1800,8 +1943,8 @@ void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows)
 
 //_____________________________________________________________________________
 
-Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr, TMatrix *m1, TMatrix *m2,
-                          Int_t *indexRange)
+Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr, 
+             AliTPCFastMatrix *m1, AliTPCFastMatrix *m2,Int_t *indexRange)
 {
 
   //---------------------------------------------------------------
@@ -1815,10 +1958,10 @@ Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr, TMatrix *m1, TMatrix *m2,
   // Modified: Marian Ivanov 
   //-----------------------------------------------------------------
 
-  TVector *tv;
-  tv = (TVector*)p1->At(ntr); // pointer to a track
-  TVector &v = *tv;
+  AliTPCFastVector *tv;
+
+  tv = (AliTPCFastVector*)p1->At(ntr); // pointer to a track
+  AliTPCFastVector &v = *tv;
   
   Float_t label = v(0);
   Int_t centralPad = (fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3])-1)/2;
@@ -1829,10 +1972,8 @@ Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr, TMatrix *m1, TMatrix *m2,
   indexRange[2]=9999; //min time
   indexRange[3]=-1; // max time
 
-  //  Float_t IneffFactor = 0.5; // inefficiency in the gain close to the edge, as above
-
-  TMatrix &signal = *m1;
-  TMatrix &total = *m2;
+  AliTPCFastMatrix &signal = *m1;
+  AliTPCFastMatrix &total = *m2;
   //
   //  Loop over all electrons
   //
@@ -1841,22 +1982,28 @@ Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr, TMatrix *m1, TMatrix *m2,
     Float_t aval =  v(idx+4);
     Float_t eltoadcfac=aval*fTPCParam->GetTotalNormFac(); 
     Float_t xyz[3]={v(idx+1),v(idx+2),v(idx+3)};
-    Int_t n = fTPCParam->CalcResponse(xyz,fCurrentIndex,fCurrentIndex[3]);
-    
-    if (n>0) for (Int_t i =0; i<n; i++){
-       Int_t *index = fTPCParam->GetResBin(i);        
+    Int_t n = ((AliTPCParamSR*)fTPCParam)->CalcResponseFast(xyz,fCurrentIndex,fCurrentIndex[3]);
+
+    Int_t *index = fTPCParam->GetResBin(0);  
+    Float_t *weight = & (fTPCParam->GetResWeight(0));
+
+    if (n>0) for (Int_t i =0; i<n; i++){       
        Int_t pad=index[1]+centralPad;  //in digit coordinates central pad has coordinate 0
-       if ( ( pad<(fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3]))) && (pad>0)) {
+
+         if (pad>=0){
         Int_t time=index[2];    
-        Float_t weight = fTPCParam->GetResWeight(i); //we normalise response to ADC channel
-        weight *= eltoadcfac;
+         Float_t qweight = *(weight)*eltoadcfac;
         
-        if (m1!=0) signal(pad,time)+=weight; 
-        total(pad,time)+=weight;
-        indexRange[0]=TMath::Min(indexRange[0],pad);
-        indexRange[1]=TMath::Max(indexRange[1],pad);
-        indexRange[2]=TMath::Min(indexRange[2],time);
-        indexRange[3]=TMath::Max(indexRange[3],time); 
+        if (m1!=0) signal.UncheckedAt(pad,time)+=qweight;
+         total.UncheckedAt(pad,time)+=qweight;
+        if (indexRange[0]>pad) indexRange[0]=pad;
+        if (indexRange[1]<pad) indexRange[1]=pad;
+        if (indexRange[2]>time) indexRange[2]=time;
+        if (indexRange[3]<time) indexRange[3]=time;
+
+        index+=3;
+        weight++;      
+
        }        
     }
   } // end of loop over electrons
@@ -1865,8 +2012,8 @@ Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr, TMatrix *m1, TMatrix *m2,
 }
 
 //_____________________________________________________________________________
-void AliTPC::GetList(Float_t label,Int_t np,TMatrix *m,Int_t *indexRange,
-                     Float_t **pList)
+void AliTPC::GetList(Float_t label,Int_t np,AliTPCFastMatrix *m,
+                     Int_t *indexRange, Float_t **pList)
 {
   //----------------------------------------------------------------------
   //  Updates the list of tracks contributing to digits for a given row
@@ -1876,7 +2023,7 @@ void AliTPC::GetList(Float_t label,Int_t np,TMatrix *m,Int_t *indexRange,
   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
   //-----------------------------------------------------------------
 
-  TMatrix &signal = *m;
+  AliTPCFastMatrix &signal = *m;
 
   // lop over nonzero digits
 
@@ -1978,18 +2125,19 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
   AliTPChit *tpcHit; // pointer to a sigle TPC hit    
   //MI change
   TBranch * branch=0;
-  if (fHitType&2) branch = TH->GetBranch("TPC2");
+  if (fHitType>1) branch = TH->GetBranch("TPC2");
   else branch = TH->GetBranch("TPC");
 
  
   //----------------------------------------------
   // Create TObjArray-s, one for each row,
-  // each TObjArray will store the TVectors
-  // of electrons, one TVector per each track.
+  // each TObjArray will store the AliTPCFastVectors
+  // of electrons, one AliTPCFastVectors per each track.
   //---------------------------------------------- 
     
   Int_t *nofElectrons = new Int_t [nrows]; // electron counter for each row
-  TVector **tracks = new TVector* [nrows]; //pointers to the track vectors
+  AliTPCFastVector **tracks = new AliTPCFastVector* [nrows]; //pointers to the track vectors
+
   for(i=0; i<nrows; i++){
     row[i] = new TObjArray;
     nofElectrons[i]=0;
@@ -2008,16 +2156,7 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
   for(Int_t track=0;track<ntracks;track++){
     Bool_t isInSector=kTRUE;
     ResetHits();
-
-    if (fHitType&2) {
-      isInSector=kFALSE;
-      TBranch * br = TH->GetBranch("fTrackHitsInfo");
-      br->GetEvent(track);
-      AliObjectArray * ar = fTrackHits->fTrackHitsInfo;
-      for (UInt_t j=0;j<ar->GetSize();j++){
-       if (  ((AliTrackHitsInfo*)ar->At(j))->fVolumeID==isec) isInSector=kTRUE;
-      }
-    }
+    isInSector = TrackInVolume(isec,track);
     if (!isInSector) continue;
     //MI change
     branch->GetEntry(track); // get next track
@@ -2034,7 +2173,6 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
     while(tpcHit){
       
       Int_t sector=tpcHit->fSector; // sector number
-      //      if(sector != isec) continue; 
       if(sector != isec){
        tpcHit = (AliTPChit*) NextHit();
        continue; 
@@ -2050,19 +2188,19 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
           for(i=0;i<nrows;i++){
              if(previousTrack != -1){
               if(nofElectrons[i]>0){
-                TVector &v = *tracks[i];
+                 AliTPCFastVector &v = *tracks[i];
                 v(0) = previousTrack;
                  tracks[i]->ResizeTo(4*nofElectrons[i]+1); // shrink if necessary
                 row[i]->Add(tracks[i]);                     
               }
                else{
-                 delete tracks[i]; // delete empty TVector
+                 delete tracks[i]; // delete empty AliTPCFastVector
                  tracks[i]=0;
               }
             }
 
              nofElectrons[i]=0;
-             tracks[i] = new TVector(481); // TVectors for the next fTrack
+             tracks[i] = new AliTPCFastVector(481); // AliTPCFastVectors for the next fTrack
 
           } // end of loop over rows
               
@@ -2119,13 +2257,11 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
            }
          }
          
-         TVector &v = *tracks[rowNumber];
+          AliTPCFastVector &v = *tracks[rowNumber];
          Int_t idx = 4*nofElectrons[rowNumber]-3;
-
-         v(idx)=  xyz[0];   // X - pad row coordinate
-         v(idx+1)=xyz[1];   // Y - pad coordinate (along the pad-row)
-          v(idx+2)=xyz[2];   // Z - time bin coordinate
-         v(idx+3)=xyz[3];   // avalanche size  
+         Real_t * position = &(((AliTPCFastVector&)v).UncheckedAt(idx)); //make code faster
+         memcpy(position,xyz,4*sizeof(Float_t));
        } // end of loop over electrons
 
         tpcHit = (AliTPChit*)NextHit();
@@ -2139,7 +2275,7 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
 
      for(i=0;i<nrows;i++){
        if(nofElectrons[i]>0){
-          TVector &v = *tracks[i];
+          AliTPCFastVector &v = *tracks[i];
          v(0) = previousTrack;
           tracks[i]->ResizeTo(4*nofElectrons[i]+1); // shrink if necessary
          row[i]->Add(tracks[i]);  
@@ -2196,7 +2332,7 @@ void AliTPC::MakeBranch(Option_t* option, const char *file)
                        branchname, &fDigits, buffersize, file);
   }    
 
-  if (fHitType&2) MakeBranch2(option,file); // MI change 14.09.2000
+  if (fHitType>1) MakeBranch2(option,file); // MI change 14.09.2000
 }
  
 //_____________________________________________________________________________
@@ -2399,7 +2535,7 @@ AliHit(shunt,track)
  
 
 //________________________________________________________________________
-// Additional code because of the AliTPCTrackHits
+// Additional code because of the AliTPCTrackHitsV2
 
 void AliTPC::MakeBranch2(Option_t *option,const char *file)
 {
@@ -2407,17 +2543,21 @@ void AliTPC::MakeBranch2(Option_t *option,const char *file)
   // Create a new branch in the current Root Tree
   // The branch of fHits is automatically split
   // MI change 14.09.2000
-  if (fHitType&2==0) return;
+  if (fHitType<2) return;
   char branchname[10];
   sprintf(branchname,"%s2",GetName());  
   //
   // Get the pointer to the header
   const char *cH = strstr(option,"H");
   //
-  if (fTrackHits   && gAlice->TreeH() && cH) {    
-    AliObjectBranch * branch = new AliObjectBranch(branchname,"AliTPCTrackHits",&fTrackHits, 
-                                                  gAlice->TreeH(),fBufferSize,99);
-    gAlice->TreeH()->GetListOfBranches()->Add(branch);
+  if (fTrackHits   && gAlice->TreeH() && cH && fHitType&4) {    
+    //    AliObjectBranch * branch = new AliObjectBranch(branchname,"AliTPCTrackHitsV2",&fTrackHits, 
+    //                                            gAlice->TreeH(),fBufferSize,99);
+    //TBranch * branch = 
+    gAlice->TreeH()->Branch(branchname,"AliTPCTrackHitsV2",&fTrackHits, 
+                                                  fBufferSize,99);
+
+    // gAlice->TreeH()->GetListOfBranches()->Add(branch);
     if (GetDebug()>1) 
       printf("* AliDetector::MakeBranch * Making Branch %s for trackhits\n",branchname);
     const char folder [] = "RunMC/Event/Data";
@@ -2437,12 +2577,39 @@ void AliTPC::MakeBranch2(Option_t *option,const char *file)
              cout << "Diverting branch " << branchname << " to file " << file << endl;  
     }
   }    
+
+  if (fTrackHitsOld   && gAlice->TreeH() && cH && fHitType&2) {    
+    AliObjectBranch * branch = new AliObjectBranch(branchname,"AliTPCTrackHits",&fTrackHitsOld, 
+                                                  gAlice->TreeH(),fBufferSize,99);
+    //TBranch * branch = gAlice->TreeH()->Branch(branchname,"AliTPCTrackHitsV2",&fTrackHits, 
+    //                                                    fBufferSize,99);
+
+    gAlice->TreeH()->GetListOfBranches()->Add(branch);
+    if (GetDebug()>1) 
+      printf("* AliDetector::MakeBranch * Making Branch %s for trackhits\n",branchname);
+    const char folder [] = "RunMC/Event/Data";
+    if (GetDebug())
+      printf("%15s: Publishing %s to %s\n",ClassName(),branchname,folder);
+    Publish(folder,&fTrackHitsOld,branchname);
+    if (file) {
+        TBranch *b = gAlice->TreeH()->GetBranch(branchname);
+        TDirectory *wd = gDirectory;
+        b->SetFile(file);
+        TIter next( b->GetListOfBranches());
+        while ((b=(TBranch*)next())) {
+         b->SetFile(file);
+        }
+        wd->cd(); 
+        if (GetDebug()>1) 
+             cout << "Diverting branch " << branchname << " to file " << file << endl;  
+    }
+  }    
 }
 
 void AliTPC::SetTreeAddress()
 {
-  if (fHitType&1) AliDetector::SetTreeAddress();
-  if (fHitType&2) SetTreeAddress2();
+  if (fHitType<=1) AliDetector::SetTreeAddress();
+  if (fHitType>1) SetTreeAddress2();
 }
 
 void AliTPC::SetTreeAddress2()
@@ -2456,15 +2623,21 @@ void AliTPC::SetTreeAddress2()
   //
   // Branch address for hit tree
   TTree *treeH = gAlice->TreeH();
-  if (treeH) {
+  if ((treeH)&&(fHitType&4)) {
     branch = treeH->GetBranch(branchname);
     if (branch) branch->SetAddress(&fTrackHits);
   }
+  if ((treeH)&&(fHitType&2)) {
+    branch = treeH->GetBranch(branchname);
+    if (branch) branch->SetAddress(&fTrackHitsOld);
+  }
+
 }
 
 void AliTPC::FinishPrimary()
 {
-  if (fTrackHits) fTrackHits->FlushHitStack();  
+  if (fTrackHits &&fHitType&4)      fTrackHits->FlushHitStack();  
+  if (fTrackHitsOld && fHitType&2)  fTrackHitsOld->FlushHitStack();  
 }
 
 
@@ -2484,32 +2657,39 @@ void AliTPC::AddHit2(Int_t track, Int_t *vol, Float_t *hits)
   //AliTPChit *hit = (AliTPChit*)fHits->UncheckedAt(fNhits-1);
   //if (hit->fTrack!=rtrack)
   //  cout<<"bad track number\n";
-  if (fTrackHits) 
+  if (fTrackHits && fHitType&4
     fTrackHits->AddHitKartez(vol[0],rtrack, hits[0],
                             hits[1],hits[2],(Int_t)hits[3]);
+  if (fTrackHitsOld &&fHitType&2 ) 
+    fTrackHitsOld->AddHitKartez(vol[0],rtrack, hits[0],
+                            hits[1],hits[2],(Int_t)hits[3]);
+  
 }
 
 void AliTPC::ResetHits()
 {
   if (fHitType&1) AliDetector::ResetHits();
-  if (fHitType&2) ResetHits2();
+  if (fHitType>1) ResetHits2();
 }
 
 void AliTPC::ResetHits2()
 {
   //
   //reset hits
-  if (fTrackHits) fTrackHits->Clear();
+  if (fTrackHits && fHitType&4) fTrackHits->Clear();
+  if (fTrackHitsOld && fHitType&2) fTrackHitsOld->Clear();
+
 }   
 
 AliHit* AliTPC::FirstHit(Int_t track)
 {
-  if (fHitType&2) return FirstHit2(track);
+  if (fHitType>1) return FirstHit2(track);
   return AliDetector::FirstHit(track);
 }
 AliHit* AliTPC::NextHit()
 {
-  if (fHitType&2) return NextHit2();
+  if (fHitType>1) return NextHit2();
+  
   return AliDetector::NextHit();
 }
 
@@ -2527,10 +2707,15 @@ AliHit* AliTPC::FirstHit2(Int_t track)
     gAlice->TreeH()->GetEvent(track);
   }
   //
-  if (fTrackHits) {
+  if (fTrackHits && fHitType&4) {
     fTrackHits->First();
     return fTrackHits->GetHit();
   }
+  if (fTrackHitsOld && fHitType&2) {
+    fTrackHitsOld->First();
+    return fTrackHitsOld->GetHit();
+  }
+
   else return 0;
 }
 
@@ -2539,6 +2724,11 @@ AliHit* AliTPC::NextHit2()
   //
   //Return the next hit for the current track
 
+
+  if (fTrackHitsOld && fHitType&2) {
+    fTrackHitsOld->Next();
+    return fTrackHitsOld->GetHit();
+  }
   if (fTrackHits) {
     fTrackHits->Next();
     return fTrackHits->GetHit();
@@ -2565,14 +2755,64 @@ void AliTPC::LoadPoints(Int_t)
 void AliTPC::RemapTrackHitIDs(Int_t *map)
 {
   if (!fTrackHits) return;
-  AliObjectArray * arr = fTrackHits->fTrackHitsInfo;
-  for (UInt_t i=0;i<arr->GetSize();i++){
-    AliTrackHitsInfo * info = (AliTrackHitsInfo *)(arr->At(i));
-    info->fTrackID = map[info->fTrackID];
-  }
   
+  if (fTrackHitsOld && fHitType&2){
+    AliObjectArray * arr = fTrackHitsOld->fTrackHitsInfo;
+    for (UInt_t i=0;i<arr->GetSize();i++){
+      AliTrackHitsInfo * info = (AliTrackHitsInfo *)(arr->At(i));
+      info->fTrackID = map[info->fTrackID];
+    }
+  }
+  if (fTrackHitsOld && fHitType&4){
+    TClonesArray * arr = fTrackHits->GetArray();;
+    for (Int_t i=0;i<arr->GetEntriesFast();i++){
+      AliTrackHitsParamV2 * info = (AliTrackHitsParamV2 *)(arr->At(i));
+      info->fTrackID = map[info->fTrackID];
+    }
+  }
 }
 
+Bool_t   AliTPC::TrackInVolume(Int_t id,Int_t track)
+{
+  //return bool information - is track in given volume
+  //load only part of the track information 
+  //return true if current track is in volume
+  //
+  //  return kTRUE;
+  if (fTrackHitsOld && fHitType&2) {
+    TBranch * br = gAlice->TreeH()->GetBranch("fTrackHitsInfo");
+    br->GetEvent(track);
+    AliObjectArray * ar = fTrackHitsOld->fTrackHitsInfo;
+    for (UInt_t j=0;j<ar->GetSize();j++){
+      if (  ((AliTrackHitsInfo*)ar->At(j))->fVolumeID==id) return kTRUE;
+    } 
+  }
+
+  if (fTrackHits && fHitType&4) {
+    TBranch * br1 = gAlice->TreeH()->GetBranch("fVolumes");
+    TBranch * br2 = gAlice->TreeH()->GetBranch("fNVolumes");    
+    br2->GetEvent(track);
+    br1->GetEvent(track);    
+    Int_t *volumes = fTrackHits->GetVolumes();
+    Int_t nvolumes = fTrackHits->GetNVolumes();
+    if (!volumes && nvolumes>0) {
+      printf("Problematic track\t%d\t%d",track,nvolumes);
+      return kFALSE;
+    }
+    for (Int_t j=0;j<nvolumes; j++)
+      if (volumes[j]==id) return kTRUE;    
+  }
+
+  if (fHitType&1) {
+    TBranch * br = gAlice->TreeH()->GetBranch("fSector");
+    br->GetEvent(track);
+    for (Int_t j=0;j<fHits->GetEntriesFast();j++){
+      if (  ((AliTPChit*)fHits->At(j))->fSector==id) return kTRUE;
+    } 
+  }
+  return kFALSE;  
+
+}
 
 //_____________________________________________________________________________
 void AliTPC::LoadPoints2(Int_t)
@@ -2580,9 +2820,12 @@ void AliTPC::LoadPoints2(Int_t)
   //
   // Store x, y, z of all hits in memory
   //
-  if (fTrackHits == 0) return;
+  if (fTrackHits == 0 && fTrackHitsOld==0) return;
   //
-  Int_t nhits = fTrackHits->GetEntriesFast();
+  Int_t nhits =0;
+  if (fHitType&4) nhits = fTrackHits->GetEntriesFast();
+  if (fHitType&2) nhits = fTrackHitsOld->GetEntriesFast();
+  
   if (nhits == 0) return;
   Int_t tracks = gAlice->GetNtrack();
   if (fPoints == 0) fPoints = new TObjArray(tracks);
@@ -2605,9 +2848,7 @@ void AliTPC::LoadPoints2(Int_t)
   // Loop over all the hits and store their position
   //
   ahit = FirstHit2(-1);
-  //for (Int_t hit=0;hit<nhits;hit++) {
   while (ahit){
-    //    ahit = (AliHit*)fHits->UncheckedAt(hit);
     trk=ahit->GetTrack();
     if(ntrk[trk]==limi[trk]) {
       //
@@ -2628,6 +2869,9 @@ void AliTPC::LoadPoints2(Int_t)
     ntrk[trk]++;
     ahit = NextHit2();
   }
+
+
+
   //
   for(trk=0; trk<tracks; ++trk) {
     if(ntrk[trk]) {
@@ -2746,8 +2990,8 @@ void AliTPC::FindTrackHitsIntersection(TClonesArray * arr)
   Int_t ipart;
   
   const Int_t kcmaxhits=30000;
-  TVector * xxxx = new TVector(kcmaxhits*4);
-  TVector & xxx = *xxxx;
+  AliTPCFastVector * xxxx = new AliTPCFastVector(kcmaxhits*4);
+  AliTPCFastVector & xxx = *xxxx;
   Int_t maxhits = kcmaxhits;
       
   //