]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPC.cxx
Renamed to TFluka.
[u/mrichter/AliRoot.git] / TPC / AliTPC.cxx
index 8d5521191c40509098d4422ea9700dcf3dba2e27..6020cde38670a2e8879f30244c0f54c3ef5eb149 100644 (file)
 #include <TSystem.h>     
 #include <TTUBS.h>
 #include <TTree.h>
-#include <TVector.h>
 #include <TVirtualMC.h>
+#include <TString.h>
+#include <TF2.h>
+#include <TStopwatch.h>
 
 #include "AliArrayBranch.h"
-#include "AliClusters.h"
-#include "AliComplexCluster.h"
 #include "AliDigits.h"
 #include "AliMagF.h"
 #include "AliPoints.h"
@@ -63,8 +63,6 @@
 #include "AliSimDigits.h"
 #include "AliTPC.h"
 #include "AliTPC.h"
-#include "AliTPCClustersArray.h"
-#include "AliTPCClustersRow.h"
 #include "AliTPCDigitsArray.h"
 #include "AliTPCLoader.h"
 #include "AliTPCPRF2D.h"
 #include "AliTPCRF1D.h"
 #include "AliTPCTrackHits.h"
 #include "AliTPCTrackHitsV2.h"
-#include "AliTPCcluster.h"
 #include "AliTrackReference.h"
 #include "AliMC.h"
+#include "AliTPCDigitizer.h"
+#include "AliTPCBuffer.h"
+#include "AliTPCDDLRawData.h"
 
 
 ClassImp(AliTPC) 
@@ -85,27 +85,20 @@ ClassImp(AliTPC)
 
 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(Int_t rowlwb, Int_t rowupb, Int_t collwb, Int_t colupb);
+#if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,1)
+  Float_t & UncheckedAt(Int_t rown, Int_t coln) const  {return fElements[(rown-fRowLwb)*fNcols+(coln-fColLwb)];} //fast acces
+  Float_t   UncheckedAtFast(Int_t rown, Int_t coln) const  {return fElements[(rown-fRowLwb)*fNcols+(coln-fColLwb)];} //fast acces
+#else
+  Float_t & UncheckedAt(Int_t rown, Int_t coln) const  {return  (fIndex[coln])[rown];} //fast acces   
+  Float_t   UncheckedAtFast(Int_t rown, Int_t coln) const  {return  (fIndex[coln])[rown];} //fast acces   
+#endif
 };
 
-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)
+AliTPCFastMatrix::AliTPCFastMatrix(Int_t rowlwb, Int_t rowupb, Int_t collwb, Int_t colupb):
+  TMatrix(rowlwb, rowupb,collwb,colupb)
    {
    };
-//_____________________________________________________________________________
-class AliTPCFastVector : public TVector {
-public :
-  AliTPCFastVector(Int_t size);
-  virtual ~AliTPCFastVector(){;}
-  inline Float_t & UncheckedAt(Int_t index) const  {return  fElements[index];} //fast acces  
-  
-};
-
-AliTPCFastVector::AliTPCFastVector(Int_t size):TVector(size){
-};
-
 //_____________________________________________________________________________
 AliTPC::AliTPC()
 {
@@ -117,7 +110,6 @@ AliTPC::AliTPC()
   fDigits   = 0;
   fNsectors = 0;
   fDigitsArray = 0;
-  fClustersArray = 0;
   fDefaults = 0;
   fTrackHits = 0; 
   fTrackHitsOld = 0;   
@@ -141,7 +133,6 @@ AliTPC::AliTPC(const char *name, const char *title)
   fHits     = new TClonesArray("AliTPChit",  176);
   gAlice->GetMCApp()->AddHitList(fHits); 
   fDigitsArray = 0;
-  fClustersArray= 0;
   fDefaults = 0;
   //
   fTrackHits = new AliTPCTrackHitsV2;  
@@ -183,6 +174,11 @@ AliTPC::AliTPC(const char *name, const char *title)
 }
 
 //_____________________________________________________________________________
+AliTPC::AliTPC(const AliTPC& t):AliDetector(t){
+  //
+  // dummy copy constructor
+  //
+}
 AliTPC::~AliTPC()
 {
   //
@@ -301,7 +297,7 @@ void AliTPC::BuildGeometry()
 }    
 
 //_____________________________________________________________________________
-Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t )
+Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t ) const
 {
   //
   // Calculate distance from TPC to mouse on the display
@@ -310,14 +306,6 @@ Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t )
   return 9999;
 }
 
-void AliTPC::Clusters2Tracks() 
- {
-  //-----------------------------------------------------------------
-  // This is a track finder.
-  //-----------------------------------------------------------------
-  Error("Clusters2Tracks",
-  "Dummy function !  Call AliTPCtracker::Clusters2Tracks(...) instead !");
- }
 
 //_____________________________________________________________________________
 void AliTPC::CreateMaterials()
@@ -435,7 +423,7 @@ void AliTPC::CreateMaterials()
   density = 0.;
   Float_t am=0;
   Int_t nc;
-  Float_t rho,absl,X0,buf[1];
+  Float_t rho,absl,x0,buf[1];
   Int_t nbuf;
   Float_t a,z;
 
@@ -444,7 +432,7 @@ void AliTPC::CreateMaterials()
     
       // retrive material constants
       
-      gMC->Gfmate((*fIdmate)[fMixtComp[nc]],namate,a,z,rho,X0,absl,buf,nbuf);
+      gMC->Gfmate((*fIdmate)[fMixtComp[nc]],namate,a,z,rho,x0,absl,buf,nbuf);
 
       amat[nc] = a;
       zmat[nc] = z;
@@ -674,7 +662,7 @@ void AliTPC::CreateMaterials()
 
   // get epoxy
 
-  gMC->Gfmate((*fIdmate)[45],namate,amat[1],zmat[1],rho,X0,absl,buf,nbuf);
+  gMC->Gfmate((*fIdmate)[45],namate,amat[1],zmat[1],rho,x0,absl,buf,nbuf);
 
   // Carbon fiber
 
@@ -687,7 +675,7 @@ void AliTPC::CreateMaterials()
 
   // get SiO2
 
-  gMC->Gfmate((*fIdmate)[38],namate,amat[0],zmat[0],rho,X0,absl,buf,nbuf); 
+  gMC->Gfmate((*fIdmate)[38],namate,amat[0],zmat[0],rho,x0,absl,buf,nbuf); 
 
   wmat[0]=0.725; // by weight!
   wmat[1]=0.275;
@@ -757,7 +745,7 @@ Float_t AliTPC::GetNoise()
 }
 
 
-Bool_t  AliTPC::IsSectorActive(Int_t sec)
+Bool_t  AliTPC::IsSectorActive(Int_t sec) const
 {
   //
   // check if the sector is active
@@ -801,7 +789,7 @@ void    AliTPC::SetActiveSectors(Int_t flag)
   else branch = TreeH()->GetBranch("TPC");
   Stat_t ntracks = TreeH()->GetEntries();
   // loop over all hits
-  cout<<"\nAliTPC::SetActiveSectors():  Got "<<ntracks<<" tracks\n";
+  if (GetDebug()) cout<<"\nAliTPC::SetActiveSectors():  Got "<<ntracks<<" tracks\n";
   
   for(Int_t track=0;track<ntracks;track++)
    {
@@ -827,406 +815,117 @@ void    AliTPC::SetActiveSectors(Int_t flag)
       } 
     }    
   }
-  
 }  
 
 
 
 
-void AliTPC::Digits2Clusters(Int_t /*eventnumber*/)
-{
-  //-----------------------------------------------------------------
-  // This is a simple cluster finder.
-  //-----------------------------------------------------------------
-  Error("Digits2Clusters",
-  "Dummy function !  Call AliTPCclusterer::Digits2Clusters(...) instead !");
-}
-
-extern Double_t SigmaY2(Double_t, Double_t, Double_t);
-extern Double_t SigmaZ2(Double_t, Double_t);
 //_____________________________________________________________________________
-void AliTPC::Hits2Clusters(Int_t /*eventn*/)
+void AliTPC::Digits2Raw()
 {
-  //--------------------------------------------------------
-  // TPC simple cluster generator from hits
-  // obtained from the TPC Fast Simulator
-  // The point errors are taken from the parametrization
-  //--------------------------------------------------------
+// convert digits of the current event to raw data
 
-  //-----------------------------------------------------------------
-  // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
-  //-----------------------------------------------------------------
-  // Adopted to Marian's cluster data structure by I.Belikov, CERN,
-  // Jouri.Belikov@cern.ch
-  //----------------------------------------------------------------
-  
-  /////////////////////////////////////////////////////////////////////////////
-  //
-  //---------------------------------------------------------------------
-  //   ALICE TPC Cluster Parameters
-  //--------------------------------------------------------------------
-       
-  
+  static const Int_t kThreshold = 0;
+  static const Bool_t kCompress = kTRUE;
 
-  // Cluster width in rphi
-  const Float_t kACrphi=0.18322;
-  const Float_t kBCrphi=0.59551e-3;
-  const Float_t kCCrphi=0.60952e-1;
-  // Cluster width in z
-  const Float_t kACz=0.19081;
-  const Float_t kBCz=0.55938e-3;
-  const Float_t kCCz=0.30428;
-
-
-  if (!fLoader) {
-     cerr<<"AliTPC::Hits2Clusters(): output file not open !\n";
-     return;
+  fLoader->LoadDigits();
+  TTree* digits = fLoader->TreeD();
+  if (!digits) {
+    Error("Digits2Raw", "no digits tree");
+    return;
   }
 
-  //if(fDefaults == 0) SetDefaults();
-
-  Float_t sigmaRphi,sigmaZ,clRphi,clZ;
-  //
-  TParticle *particle; // pointer to a given particle
-  AliTPChit *tpcHit; // pointer to a sigle TPC hit
-  Int_t sector;
-  Int_t ipart;
-  Float_t xyz[5];
-  Float_t pl,pt,tanth,rpad,ratio;
-  Float_t cph,sph;
-  
-  //---------------------------------------------------------------
-  //  Get the access to the tracks 
-  //---------------------------------------------------------------
-  
-  TTree *tH = TreeH();
-  if (tH == 0x0)
-   {
-     Fatal("Hits2Clusters","Can not find TreeH in folder");
-     return;
-   }
-  SetTreeAddress();
-  
-  Stat_t ntracks = tH->GetEntries();
-
-  //Switch to the output file
-  
-  if (fLoader->TreeR() == 0x0) fLoader->MakeTree("R");
-  
-  cout<<"fTPCParam->GetTitle() = "<<fTPCParam->GetTitle()<<endl;
-  
-  AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::fgkRunLoaderName);
-  rl->CdGAFile();
-  //fTPCParam->Write(fTPCParam->GetTitle());
-
-  AliTPCClustersArray carray;
-  carray.Setup(fTPCParam);
-  carray.SetClusterType("AliTPCcluster");
-  carray.MakeTree(fLoader->TreeR());
-
-  Int_t nclusters=0; //cluster counter
-  
-  //------------------------------------------------------------
-  // Loop over all sectors (72 sectors for 20 deg
-  // segmentation for both lower and upper sectors)
-  // Sectors 0-35 are lower sectors, 0-17 z>0, 17-35 z<0
-  // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
-  //
-  // First cluster for sector 0 starts at "0"
-  //------------------------------------------------------------
-   
-  for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++){
-    //MI change
-    fTPCParam->AdjustCosSin(isec,cph,sph);
-    
-    //------------------------------------------------------------
-    // Loop over tracks
-    //------------------------------------------------------------
-    
-    for(Int_t track=0;track<ntracks;track++){
-      ResetHits();
-      tH->GetEvent(track);
-      //
-      //  Get number of the TPC hits
-      //     
-       tpcHit = (AliTPChit*)FirstHit(-1);
-
-      // Loop over hits
-      //
-       while(tpcHit){
-        if (tpcHit->fQ == 0.) {
-           tpcHit = (AliTPChit*) NextHit();
-           continue; //information about track (I.Belikov)
-        }
-       sector=tpcHit->fSector; // sector number
-
-       if(sector != isec){
-        tpcHit = (AliTPChit*) NextHit();
-        continue; 
-       }
-       ipart=tpcHit->Track();
-       particle=gAlice->GetMCApp()->Particle(ipart);
-       pl=particle->Pz();
-       pt=particle->Pt();
-       if(pt < 1.e-9) pt=1.e-9;
-       tanth=pl/pt;
-       tanth = TMath::Abs(tanth);
-       rpad=TMath::Sqrt(tpcHit->X()*tpcHit->X() + tpcHit->Y()*tpcHit->Y());
-       ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
-
-       //   space-point resolutions
-       
-       sigmaRphi=SigmaY2(rpad,tanth,pt);
-       sigmaZ   =SigmaZ2(rpad,tanth   );
-       
-       //   cluster widths
-       
-       clRphi=kACrphi-kBCrphi*rpad*tanth+kCCrphi*ratio*ratio;
-       clZ=kACz-kBCz*rpad*tanth+kCCz*tanth*tanth;
-       
-       // temporary protection
-       
-       if(sigmaRphi < 0.) sigmaRphi=0.4e-3;
-       if(sigmaZ < 0.) sigmaZ=0.4e-3;
-       if(clRphi < 0.) clRphi=2.5e-3;
-       if(clZ < 0.) clZ=2.5e-5;
-       
-       //
-       
-       //
-       // smearing --> rotate to the 1 (13) or to the 25 (49) sector,
-       // then the inaccuracy in a X-Y plane is only along Y (pad row)!
-       //
-        Float_t xprim= tpcHit->X()*cph + tpcHit->Y()*sph;
-       Float_t yprim=-tpcHit->X()*sph + tpcHit->Y()*cph;
-       xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigmaRphi));   // y
-          Float_t alpha=(isec < fTPCParam->GetNInnerSector()) ?
-         fTPCParam->GetInnerAngle() : fTPCParam->GetOuterAngle();
-          Float_t ymax=xprim*TMath::Tan(0.5*alpha);
-          if (TMath::Abs(xyz[0])>ymax) xyz[0]=yprim; 
-       xyz[1]=gRandom->Gaus(tpcHit->Z(),TMath::Sqrt(sigmaZ)); // z
-          if (TMath::Abs(xyz[1])>fTPCParam->GetZLength()) xyz[1]=tpcHit->Z(); 
-       xyz[2]=sigmaRphi;                                     // fSigmaY2
-       xyz[3]=sigmaZ;                                        // fSigmaZ2
-       xyz[4]=tpcHit->fQ;                                    // q
-
-        AliTPCClustersRow *clrow=carray.GetRow(sector,tpcHit->fPadRow);
-        if (!clrow) clrow=carray.CreateRow(sector,tpcHit->fPadRow);    
-
-        Int_t tracks[3]={tpcHit->Track(), -1, -1};
-       AliTPCcluster cluster(tracks,xyz);
-
-        clrow->InsertCluster(&cluster); nclusters++;
-
-        tpcHit = (AliTPChit*)NextHit();
-        
-
-      } // end of loop over hits
-
-    }   // end of loop over tracks
-
-    Int_t nrows=fTPCParam->GetNRow(isec);
-    for (Int_t irow=0; irow<nrows; irow++) {
-        AliTPCClustersRow *clrow=carray.GetRow(isec,irow);
-        if (!clrow) continue;
-        carray.StoreRow(isec,irow);
-        carray.ClearRow(isec,irow);
+  AliSimDigits digarr;
+  AliSimDigits* digrow = &digarr;
+  digits->GetBranch("Segment")->SetAddress(&digrow);
+
+  const char* fileName = "AliTPCDDL.dat";
+  AliTPCBuffer* buffer  = new AliTPCBuffer(fileName);
+  //Verbose level
+  // 0: Silent
+  // 1: cout messages
+  // 2: txt files with digits 
+  //BE CAREFUL, verbose level 2 MUST be used only for debugging and
+  //it is highly suggested to use this mode only for debugging digits files
+  //reasonably small, because otherwise the size of the txt files can reach
+  //quickly several MB wasting time and disk space.
+  buffer->SetVerbose(0);
+
+  Int_t nEntries = Int_t(digits->GetEntries());
+  Int_t previousSector = -1;
+  Int_t subSector = 0;
+  for (Int_t i = 0; i < nEntries; i++) {
+    digits->GetEntry(i);
+    Int_t sector, row;
+    fTPCParam->AdjustSectorRow(digarr.GetID(), sector, row);
+    if(previousSector != sector) {
+      subSector = 0;
+      previousSector = sector;
     }
 
-  } // end of loop over sectors  
-
-  cerr<<"Number of made clusters : "<<nclusters<<"                        \n";
-  fLoader->WriteRecPoints("OVERWRITE");
-  
-  
-} // end of function
-
-//_________________________________________________________________
-void AliTPC::Hits2ExactClustersSector(Int_t isec)
-{
-  //--------------------------------------------------------
-  //calculate exact cross point of track and given pad row
-  //resulting values are expressed in "digit" coordinata
-  //--------------------------------------------------------
-
-  //-----------------------------------------------------------------
-  // Origin: Marian Ivanov  GSI Darmstadt, m.ivanov@gsi.de
-  //-----------------------------------------------------------------
-  //
-  if (fClustersArray==0){    
-    return;
-  }
-  //
-  TParticle *particle; // pointer to a given particle
-  AliTPChit *tpcHit; // pointer to a sigle TPC hit
-  //  Int_t sector,nhits;
-  Int_t ipart;
-  const Int_t kcmaxhits=30000;
-  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++) 
-    fClustersArray->CreateRow(isec,i);
-  
-  //---------------------------------------------------------------
-  //  Get the access to the tracks 
-  //---------------------------------------------------------------
-  
-  TTree *tH = TreeH();
-  if (tH == 0x0)
-   {
-     Fatal("Hits2Clusters","Can not find TreeH in folder");
-     return;
-   }
-  SetTreeAddress();
-
-  Stat_t ntracks = tH->GetEntries();
-  Int_t npart = gAlice->GetMCApp()->GetNtrack();
-  //MI change
-  TBranch * branch=0;
-  if (fHitType>1) branch = tH->GetBranch("TPC2");
-  else branch = tH->GetBranch("TPC");
+    if (sector < 36) { //inner sector [0;35]
+      if (row != 30) {
+       //the whole row is written into the output file
+       buffer->WriteRowBinary(kThreshold, digrow, 0, 0, 0, 
+                              sector, subSector, row);
+      } else {
+       //only the pads in the range [37;48] are written into the output file
+       buffer->WriteRowBinary(kThreshold, digrow, 37, 48, 1, 
+                              sector, subSector, row);
+       subSector = 1;
+       //only the pads outside the range [37;48] are written into the output file
+       buffer->WriteRowBinary(kThreshold, digrow, 37, 48, 2, 
+                              sector, subSector, row);
+      }//end else
+
+    } else { //outer sector [36;71]
+      if (row == 54) subSector = 2;
+      if ((row != 27) && (row != 76)) {
+       buffer->WriteRowBinary(kThreshold, digrow, 0, 0, 0,
+                              sector, subSector, row);
+      } else if (row == 27) {
+         //only the pads outside the range [43;46] are written into the output file
+         buffer->WriteRowBinary(kThreshold, digrow, 43, 46, 2,
+                                sector, subSector, row);
+         subSector = 1;
+         //only the pads in the range [43;46] are written into the output file
+         buffer->WriteRowBinary(kThreshold, digrow, 43, 46, 1,
+                                sector, subSector, row);
+      } else if (row == 76) {
+         //only the pads outside the range [33;88] are written into the output file
+         buffer->WriteRowBinary(kThreshold, digrow, 33, 88, 2,
+                                sector, subSector, row);
+         subSector = 3;
+         //only the pads in the range [33;88] are written into the output file
+         buffer->WriteRowBinary(kThreshold, digrow, 33, 88, 1,
+                                sector, subSector, row);
+      }
+    }//end else
+  }//end for
 
-  //------------------------------------------------------------
-  // Loop over tracks
-  //------------------------------------------------------------
+  delete buffer;
+  fLoader->UnloadDigits();
 
-  for(Int_t track=0;track<ntracks;track++){ 
-    Bool_t isInSector=kTRUE;
-    ResetHits();
-     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
-    // Loop over hits
-    //
-    Int_t currentIndex=0;
-    Int_t lastrow=-1;  //last writen row
+  AliTPCDDLRawData rawWriter;
+  rawWriter.SetVerbose(0);
 
-    //M.I. changes
+  rawWriter.RawData(fileName);
+  gSystem->Unlink(fileName);
 
-    tpcHit = (AliTPChit*)FirstHit(-1);
-    while(tpcHit){
-      
-      Int_t sector=tpcHit->fSector; // sector number
-      if(sector != isec){
-       tpcHit = (AliTPChit*) NextHit();
-       continue; 
-      }
-
-      ipart=tpcHit->Track();
-      if (ipart<npart) particle=gAlice->GetMCApp()->Particle(ipart);
-      
-      //find row number
-
-      Float_t  x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
-      Int_t    index[3]={1,isec,0};
-      Int_t    currentrow = fTPCParam->GetPadRow(x,index) ;    
-      if (currentrow<0) {tpcHit = (AliTPChit*)NextHit(); continue;}
-      if (lastrow<0) lastrow=currentrow;
-      if (currentrow==lastrow){
-       if ( currentIndex>=maxhits){
-         maxhits+=kcmaxhits;
-         xxx.ResizeTo(4*maxhits);
-       }     
-       xxx(currentIndex*4)=x[0];
-       xxx(currentIndex*4+1)=x[1];
-       xxx(currentIndex*4+2)=x[2];     
-       xxx(currentIndex*4+3)=tpcHit->fQ;
-       currentIndex++; 
-      }
-      else 
-       if (currentIndex>2){
-         Float_t sumx=0;
-         Float_t sumx2=0;
-         Float_t sumx3=0;
-         Float_t sumx4=0;
-         Float_t sumy=0;
-         Float_t sumxy=0;
-         Float_t sumx2y=0;
-         Float_t sumz=0;
-         Float_t sumxz=0;
-         Float_t sumx2z=0;
-         Float_t sumq=0;
-         for (Int_t index=0;index<currentIndex;index++){
-           Float_t x,x2,x3,x4;
-           x=x2=x3=x4=xxx(index*4);
-           x2*=x;
-           x3*=x2;
-           x4*=x3;
-           sumx+=x;
-           sumx2+=x2;
-           sumx3+=x3;
-           sumx4+=x4;
-           sumy+=xxx(index*4+1);
-           sumxy+=xxx(index*4+1)*x;
-           sumx2y+=xxx(index*4+1)*x2;
-           sumz+=xxx(index*4+2);
-           sumxz+=xxx(index*4+2)*x;
-           sumx2z+=xxx(index*4+2)*x2;   
-           sumq+=xxx(index*4+3);
-         }
-         Float_t centralPad = (fTPCParam->GetNPads(isec,lastrow)-1)/2;
-         Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
-           sumx2*(sumx*sumx3-sumx2*sumx2);
-         
-         Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
-           sumx2*(sumxy*sumx3-sumx2y*sumx2);
-         Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
-           sumx2*(sumxz*sumx3-sumx2z*sumx2);
-         
-         Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
-           sumx2*(sumx*sumx2y-sumx2*sumxy);
-         Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
-           sumx2*(sumx*sumx2z-sumx2*sumxz);
-         
-         if (TMath::Abs(det)<0.00001){
-            tpcHit = (AliTPChit*)NextHit();
-           continue;
-         }
-       
-         Float_t y=detay/det+centralPad;
-         Float_t z=detaz/det;  
-         Float_t by=detby/det; //y angle
-         Float_t bz=detbz/det; //z angle
-         sumy/=Float_t(currentIndex);
-         sumz/=Float_t(currentIndex);
-
-         AliTPCClustersRow * row = (fClustersArray->GetRow(isec,lastrow));
-         if (row!=0) {
-           AliComplexCluster* cl = new((AliComplexCluster*)row->Append()) AliComplexCluster ;
-           //    AliComplexCluster cl;
-           cl->fX=z;
-           cl->fY=y;
-           cl->fQ=sumq;
-           cl->fSigmaX2=bz;
-           cl->fSigmaY2=by;
-           cl->fTracks[0]=ipart;
-         }
-         currentIndex=0;
-         lastrow=currentrow;
-       } //end of calculating cluster for given row
-       
-       
-      tpcHit = (AliTPChit*)NextHit();
-    } // end of loop over hits
-  }   // end of loop over tracks 
-  //write padrows to tree 
-  for (Int_t ii=0; ii<fTPCParam->GetNRow(isec);ii++) {
-    fClustersArray->StoreRow(isec,ii);    
-    fClustersArray->ClearRow(isec,ii);        
+  if (kCompress) {
+    Info("Digits2Raw", "compressing raw data");
+    rawWriter.RawDataCompDecompress(kTRUE);
+    gSystem->Unlink("Statistics");
   }
-  xxxx->Delete();
 }
 
 
 
+//______________________________________________________________________
+AliDigitizer* AliTPC::CreateDigitizer(AliRunDigitizer* manager) const
+{
+  return new AliTPCDigitizer(manager);
+}
 //__
 void AliTPC::SDigits2Digits2(Int_t /*eventnumber*/)  
 {
@@ -1351,14 +1050,16 @@ void AliTPC::SDigits2Digits2(Int_t /*eventnumber*/)
 }
 //__________________________________________________________________
 void AliTPC::SetDefaults(){
-
+  //
+  // setting the defaults
+  //
    
-   cerr<<"Setting default parameters...\n";
+  //   cerr<<"Setting default parameters...\n";
 
   // Set response functions
 
   //
-  AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::fgkRunLoaderName);
+  AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::GetRunLoaderName());
   rl->CdGAFile();
   AliTPCParamSR *param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60");
   if(param){
@@ -1389,9 +1090,24 @@ void AliTPC::SetDefaults(){
     cerr<<"Can't open $ALICE_ROOT/TPC/AliTPCprf2d.root !\n" ;
      exit(3);
   }
+
+  TString s;
   prfinner->Read("prf_07504_Gati_056068_d02");
+  //PH Set different names
+  s=prfinner->GetGRF()->GetName();
+  s+="in";
+  prfinner->GetGRF()->SetName(s.Data());
+
   prfouter1->Read("prf_10006_Gati_047051_d03");
+  s=prfouter1->GetGRF()->GetName();
+  s+="out1";
+  prfouter1->GetGRF()->SetName(s.Data());
+
   prfouter2->Read("prf_15006_Gati_047051_d03");  
+  s=prfouter2->GetGRF()->GetName();
+  s+="out2";
+  prfouter2->GetGRF()->SetName(s.Data());
+
   f->Close();
   savedir->cd();
 
@@ -1409,12 +1125,32 @@ void AliTPC::SetDefaults(){
 
 }
 //__________________________________________________________________  
+void AliTPC::Hits2Digits()  
+{
+  //
+  // creates digits from hits
+  //
+
+  fLoader->LoadHits("read");
+  fLoader->LoadDigits("recreate");
+  AliRunLoader* runLoader = fLoader->GetRunLoader(); 
+
+  for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
+    runLoader->GetEvent(iEvent);
+    SetActiveSectors();   
+    Hits2Digits(iEvent);
+  }
+
+  fLoader->UnloadHits();
+  fLoader->UnloadDigits();
+} 
+//__________________________________________________________________  
 void AliTPC::Hits2Digits(Int_t eventnumber)  
 { 
  //----------------------------------------------------
  // Loop over all sectors for a single event
  //----------------------------------------------------
-  AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::fgkRunLoaderName);
+  AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::GetRunLoaderName());
   rl->GetEvent(eventnumber);
   if (fLoader->TreeH() == 0x0)
    {
@@ -1451,7 +1187,7 @@ void AliTPC::Hits2Digits(Int_t eventnumber)
 
   fDigitsSwitch=0; // standard digits
 
-  cerr<<"Digitizing TPC -- normal digits...\n";
+  //  cerr<<"Digitizing TPC -- normal digits...\n";
 
  for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) 
   if (IsSectorActive(isec)) 
@@ -1522,7 +1258,7 @@ void AliTPC::Hits2SDigits2(Int_t eventnumber)
 
   SetDigitsArray(arr);
 
-  cerr<<"Digitizing TPC -- summable digits...\n"; 
+  //  cerr<<"Digitizing TPC -- summable digits...\n"; 
 
   fDigitsSwitch=1; // summable digits
   
@@ -1555,53 +1291,19 @@ void AliTPC::Hits2SDigits()
   //   summable digits - 16 bit "ADC", no noise, no saturation
   //-----------------------------------------------------------
 
- //----------------------------------------------------
- // Loop over all sectors for a single event
- //----------------------------------------------------
-  //MI change - for pp run
-//  Int_t eventnumber = gAlice->GetEvNumber();
-//  AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::fgkRunLoaderName);
-//  rl->GetEvent(eventnumber);
-// 12/05/2003 This method is obsolete and not used. It should be redesingned
-// M.Kowalski
-
-  if (fLoader->TreeH() == 0x0)
-   {
-     if(fLoader->LoadHits())
-      {
-        Error("Hits2Digits","Can not load hits.");
-      }
-   }
-  SetTreeAddress();
-
-  if(fDefaults == 0) SetDefaults();
-  GenerNoise(500000); //create table with noise
-
-  //setup TPCDigitsArray 
-
-  if(GetDigitsArray()) delete GetDigitsArray();
+  fLoader->LoadHits("read");
+  fLoader->LoadSDigits("recreate");
+  AliRunLoader* runLoader = fLoader->GetRunLoader(); 
 
-  if (fLoader->TreeS() == 0x0 ) 
-   {
-     fLoader->MakeTree("S");
-   }
-  
-  AliTPCDigitsArray *arr = new AliTPCDigitsArray; 
-  arr->SetClass("AliSimDigits");
-  arr->Setup(fTPCParam);
-  arr->MakeTree(fLoader->TreeS());
-  SetDigitsArray(arr);
-
-  cerr<<"Digitizing TPC -- summable digits...\n"; 
-
-  //  fDigitsSwitch=1; // summable digits  -for the moment direct
-
-  for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) if (IsSectorActive(isec)) Hits2DigitsSector(isec);
+  for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
+    runLoader->GetEvent(iEvent);
+    SetTreeAddress();
+    SetActiveSectors();
+    Hits2SDigits2(iEvent);
+  }
 
-  // write results
-  //
-  cout<<"Why method TPC::Hits2SDigits writes digits and not sdigits? skowron\n";
-  fLoader->WriteDigits("OVERWRITE");
+  fLoader->UnloadHits();
+  fLoader->UnloadSDigits();
 }
 //_____________________________________________________________________________
 
@@ -2396,23 +2098,6 @@ void AliTPC::TransportElectron(Float_t *xyz, Int_t *index)
   //add nonisochronity (not implemented yet)  
 }
   
-ClassImp(AliTPCdigit)
-//_____________________________________________________________________________
-AliTPCdigit::AliTPCdigit(Int_t *tracks, Int_t *digits):
-  AliDigit(tracks)
-{
-  //
-  // Creates a TPC digit object
-  //
-  fSector     = digits[0];
-  fPadRow     = digits[1];
-  fPad        = digits[2];
-  fTime       = digits[3];
-  fSignal     = digits[4];
-}
-
 ClassImp(AliTPChit)
  
 //_____________________________________________________________________________
@@ -2571,6 +2256,9 @@ AliHit* AliTPC::FirstHit(Int_t track)
 }
 AliHit* AliTPC::NextHit()
 {
+  //
+  // gets next hit
+  //
   if (fHitType>1) return NextHit2();
   
   return AliDetector::NextHit();
@@ -2637,6 +2325,9 @@ void AliTPC::LoadPoints(Int_t)
 
 void AliTPC::RemapTrackHitIDs(Int_t *map)
 {
+  //
+  // remapping
+  //
   if (!fTrackHits) return;
   
   if (fTrackHitsOld && fHitType&2){
@@ -2863,177 +2554,6 @@ void AliTPC::LoadPoints3(Int_t)
 
 
 
-void AliTPC::FindTrackHitsIntersection(TClonesArray * /*arr*/)
-{
-
-  //
-  //fill clones array with intersection of current point with the
-  //middle of the row
-  Int_t sector;
-  Int_t ipart;
-  
-  const Int_t kcmaxhits=30000;
-  AliTPCFastVector * xxxx = new AliTPCFastVector(kcmaxhits*4);
-  AliTPCFastVector & xxx = *xxxx;
-  Int_t maxhits = kcmaxhits;
-      
-  //
-  AliTPChit * tpcHit=0;
-  tpcHit = (AliTPChit*)FirstHit2(-1);
-  Int_t currentIndex=0;
-  Int_t lastrow=-1;  //last writen row
-
-  while (tpcHit){
-    if (tpcHit==0) continue;
-    sector=tpcHit->fSector; // sector number
-    ipart=tpcHit->Track();
-    
-    //find row number
-    
-    Float_t  x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
-    Int_t    index[3]={1,sector,0};
-    Int_t    currentrow = fTPCParam->GetPadRow(x,index) ;      
-    if (currentrow<0) continue;
-    if (lastrow<0) lastrow=currentrow;
-    if (currentrow==lastrow){
-      if ( currentIndex>=maxhits){
-       maxhits+=kcmaxhits;
-       xxx.ResizeTo(4*maxhits);
-      }     
-      xxx(currentIndex*4)=x[0];
-      xxx(currentIndex*4+1)=x[1];
-      xxx(currentIndex*4+2)=x[2];      
-      xxx(currentIndex*4+3)=tpcHit->fQ;
-      currentIndex++;  
-    }
-    else 
-      if (currentIndex>2){
-       Float_t sumx=0;
-       Float_t sumx2=0;
-       Float_t sumx3=0;
-       Float_t sumx4=0;
-       Float_t sumy=0;
-       Float_t sumxy=0;
-       Float_t sumx2y=0;
-       Float_t sumz=0;
-       Float_t sumxz=0;
-       Float_t sumx2z=0;
-       Float_t sumq=0;
-       for (Int_t index=0;index<currentIndex;index++){
-         Float_t x,x2,x3,x4;
-         x=x2=x3=x4=xxx(index*4);
-         x2*=x;
-         x3*=x2;
-         x4*=x3;
-         sumx+=x;
-         sumx2+=x2;
-         sumx3+=x3;
-         sumx4+=x4;
-         sumy+=xxx(index*4+1);
-         sumxy+=xxx(index*4+1)*x;
-         sumx2y+=xxx(index*4+1)*x2;
-         sumz+=xxx(index*4+2);
-         sumxz+=xxx(index*4+2)*x;
-         sumx2z+=xxx(index*4+2)*x2;     
-         sumq+=xxx(index*4+3);
-       }
-       Float_t centralPad = (fTPCParam->GetNPads(sector,lastrow)-1)/2;
-       Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
-         sumx2*(sumx*sumx3-sumx2*sumx2);
-       
-       Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
-         sumx2*(sumxy*sumx3-sumx2y*sumx2);
-       Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
-         sumx2*(sumxz*sumx3-sumx2z*sumx2);
-       
-       Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
-         sumx2*(sumx*sumx2y-sumx2*sumxy);
-       Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
-         sumx2*(sumx*sumx2z-sumx2*sumxz);
-       
-       Float_t y=detay/det+centralPad;
-       Float_t z=detaz/det;    
-       Float_t by=detby/det; //y angle
-       Float_t bz=detbz/det; //z angle
-       sumy/=Float_t(currentIndex);
-       sumz/=Float_t(currentIndex);
-       
-       AliComplexCluster cl;
-       cl.fX=z;
-       cl.fY=y;
-       cl.fQ=sumq;
-       cl.fSigmaX2=bz;
-       cl.fSigmaY2=by;
-       cl.fTracks[0]=ipart;
-       
-       AliTPCClustersRow * row = (fClustersArray->GetRow(sector,lastrow));
-       if (row!=0) row->InsertCluster(&cl);
-       currentIndex=0;
-       lastrow=currentrow;
-      } //end of calculating cluster for given row
-               
-  } // end of loop over hits
-  xxxx->Delete();
-
-}
-//_______________________________________________________________________________
-void AliTPC::Digits2Reco(Int_t firstevent,Int_t lastevent)
-{
-  // produces rec points from digits and writes them on the root file
-  // AliTPCclusters.root
-
-  TDirectory *cwd = gDirectory;
-
-
-  AliTPCParamSR *dig=(AliTPCParamSR *)gDirectory->Get("75x40_100x60");
-  if(dig){
-    printf("You are running 2 pad-length geom hits with 3 pad-length geom digits\n");
-    delete dig;
-    dig = new AliTPCParamSR();
-  }
-  else
-  {
-   dig=(AliTPCParamSR *)gDirectory->Get("75x40_100x60_150x60"); 
-  }
-  if(!dig){
-   printf("No TPC parameters found\n");
-   exit(3);
-  }
-   
-  SetParam(dig);
-  cout<<"AliTPC::Digits2Reco: TPC parameteres have been set"<<endl; 
-  TFile *out;
-  if(!gSystem->Getenv("CONFIG_FILE")){
-    out=TFile::Open("AliTPCclusters.root","recreate");
-  }
-  else {
-    const char *tmp1;
-    const char *tmp2;
-    char tmp3[80];
-    tmp1=gSystem->Getenv("CONFIG_FILE_PREFIX");
-    tmp2=gSystem->Getenv("CONFIG_OUTDIR");
-    sprintf(tmp3,"%s%s/AliTPCclusters.root",tmp1,tmp2);
-    out=TFile::Open(tmp3,"recreate");
-  }
-
-  TStopwatch timer;
-  cout<<"AliTPC::Digits2Reco - determination of rec points begins"<<endl;
-  timer.Start();
-  cwd->cd();
-  for(Int_t iev=firstevent;iev<lastevent+1;iev++){
-
-    printf("Processing event %d\n",iev);
-    Digits2Clusters(iev);
-  }
-  cout<<"AliTPC::Digits2Reco - determination of rec points ended"<<endl;
-  timer.Stop();
-  timer.Print();
-  out->Close();
-  cwd->cd(); 
-
-
-}
-
 AliLoader* AliTPC::MakeLoader(const char* topfoldername)
 {
 //Makes TPC loader
@@ -3087,3 +2607,4 @@ AliTPCParam* AliTPC::LoadTPCParam(TFile *file) {
 
 }
 
+