]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPC.cxx
- AddTrackReference moved to AliModule
[u/mrichter/AliRoot.git] / TPC / AliTPC.cxx
index a840e735f7b1255352a86c175588d033fddd9fea..24a1d5a490acf4a11d9849a15225556e7fa29ae1 100644 (file)
 
 /*
 $Log$
+Revision 1.67  2003/03/03 16:36:16  kowal2
+Added LoadTPCParam method
+
+Revision 1.66  2003/02/11 16:54:07  hristov
+Updated AliTrackReference class (S.Radomski)
+
+Revision 1.65  2002/11/21 22:43:32  alibrary
+Removing AliMC and AliMCProcess
+
+Revision 1.64  2002/10/23 07:17:33  alibrary
+Introducing Riostream.h
+
+Revision 1.63  2002/10/14 14:57:42  hristov
+Merging the VirtualMC branch to the main development branch (HEAD)
+
+Revision 1.54.4.3  2002/10/11 08:34:48  hristov
+Updating VirtualMC to v3-09-02
+
+Revision 1.62  2002/09/23 09:22:56  hristov
+The address of the TrackReferences is set (M.Ivanov)
+
+Revision 1.61  2002/09/10 07:06:42  kowal2
+Corrected for the memory leak. Thanks to J. Chudoba and M. Ivanov
+
+Revision 1.60  2002/06/12 14:56:56  kowal2
+Added track length to the reference hits
+
+Revision 1.59  2002/06/05 15:37:31  kowal2
+Added cross-talk from the wires beyond the first and the last rows
+
+Revision 1.58  2002/05/27 14:33:14  hristov
+The new class AliTrackReference used (M.Ivanov)
+
+Revision 1.57  2002/05/07 17:23:11  kowal2
+Linear gain inefficiency instead of the step one at the wire edges.
+
+Revision 1.56  2002/04/04 16:26:33  kowal2
+Digits (Sdigits) go to separate files now.
+
+Revision 1.55  2002/03/29 06:57:45  kowal2
+Restored backward compatibility to use the hits from Dec. 2000 production.
+
+Revision 1.54  2002/03/18 17:59:13  kowal2
+Chnges in the pad geometry - 3 pad lengths introduced.
+
+Revision 1.53  2002/02/25 11:02:56  kowal2
+Changes towards speeding up the code. Thanks to Marian Ivanov.
+
 Revision 1.52  2002/02/18 09:26:09  kowal2
 Removed compiler warning
 
@@ -188,11 +236,11 @@ Introduction of the Copyright and cvs Log
 #include <TROOT.h>
 #include <TSystem.h>     
 #include "AliRun.h"
-#include <iostream.h>
+#include <Riostream.h>
 #include <stdlib.h>
-#include <fstream.h>
-#include "AliMC.h"
+#include <Riostream.h>
 #include "AliMagF.h"
+#include "AliTrackReference.h"
 
 
 #include "AliTPCParamSR.h"
@@ -242,7 +290,9 @@ AliTPCFastMatrix::AliTPCFastMatrix(Int_t row_lwb, Int_t row_upb, Int_t col_lwb,
 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){
@@ -263,7 +313,7 @@ AliTPC::AliTPC()
   fDefaults = 0;
   fTrackHits = 0; 
   fTrackHitsOld = 0;   
-  fHitType = 4; //default CONTAINERS - based on ROOT structure 
+  fHitType = 2; //default CONTAINERS - based on ROOT structure 
   fTPCParam = 0;    
   fNoiseTable = 0;
   fActiveSectors =0;
@@ -298,7 +348,7 @@ AliTPC::AliTPC(const char *name, const char *title)
 
   fNoiseTable =0;
 
-  fHitType = 4;
+  fHitType = 2;
   fActiveSectors = 0;
   //
   // Initialise counters
@@ -917,7 +967,7 @@ void    AliTPC::SetActiveSectors(Int_t * sectors, Int_t n)
     
 }
 
-void    AliTPC::SetActiveSectors()
+void    AliTPC::SetActiveSectors(Int_t flag)
 {
   //
   // activate sectors which were hitted by tracks 
@@ -925,6 +975,10 @@ void    AliTPC::SetActiveSectors()
   if (fHitType==0) return;  // if Clones hit - not short volume ID information
   if (fActiveSectors) delete [] fActiveSectors;
   fActiveSectors = new Bool_t[fTPCParam->GetNSector()];
+  if (flag) {
+    for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kTRUE;
+    return;
+  }
   for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kFALSE;
   TBranch * branch=0;
   if (fHitType>1) branch = gAlice->TreeH()->GetBranch("TPC2");
@@ -1298,6 +1352,11 @@ void AliTPC::Hits2ExactClustersSector(Int_t isec)
          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
@@ -1346,7 +1405,16 @@ void AliTPC::SDigits2Digits2(Int_t eventnumber)
   sprintf(dname,"TreeD_%s_%d",fTPCParam->GetTitle(),eventnumber);
 
   //conect tree with sSDigits
-  TTree *t = (TTree *)gDirectory->Get(sname); 
+  TTree *t;
+  if (gAlice->GetTreeDFile()) {
+    t = (TTree *)gAlice->GetTreeDFile()->Get(sname); 
+  } else {
+    t = (TTree *)gDirectory->Get(sname); 
+  }
+  if (!t) {
+    cerr<<"TPC tree with sdigits not found"<<endl;
+    return;
+  }
   AliSimDigits digarr, *dummy=&digarr;
   t->GetBranch("Segment")->SetAddress(&dummy);
   Stat_t nentries = t->GetEntries();
@@ -1364,7 +1432,12 @@ void AliTPC::SDigits2Digits2(Int_t eventnumber)
   AliTPCDigitsArray *arr = new AliTPCDigitsArray; 
   arr->SetClass("AliSimDigits");
   arr->Setup(fTPCParam);
-  arr->MakeTree(fDigitsFile);
+// Note that methods arr->MakeTree have different signatures
+  if (gAlice->GetTreeDFile()) {
+    arr->MakeTree(gAlice->GetTreeDFile());
+  } else {
+    arr->MakeTree(fDigitsFile);
+  }
   
   AliTPCParam * par =fTPCParam;
 
@@ -1430,7 +1503,7 @@ void AliTPC::SDigits2Digits2(Int_t eventnumber)
 
   
   arr->GetTree()->SetName(dname);  
-  arr->GetTree()->Write();  
+  arr->GetTree()->AutoSave();  
   delete arr;
 }
 //_________________________________________
@@ -1571,8 +1644,23 @@ void AliTPC::SetDefaults(){
   // Set response functions
 
   AliTPCParamSR *param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60");
+  if(param){
+    printf("You are using 2 pad-length geom hits with 3 pad-lenght geom digits...\n");
+    delete param;
+    param = new AliTPCParamSR();
+  }
+  else {
+    param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60_150x60");
+  }
+  if(!param){
+    printf("No TPC parameters found\n");
+    exit(4);
+  }
+
+
   AliTPCPRF2D    * prfinner   = new AliTPCPRF2D;
-  AliTPCPRF2D    * prfouter   = new AliTPCPRF2D;
+  AliTPCPRF2D    * prfouter1   = new AliTPCPRF2D;
+  AliTPCPRF2D    * prfouter2   = new AliTPCPRF2D;  
   AliTPCRF1D     * rf    = new AliTPCRF1D(kTRUE);
   rf->SetGauss(param->GetZSigma(),param->GetZWidth(),1.);
   rf->SetOffset(3*param->GetZSigma());
@@ -1585,12 +1673,14 @@ void AliTPC::SetDefaults(){
      exit(3);
   }
   prfinner->Read("prf_07504_Gati_056068_d02");
-  prfouter->Read("prf_10006_Gati_047051_d03");
+  prfouter1->Read("prf_10006_Gati_047051_d03");
+  prfouter2->Read("prf_15006_Gati_047051_d03");  
   f->Close();
   savedir->cd();
 
   param->SetInnerPRF(prfinner);
-  param->SetOuterPRF(prfouter); 
+  param->SetOuter1PRF(prfouter1); 
+  param->SetOuter2PRF(prfouter2);
   param->SetTimeRF(rf);
 
   // set fTPCParam
@@ -1619,15 +1709,20 @@ void AliTPC::Hits2Digits(Int_t eventnumber)
   AliTPCDigitsArray *arr = new AliTPCDigitsArray; 
   arr->SetClass("AliSimDigits");
   arr->Setup(fTPCParam);
-  arr->MakeTree(fDigitsFile);
+// Note that methods arr->MakeTree have different signatures
+  if (gAlice->GetTreeDFile()) {
+    arr->MakeTree(gAlice->GetTreeDFile());
+  } else {
+    arr->MakeTree(fDigitsFile);
+  }
   SetDigitsArray(arr);
 
   fDigitsSwitch=0; // standard digits
 
   cerr<<"Digitizing TPC -- normal digits...\n";
 
- for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) if (IsSectorActive(isec)) Hits2DigitsSector(isec);
-
+  for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) if (IsSectorActive(isec)) Hits2DigitsSector(isec); 
+   
   // write results
 
   char treeName[100];
@@ -1635,7 +1730,7 @@ void AliTPC::Hits2Digits(Int_t eventnumber)
   sprintf(treeName,"TreeD_%s_%d",fTPCParam->GetTitle(),eventnumber);
   
   GetDigitsArray()->GetTree()->SetName(treeName);  
-  GetDigitsArray()->GetTree()->Write();  
+  GetDigitsArray()->GetTree()->AutoSave();  
 
 
 }
@@ -1664,7 +1759,12 @@ void AliTPC::Hits2SDigits2(Int_t eventnumber)
   AliTPCDigitsArray *arr = new AliTPCDigitsArray; 
   arr->SetClass("AliSimDigits");
   arr->Setup(fTPCParam);
-  arr->MakeTree(fDigitsFile);
+// Note that methods arr->MakeTree have different signatures
+  if (gAlice->GetTreeSFile()) {
+    arr->MakeTree(gAlice->GetTreeSFile());
+  } else {
+    arr->MakeTree(fDigitsFile);
+  }
   SetDigitsArray(arr);
 
   cerr<<"Digitizing TPC -- summable digits...\n"; 
@@ -1685,7 +1785,7 @@ void AliTPC::Hits2SDigits2(Int_t eventnumber)
   sprintf(treeName,"TreeS_%s_%d",fTPCParam->GetTitle(),eventnumber);
   
   GetDigitsArray()->GetTree()->SetName(treeName); 
-  GetDigitsArray()->GetTree()->Write(); 
+  GetDigitsArray()->GetTree()->AutoSave(); 
 
 }
 
@@ -1715,7 +1815,12 @@ void AliTPC::Hits2SDigits()
   AliTPCDigitsArray *arr = new AliTPCDigitsArray; 
   arr->SetClass("AliSimDigits");
   arr->Setup(fTPCParam);
-  arr->MakeTree(fDigitsFile);
+// Note that methods arr->MakeTree have different signatures
+  if (gAlice->GetTreeSFile()) {
+    arr->MakeTree(gAlice->GetTreeSFile());
+  } else {
+    arr->MakeTree(fDigitsFile);
+  }
   SetDigitsArray(arr);
 
   cerr<<"Digitizing TPC -- summable digits...\n"; 
@@ -1731,7 +1836,7 @@ void AliTPC::Hits2SDigits()
   sprintf(treeName,"TreeD_%s_%d",fTPCParam->GetTitle(),eventnumber);
   
   GetDigitsArray()->GetTree()->SetName(treeName); 
-  GetDigitsArray()->GetTree()->Write(); 
+  GetDigitsArray()->GetTree()->AutoSave(); 
 
 }
 
@@ -1771,7 +1876,7 @@ void AliTPC::Hits2DigitsSector(Int_t isec)
 
       Int_t nrows =fTPCParam->GetNRow(isec);
 
-      row= new TObjArray* [nrows];
+      row= new TObjArray* [nrows+2]; // 2 extra rows for cross talk
     
       MakeSector(isec,nrows,tH,ntracks,row);
 
@@ -1802,7 +1907,7 @@ void AliTPC::Hits2DigitsSector(Int_t isec)
    
        } // end of the sector digitization
 
-      for(i=0;i<nrows;i++){
+      for(i=0;i<nrows+2;i++){
         row[i]->Delete();  
         delete row[i];   
       }
@@ -1829,7 +1934,7 @@ void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows)
  
 
   Float_t zerosup = fTPCParam->GetZeroSup();
-  Int_t nrows =fTPCParam->GetNRow(isec);
+  //  Int_t nrows =fTPCParam->GetNRow(isec);
   fCurrentIndex[1]= isec;
   
 
@@ -1857,14 +1962,16 @@ void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows)
   //
   //calculate signal 
   //
-  Int_t row1 = TMath::Max(irow-fTPCParam->GetNCrossRows(),0);
-  Int_t row2 = TMath::Min(irow+fTPCParam->GetNCrossRows(),nrows-1);
+  //Int_t row1 = TMath::Max(irow-fTPCParam->GetNCrossRows(),0);
+  //Int_t row2 = TMath::Min(irow+fTPCParam->GetNCrossRows(),nrows-1);
+  Int_t row1=irow;
+  Int_t row2=irow+2; 
   for (Int_t row= row1;row<=row2;row++){
     Int_t nTracks= rows[row]->GetEntries();
     for (i1=0;i1<nTracks;i1++){
       fCurrentIndex[2]= row;
-      fCurrentIndex[3]=irow;
-      if (row==irow){
+      fCurrentIndex[3]=irow+1;
+      if (row==irow+1){
        m2->Zero();  // clear single track signal matrix
        Float_t trackLabel = GetSignal(rows[row],i1,m2,m1,indexRange); 
        GetList(trackLabel,nofPads,m2,indexRange,pList);
@@ -1964,7 +2071,7 @@ Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr,
   AliTPCFastVector &v = *tv;
   
   Float_t label = v(0);
-  Int_t centralPad = (fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3])-1)/2;
+  Int_t centralPad = (fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3]-1)-1)/2;
 
   Int_t nElectrons = (tv->GetNrows()-1)/4;
   indexRange[0]=9999; // min pad
@@ -2135,10 +2242,10 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
   // of electrons, one AliTPCFastVectors per each track.
   //---------------------------------------------- 
     
-  Int_t *nofElectrons = new Int_t [nrows]; // electron counter for each row
-  AliTPCFastVector **tracks = new AliTPCFastVector* [nrows]; //pointers to the track vectors
+  Int_t *nofElectrons = new Int_t [nrows+2]; // electron counter for each row
+  AliTPCFastVector **tracks = new AliTPCFastVector* [nrows+2]; //pointers to the track vectors
 
-  for(i=0; i<nrows; i++){
+  for(i=0; i<nrows+2; i++){
     row[i] = new TObjArray;
     nofElectrons[i]=0;
     tracks[i]=0;
@@ -2185,7 +2292,7 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
                           
            // store already filled fTrack
               
-          for(i=0;i<nrows;i++){
+          for(i=0;i<nrows+2;i++){
              if(previousTrack != -1){
               if(nofElectrons[i]>0){
                  AliTPCFastVector &v = *tracks[i];
@@ -2238,13 +2345,29 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
          xyz[3]= (Float_t) (-gasgain*TMath::Log(rn)); 
          index[0]=1;
          
-         TransportElectron(xyz,index); //MI change -august       
+         TransportElectron(xyz,index);    
          Int_t rowNumber;
-         fTPCParam->GetPadRow(xyz,index); //MI change august
-         rowNumber = index[2];
+         fTPCParam->GetPadRow(xyz,index); 
+         // row 0 - cross talk from the innermost row
+         // row fNRow+1 cross talk from the outermost row
+         rowNumber = index[2]+1; 
          //transform position to local digit coordinates
          //relative to nearest pad row 
-         if ((rowNumber<0)||rowNumber>=fTPCParam->GetNRow(isec)) continue;       
+         if ((rowNumber<0)||rowNumber>fTPCParam->GetNRow(isec)+1) continue;
+          Float_t x1,y1;
+         if (isec <fTPCParam->GetNInnerSector()) {
+           x1 = xyz[1]*fTPCParam->GetInnerPadPitchWidth();
+           y1 = fTPCParam->GetYInner(rowNumber);
+         }
+         else{
+           x1=xyz[1]*fTPCParam->GetOuterPadPitchWidth();
+           y1 = fTPCParam->GetYOuter(rowNumber);
+         }
+         // gain inefficiency at the wires edges - linear
+         x1=TMath::Abs(x1);
+         y1-=1.;
+          if(x1>y1) xyz[3]*=TMath::Max(1.e-6,(y1-x1+1.));      
+       
          nofElectrons[rowNumber]++;      
          //----------------------------------
          // Expand vector if necessary
@@ -2273,7 +2396,7 @@ void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
     //   store remaining track (the last one) if not empty
     //
 
-     for(i=0;i<nrows;i++){
+     for(i=0;i<nrows+2;i++){
        if(nofElectrons[i]>0){
           AliTPCFastVector &v = *tracks[i];
          v(0) = previousTrack;
@@ -2490,13 +2613,10 @@ void AliTPC::TransportElectron(Float_t *xyz, Int_t *index)
   // ExB
   
   if (fTPCParam->GetMWPCReadout()==kTRUE){
-    Float_t x1=xyz[0];
-    fTPCParam->Transform2to2NearestWire(xyz,index);
-    Float_t dx=xyz[0]-x1;
+    Float_t dx = fTPCParam->Transform2to2NearestWire(xyz,index);
     xyz[1]+=dx*(fTPCParam->GetOmegaTau());
   }
-  //add nonisochronity (not implemented yet)
-  
+  //add nonisochronity (not implemented yet)  
 }
   
 ClassImp(AliTPCdigit)
@@ -2631,6 +2751,13 @@ void AliTPC::SetTreeAddress2()
     branch = treeH->GetBranch(branchname);
     if (branch) branch->SetAddress(&fTrackHitsOld);
   }
+  //set address to TREETR
+  TTree *treeTR = gAlice->TreeTR();
+  if (treeTR && fTrackReferences) {
+    branch = treeTR->GetBranch(GetName());
+    if (branch) branch->SetAddress(&fTrackReferences);
+  }
+
 
 }
 
@@ -3102,7 +3229,21 @@ void AliTPC::Digits2Reco(Int_t firstevent,Int_t lastevent)
   TDirectory *cwd = gDirectory;
 
 
-  AliTPCParam *dig=(AliTPCParam *)gDirectory->Get("75x40_100x60");
+  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;
@@ -3136,3 +3277,48 @@ void AliTPC::Digits2Reco(Int_t firstevent,Int_t lastevent)
 
 
 }
+
+////////////////////////////////////////////////////////////////////////
+AliTPCParam* AliTPC::LoadTPCParam(TFile *file) {
+//
+// load TPC paarmeters from a given file or create new if the object
+// is not found there
+//
+  char paramName[50];
+  sprintf(paramName,"75x40_100x60_150x60");
+  AliTPCParam *paramTPC=(AliTPCParam*)file->Get(paramName);
+  if (paramTPC) {
+    cout<<"TPC parameters "<<paramName<<" found."<<endl;
+  } else {
+    cerr<<"TPC parameters not found. Create new (they may be incorrect)."
+       <<endl;    
+    paramTPC = new AliTPCParamSR;
+  }
+  return paramTPC;
+
+// the older version of parameters can be accessed with this code.
+// In some cases, we have old parameters saved in the file but 
+// digits were created with new parameters, it can be distinguish 
+// by the name of TPC TreeD. The code here is just for the case 
+// we would need to compare with old data, uncomment it if needed.
+//
+//  char paramName[50];
+//  sprintf(paramName,"75x40_100x60");
+//  AliTPCParam *paramTPC=(AliTPCParam*)in->Get(paramName);
+//  if (paramTPC) {
+//    cout<<"TPC parameters "<<paramName<<" found."<<endl;
+//  } else {
+//    sprintf(paramName,"75x40_100x60_150x60");
+//    paramTPC=(AliTPCParam*)in->Get(paramName);
+//    if (paramTPC) {
+//     cout<<"TPC parameters "<<paramName<<" found."<<endl;
+//    } else {
+//     cerr<<"TPC parameters not found. Create new (they may be incorrect)."
+//         <<endl;    
+//     paramTPC = new AliTPCParamSR;
+//    }
+//  }
+//  return paramTPC;
+
+}
+