TPC/AliTPC.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /*
  17 $Log$
  18 Revision 1.65  2002/11/21 22:43:32  alibrary
  19 Removing AliMC and AliMCProcess
  20
  21 Revision 1.64  2002/10/23 07:17:33  alibrary
  22 Introducing Riostream.h
  23
  24 Revision 1.63  2002/10/14 14:57:42  hristov
  25 Merging the VirtualMC branch to the main development branch (HEAD)
  26
  27 Revision 1.54.4.3  2002/10/11 08:34:48  hristov
  28 Updating VirtualMC to v3-09-02
  29
  30 Revision 1.62  2002/09/23 09:22:56  hristov
  31 The address of the TrackReferences is set (M.Ivanov)
  32
  33 Revision 1.61  2002/09/10 07:06:42  kowal2
  34 Corrected for the memory leak. Thanks to J. Chudoba and M. Ivanov
  35
  36 Revision 1.60  2002/06/12 14:56:56  kowal2
  37 Added track length to the reference hits
  38
  39 Revision 1.59  2002/06/05 15:37:31  kowal2
  40 Added cross-talk from the wires beyond the first and the last rows
  41
  42 Revision 1.58  2002/05/27 14:33:14  hristov
  43 The new class AliTrackReference used (M.Ivanov)
  44
  45 Revision 1.57  2002/05/07 17:23:11  kowal2
  46 Linear gain inefficiency instead of the step one at the wire edges.
  47
  48 Revision 1.56  2002/04/04 16:26:33  kowal2
  49 Digits (Sdigits) go to separate files now.
  50
  51 Revision 1.55  2002/03/29 06:57:45  kowal2
  52 Restored backward compatibility to use the hits from Dec. 2000 production.
  53
  54 Revision 1.54  2002/03/18 17:59:13  kowal2
  55 Chnges in the pad geometry - 3 pad lengths introduced.
  56
  57 Revision 1.53  2002/02/25 11:02:56  kowal2
  58 Changes towards speeding up the code. Thanks to Marian Ivanov.
  59
  60 Revision 1.52  2002/02/18 09:26:09  kowal2
  61 Removed compiler warning
  62
  63 Revision 1.51  2002/01/21 17:13:21  kowal2
  64 New track hits using root containers. Setting active sectors added.
  65
  66 Revision 1.50  2001/12/06 14:16:19  kowal2
  67 meaningfull printouts
  68
  69 Revision 1.49  2001/11/30 11:55:37  hristov
  70 Noise table created in Hits2SDigits (M.Ivanov)
  71
  72 Revision 1.48  2001/11/24 16:10:21  kowal2
  73 Faster algorithms.
  74
  75 Revision 1.47  2001/11/19 10:25:34  kowal2
  76 Nearest integer instead of integer when converting to ADC counts
  77
  78 Revision 1.46  2001/11/07 06:47:12  kowal2
  79 Removed printouts
  80
  81 Revision 1.45  2001/11/03 13:33:48  kowal2
  82 Updated algorithms in Hits2SDigits, SDigits2Digits,
  83 Hits2ExactClusters.
  84 Added method Merge
  85
  86 Revision 1.44  2001/08/30 09:28:48  hristov
  87 TTree names are explicitly set via SetName(name) and then Write() is called
  88
  89 Revision 1.43  2001/07/28 12:02:54  hristov
  90 Branch split level set to 99
  91
  92 Revision 1.42  2001/07/28 11:38:52  hristov
  93 Loop variable declared once
  94
  95 Revision 1.41  2001/07/28 10:53:50  hristov
  96 Digitisation done according to the general scheme (M.Ivanov)
  97
  98 Revision 1.40  2001/07/27 13:03:14  hristov
  99 Default Branch split level set to 99
 100
 101 Revision 1.39  2001/07/26 09:09:34  kowal2
 102 Hits2Reco method added
 103
 104 Revision 1.38  2001/07/20 14:32:43  kowal2
 105 Processing of many events possible now
 106
 107 Revision 1.37  2001/06/12 07:17:18  kowal2
 108 Hits2SDigits method implemented (summable digits)
 109
 110 Revision 1.36  2001/05/16 14:57:25  alibrary
 111 New files for folders and Stack
 112
 113 Revision 1.35  2001/05/08 16:02:22  kowal2
 114 Updated material specifications
 115
 116 Revision 1.34  2001/05/08 15:00:15  hristov
 117 Corrections for tracking in arbitrary magnenetic field. Changes towards a concept of global Alice track. Back propagation of reconstructed tracks (Yu.Belikov)
 118
 119 Revision 1.33  2001/04/03 12:40:43  kowal2
 120 Removed printouts
 121
 122 Revision 1.32  2001/03/12 17:47:36  hristov
 123 Changes needed on Sun with CC 5.0
 124
 125 Revision 1.31  2001/03/12 08:21:50  kowal2
 126 Corrected C++ bug in the material definitions
 127
 128 Revision 1.30  2001/03/01 17:34:47  kowal2
 129 Correction due to the accuracy problem
 130
 131 Revision 1.29  2001/02/28 16:34:40  kowal2
 132 Protection against nonphysical values of the avalanche size,
 133 10**6 is the maximum
 134
 135 Revision 1.28  2001/01/26 19:57:19  hristov
 136 Major upgrade of AliRoot code
 137
 138 Revision 1.27  2001/01/13 17:29:33  kowal2
 139 Sun compiler correction
 140
 141 Revision 1.26  2001/01/10 07:59:43  kowal2
 142 Corrections to load points from the noncompressed hits.
 143
 144 Revision 1.25  2000/11/02 07:25:31  kowal2
 145 Changes due to the new hit structure.
 146 Memory leak removed.
 147
 148 Revision 1.24  2000/10/05 16:06:09  kowal2
 149 Forward declarations. Changes due to a new class AliComplexCluster.
 150
 151 Revision 1.23  2000/10/02 21:28:18  fca
 152 Removal of useless dependecies via forward declarations
 153
 154 Revision 1.22  2000/07/10 20:57:39  hristov
 155 Update of TPC code and macros by M.Kowalski
 156
 157 Revision 1.19.2.4  2000/06/26 07:39:42  kowal2
 158 Changes to obey the coding rules
 159
 160 Revision 1.19.2.3  2000/06/25 08:38:41  kowal2
 161 Splitted from AliTPCtracking
 162
 163 Revision 1.19.2.2  2000/06/16 12:59:28  kowal2
 164 Changed parameter settings
 165
 166 Revision 1.19.2.1  2000/06/09 07:15:07  kowal2
 167
 168 Defaults loaded automatically (hard-wired)
 169 Optional parameters can be set via macro called in the constructor
 170
 171 Revision 1.19  2000/04/18 19:00:59  fca
 172 Small bug fixes to TPC files
 173
 174 Revision 1.18  2000/04/17 09:37:33  kowal2
 175 removed obsolete AliTPCDigitsDisplay.C
 176
 177 Revision 1.17.2.2  2000/04/10 08:15:12  kowal2
 178
 179 New, experimental data structure from M. Ivanov
 180 New tracking algorithm
 181 Different pad geometry for different sectors
 182 Digitization rewritten
 183
 184 Revision 1.17.2.1  2000/04/10 07:56:53  kowal2
 185 Not used anymore - removed
 186
 187 Revision 1.17  2000/01/19 17:17:30  fca
 188 Introducing a list of lists of hits -- more hits allowed for detector now
 189
 190 Revision 1.16  1999/11/05 09:29:23  fca
 191 Accept only signals > 0
 192
 193 Revision 1.15  1999/10/08 06:26:53  fca
 194 Removed ClustersIndex - not used anymore
 195
 196 Revision 1.14  1999/09/29 09:24:33  fca
 197 Introduction of the Copyright and cvs Log
 198
 199 */
 200
 201 ///////////////////////////////////////////////////////////////////////////////
 202 //                                                                           //
 203 //  Time Projection Chamber                                                  //
 204 //  This class contains the basic functions for the Time Projection Chamber  //
 205 //  detector. Functions specific to one particular geometry are              //
 206 //  contained in the derived classes                                         //
 207 //                                                                           //
 208 //Begin_Html
 209 /*
 210 <img src="picts/AliTPCClass.gif">
 211 */
 212 //End_Html
 213 //                                                                           //
 214 //                                                                          //
 215 ///////////////////////////////////////////////////////////////////////////////
 216
 217 //
 218
 219 #include <TMath.h>
 220 #include <TRandom.h>
 221 #include <TVector.h>
 222 #include <TMatrix.h>
 223 #include <TGeometry.h>
 224 #include <TNode.h>
 225 #include <TTUBS.h>
 226 #include <TObjectTable.h>
 227 #include "TParticle.h"
 228 #include "AliTPC.h"
 229 #include <TFile.h>
 230 #include <TROOT.h>
 231 #include <TSystem.h>
 232 #include "AliRun.h"
 233 #include <Riostream.h>
 234 #include <stdlib.h>
 235 #include <Riostream.h>
 236 #include "AliMagF.h"
 237 #include "AliTrackReference.h"
 238
 239
 240 #include "AliTPCParamSR.h"
 241 #include "AliTPCPRF2D.h"
 242 #include "AliTPCRF1D.h"
 243 #include "AliDigits.h"
 244 #include "AliSimDigits.h"
 245 #include "AliTPCTrackHits.h"
 246 #include "AliTPCTrackHitsV2.h"
 247 #include "AliPoints.h"
 248 #include "AliArrayBranch.h"
 249
 250
 251 #include "AliTPCDigitsArray.h"
 252 #include "AliComplexCluster.h"
 253 #include "AliClusters.h"
 254 #include "AliTPCClustersRow.h"
 255 #include "AliTPCClustersArray.h"
 256
 257 #include "AliTPCcluster.h"
 258 #include "AliTPCclusterer.h"
 259 #include "AliTPCtracker.h"
 260
 261 #include <TInterpreter.h>
 262 #include <TTree.h>
 263
 264
 265
 266 ClassImp(AliTPC)
 267
 268 //_____________________________________________________________________________
 269 // helper class for fast matrix and vector manipulation - no range checking
 270 // origin - Marian Ivanov
 271
 272 class AliTPCFastMatrix : public TMatrix {
 273 public :
 274   AliTPCFastMatrix(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb);
 275   inline Float_t & UncheckedAt(Int_t rown, Int_t coln) const  {return  (fIndex[coln])[rown];} //fast acces
 276   inline Float_t   UncheckedAtFast(Int_t rown, Int_t coln) const  {return  (fIndex[coln])[rown];} //fast acces
 277 };
 278
 279 AliTPCFastMatrix::AliTPCFastMatrix(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb):
 280   TMatrix(row_lwb, row_upb,col_lwb,col_upb)
 281    {
 282    };
 283 //_____________________________________________________________________________
 284 class AliTPCFastVector : public TVector {
 285 public :
 286   AliTPCFastVector(Int_t size);
 287   virtual ~AliTPCFastVector(){;}
 288   inline Float_t & UncheckedAt(Int_t index) const  {return  fElements[index];} //fast acces
 289
 290 };
 291
 292 AliTPCFastVector::AliTPCFastVector(Int_t size):TVector(size){
 293 };
 294
 295 //_____________________________________________________________________________
 296 AliTPC::AliTPC()
 297 {
 298   //
 299   // Default constructor
 300   //
 301   fIshunt   = 0;
 302   fHits     = 0;
 303   fDigits   = 0;
 304   fNsectors = 0;
 305   fDigitsArray = 0;
 306   fClustersArray = 0;
 307   fDefaults = 0;
 308   fTrackHits = 0;
 309   fTrackHitsOld = 0;
 310   fHitType = 2; //default CONTAINERS - based on ROOT structure
 311   fTPCParam = 0;
 312   fNoiseTable = 0;
 313   fActiveSectors =0;
 314
 315 }
 316
 317 //_____________________________________________________________________________
 318 AliTPC::AliTPC(const char *name, const char *title)
 319       : AliDetector(name,title)
 320 {
 321   //
 322   // Standard constructor
 323   //
 324
 325   //
 326   // Initialise arrays of hits and digits
 327   fHits     = new TClonesArray("AliTPChit",  176);
 328   gAlice->AddHitList(fHits);
 329   fDigitsArray = 0;
 330   fClustersArray= 0;
 331   fDefaults = 0;
 332   //
 333   fTrackHits = new AliTPCTrackHitsV2;
 334   fTrackHits->SetHitPrecision(0.002);
 335   fTrackHits->SetStepPrecision(0.003);
 336   fTrackHits->SetMaxDistance(100);
 337
 338   fTrackHitsOld = new AliTPCTrackHits;  //MI - 13.09.2000
 339   fTrackHitsOld->SetHitPrecision(0.002);
 340   fTrackHitsOld->SetStepPrecision(0.003);
 341   fTrackHitsOld->SetMaxDistance(100);
 342
 343   fNoiseTable =0;
 344
 345   fHitType = 2;
 346   fActiveSectors = 0;
 347   //
 348   // Initialise counters
 349   fNsectors = 0;
 350
 351   //
 352   fIshunt     =  0;
 353   //
 354   // Initialise color attributes
 355   SetMarkerColor(kYellow);
 356
 357   //
 358   //  Set TPC parameters
 359   //
 360
 361
 362   if (!strcmp(title,"Default")) {
 363     fTPCParam = new AliTPCParamSR;
 364   } else {
 365     cerr<<"AliTPC warning: in Config.C you must set non-default parameters\n";
 366     fTPCParam=0;
 367   }
 368
 369 }
 370
 371 //_____________________________________________________________________________
 372 AliTPC::~AliTPC()
 373 {
 374   //
 375   // TPC destructor
 376   //
 377
 378   fIshunt   = 0;
 379   delete fHits;
 380   delete fDigits;
 381   delete fTPCParam;
 382   delete fTrackHits; //MI 15.09.2000
 383   delete fTrackHitsOld; //MI 10.12.2001
 384   if (fNoiseTable) delete [] fNoiseTable;
 385
 386 }
 387
 388 //_____________________________________________________________________________
 389 void AliTPC::AddHit(Int_t track, Int_t *vol, Float_t *hits)
 390 {
 391   //
 392   // Add a hit to the list
 393   //
 394   //  TClonesArray &lhits = *fHits;
 395   //  new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits);
 396   if (fHitType&1){
 397     TClonesArray &lhits = *fHits;
 398     new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits);
 399   }
 400   if (fHitType>1)
 401    AddHit2(track,vol,hits);
 402 }
 403
 404 void  AliTPC::AddTrackReference(Int_t label, TVirtualMC *vMC){
 405   //
 406   // add a trackrefernce to the list
 407   if (!fTrackReferences) {
 408     cerr<<"Container trackrefernce not active\n";
 409     return;
 410   }
 411   Int_t nref = fTrackReferences->GetEntriesFast();
 412   TClonesArray &lref = *fTrackReferences;
 413   new(lref[nref]) AliTrackReference(label, vMC);
 414 }
 415
 416 //_____________________________________________________________________________
 417 void AliTPC::BuildGeometry()
 418 {
 419
 420   //
 421   // Build TPC ROOT TNode geometry for the event display
 422   //
 423   TNode *nNode, *nTop;
 424   TTUBS *tubs;
 425   Int_t i;
 426   const int kColorTPC=19;
 427   char name[5], title[25];
 428   const Double_t kDegrad=TMath::Pi()/180;
 429   const Double_t kRaddeg=180./TMath::Pi();
 430
 431
 432   Float_t innerOpenAngle = fTPCParam->GetInnerAngle();
 433   Float_t outerOpenAngle = fTPCParam->GetOuterAngle();
 434
 435   Float_t innerAngleShift = fTPCParam->GetInnerAngleShift();
 436   Float_t outerAngleShift = fTPCParam->GetOuterAngleShift();
 437
 438   Int_t nLo = fTPCParam->GetNInnerSector()/2;
 439   Int_t nHi = fTPCParam->GetNOuterSector()/2;
 440
 441   const Double_t kloAng = (Double_t)TMath::Nint(innerOpenAngle*kRaddeg);
 442   const Double_t khiAng = (Double_t)TMath::Nint(outerOpenAngle*kRaddeg);
 443   const Double_t kloAngSh = (Double_t)TMath::Nint(innerAngleShift*kRaddeg);
 444   const Double_t khiAngSh = (Double_t)TMath::Nint(outerAngleShift*kRaddeg);
 445
 446
 447   const Double_t kloCorr = 1/TMath::Cos(0.5*kloAng*kDegrad);
 448   const Double_t khiCorr = 1/TMath::Cos(0.5*khiAng*kDegrad);
 449
 450   Double_t rl,ru;
 451
 452
 453   //
 454   // Get ALICE top node
 455   //
 456
 457   nTop=gAlice->GetGeometry()->GetNode("alice");
 458
 459   //  inner sectors
 460
 461   rl = fTPCParam->GetInnerRadiusLow();
 462   ru = fTPCParam->GetInnerRadiusUp();
 463
 464
 465   for(i=0;i<nLo;i++) {
 466     sprintf(name,"LS%2.2d",i);
 467     name[4]='\0';
 468     sprintf(title,"TPC low sector %3d",i);
 469     title[24]='\0';
 470
 471     tubs = new TTUBS(name,title,"void",rl*kloCorr,ru*kloCorr,250.,
 472                      kloAng*(i-0.5)+kloAngSh,kloAng*(i+0.5)+kloAngSh);
 473     tubs->SetNumberOfDivisions(1);
 474     nTop->cd();
 475     nNode = new TNode(name,title,name,0,0,0,"");
 476     nNode->SetLineColor(kColorTPC);
 477     fNodes->Add(nNode);
 478   }
 479
 480   // Outer sectors
 481
 482   rl = fTPCParam->GetOuterRadiusLow();
 483   ru = fTPCParam->GetOuterRadiusUp();
 484
 485   for(i=0;i<nHi;i++) {
 486     sprintf(name,"US%2.2d",i);
 487     name[4]='\0';
 488     sprintf(title,"TPC upper sector %d",i);
 489     title[24]='\0';
 490     tubs = new TTUBS(name,title,"void",rl*khiCorr,ru*khiCorr,250,
 491                      khiAng*(i-0.5)+khiAngSh,khiAng*(i+0.5)+khiAngSh);
 492     tubs->SetNumberOfDivisions(1);
 493     nTop->cd();
 494     nNode = new TNode(name,title,name,0,0,0,"");
 495     nNode->SetLineColor(kColorTPC);
 496     fNodes->Add(nNode);
 497   }
 498
 499 }
 500
 501 //_____________________________________________________________________________
 502 Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t )
 503 {
 504   //
 505   // Calculate distance from TPC to mouse on the display
 506   // Dummy procedure
 507   //
 508   return 9999;
 509 }
 510
 511 void AliTPC::Clusters2Tracks(TFile *of) {
 512   //-----------------------------------------------------------------
 513   // This is a track finder.
 514   //-----------------------------------------------------------------
 515   AliTPCtracker tracker(fTPCParam);
 516   tracker.Clusters2Tracks(gFile,of);
 517 }
 518
 519 //_____________________________________________________________________________
 520 void AliTPC::CreateMaterials()
 521 {
 522   //-----------------------------------------------
 523   // Create Materials for for TPC simulations
 524   //-----------------------------------------------
 525
 526   //-----------------------------------------------------------------
 527   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
 528   //-----------------------------------------------------------------
 529
 530   Int_t iSXFLD=gAlice->Field()->Integ();
 531   Float_t sXMGMX=gAlice->Field()->Max();
 532
 533   Float_t amat[5]; // atomic numbers
 534   Float_t zmat[5]; // z
 535   Float_t wmat[5]; // proportions
 536
 537   Float_t density;
 538   Float_t apure[2];
 539
 540
 541   //***************** Gases *************************
 542
 543   //-------------------------------------------------
 544   // pure gases
 545   //-------------------------------------------------
 546
 547   // Neon
 548
 549
 550   amat[0]= 20.18;
 551   zmat[0]= 10.;
 552   density = 0.0009;
 553
 554   apure[0]=amat[0];
 555
 556   AliMaterial(20,"Ne",amat[0],zmat[0],density,999.,999.);
 557
 558   // Argon
 559
 560   amat[0]= 39.948;
 561   zmat[0]= 18.;
 562   density = 0.001782;
 563
 564   apure[1]=amat[0];
 565
 566   AliMaterial(21,"Ar",amat[0],zmat[0],density,999.,999.);
 567
 568
 569   //--------------------------------------------------------------
 570   // gases - compounds
 571   //--------------------------------------------------------------
 572
 573   Float_t amol[3];
 574
 575   // CO2
 576
 577   amat[0]=12.011;
 578   amat[1]=15.9994;
 579
 580   zmat[0]=6.;
 581   zmat[1]=8.;
 582
 583   wmat[0]=1.;
 584   wmat[1]=2.;
 585
 586   density=0.001977;
 587
 588   amol[0] = amat[0]*wmat[0]+amat[1]*wmat[1];
 589
 590   AliMixture(10,"CO2",amat,zmat,density,-2,wmat);
 591
 592   // CF4
 593
 594   amat[0]=12.011;
 595   amat[1]=18.998;
 596
 597   zmat[0]=6.;
 598   zmat[1]=9.;
 599
 600   wmat[0]=1.;
 601   wmat[1]=4.;
 602
 603   density=0.003034;
 604
 605   amol[1] = amat[0]*wmat[0]+amat[1]*wmat[1];
 606
 607   AliMixture(11,"CF4",amat,zmat,density,-2,wmat);
 608
 609
 610   // CH4
 611
 612   amat[0]=12.011;
 613   amat[1]=1.;
 614
 615   zmat[0]=6.;
 616   zmat[1]=1.;
 617
 618   wmat[0]=1.;
 619   wmat[1]=4.;
 620
 621   density=0.000717;
 622
 623   amol[2] = amat[0]*wmat[0]+amat[1]*wmat[1];
 624
 625   AliMixture(12,"CH4",amat,zmat,density,-2,wmat);
 626
 627   //----------------------------------------------------------------
 628   // gases - mixtures, ID >= 20 pure gases, <= 10 ID < 20 -compounds
 629   //----------------------------------------------------------------
 630
 631   char namate[21];
 632   density = 0.;
 633   Float_t am=0;
 634   Int_t nc;
 635   Float_t rho,absl,X0,buf[1];
 636   Int_t nbuf;
 637   Float_t a,z;
 638
 639   for(nc = 0;nc<fNoComp;nc++)
 640     {
 641
 642       // retrive material constants
 643
 644       gMC->Gfmate((*fIdmate)[fMixtComp[nc]],namate,a,z,rho,X0,absl,buf,nbuf);
 645
 646       amat[nc] = a;
 647       zmat[nc] = z;
 648
 649       Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10;
 650
 651       am += fMixtProp[nc]*((fMixtComp[nc]>=20) ? apure[nnc] : amol[nnc]);
 652       density += fMixtProp[nc]*rho;  // density of the mixture
 653
 654     }
 655
 656   // mixture proportions by weight!
 657
 658   for(nc = 0;nc<fNoComp;nc++)
 659     {
 660
 661       Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10;
 662
 663       wmat[nc] = fMixtProp[nc]*((fMixtComp[nc]>=20) ?
 664                  apure[nnc] : amol[nnc])/am;
 665
 666     }
 667
 668   // Drift gases 1 - nonsensitive, 2 - sensitive
 669
 670   AliMixture(31,"Drift gas 1",amat,zmat,density,fNoComp,wmat);
 671   AliMixture(32,"Drift gas 2",amat,zmat,density,fNoComp,wmat);
 672
 673
 674   // Air
 675
 676   amat[0] = 14.61;
 677   zmat[0] = 7.3;
 678   density = 0.001205;
 679
 680   AliMaterial(24,"Air",amat[0],zmat[0],density,999.,999.);
 681
 682
 683   //----------------------------------------------------------------------
 684   //               solid materials
 685   //----------------------------------------------------------------------
 686
 687
 688   // Kevlar C14H22O2N2
 689
 690   amat[0] = 12.011;
 691   amat[1] = 1.;
 692   amat[2] = 15.999;
 693   amat[3] = 14.006;
 694
 695   zmat[0] = 6.;
 696   zmat[1] = 1.;
 697   zmat[2] = 8.;
 698   zmat[3] = 7.;
 699
 700   wmat[0] = 14.;
 701   wmat[1] = 22.;
 702   wmat[2] = 2.;
 703   wmat[3] = 2.;
 704
 705   density = 1.45;
 706
 707   AliMixture(34,"Kevlar",amat,zmat,density,-4,wmat);
 708
 709   // NOMEX
 710
 711   amat[0] = 12.011;
 712   amat[1] = 1.;
 713   amat[2] = 15.999;
 714   amat[3] = 14.006;
 715
 716   zmat[0] = 6.;
 717   zmat[1] = 1.;
 718   zmat[2] = 8.;
 719   zmat[3] = 7.;
 720
 721   wmat[0] = 14.;
 722   wmat[1] = 22.;
 723   wmat[2] = 2.;
 724   wmat[3] = 2.;
 725
 726   density = 0.03;
 727
 728
 729   AliMixture(35,"NOMEX",amat,zmat,density,-4,wmat);
 730
 731   // Makrolon C16H18O3
 732
 733   amat[0] = 12.011;
 734   amat[1] = 1.;
 735   amat[2] = 15.999;
 736
 737   zmat[0] = 6.;
 738   zmat[1] = 1.;
 739   zmat[2] = 8.;
 740
 741   wmat[0] = 16.;
 742   wmat[1] = 18.;
 743   wmat[2] = 3.;
 744
 745   density = 1.2;
 746
 747   AliMixture(36,"Makrolon",amat,zmat,density,-3,wmat);
 748
 749   // Mylar C5H4O2
 750
 751   amat[0]=12.011;
 752   amat[1]=1.;
 753   amat[2]=15.9994;
 754
 755   zmat[0]=6.;
 756   zmat[1]=1.;
 757   zmat[2]=8.;
 758
 759   wmat[0]=5.;
 760   wmat[1]=4.;
 761   wmat[2]=2.;
 762
 763   density = 1.39;
 764
 765   AliMixture(37, "Mylar",amat,zmat,density,-3,wmat);
 766
 767   // SiO2 - used later for the glass fiber
 768
 769   amat[0]=28.086;
 770   amat[1]=15.9994;
 771
 772   zmat[0]=14.;
 773   zmat[1]=8.;
 774
 775   wmat[0]=1.;
 776   wmat[1]=2.;
 777
 778
 779   AliMixture(38,"SiO2",amat,zmat,2.2,-2,wmat); //SiO2 - quartz (rho=2.2)
 780
 781   // Al
 782
 783   amat[0] = 26.98;
 784   zmat[0] = 13.;
 785
 786   density = 2.7;
 787
 788   AliMaterial(40,"Al",amat[0],zmat[0],density,999.,999.);
 789
 790   // Si
 791
 792   amat[0] = 28.086;
 793   zmat[0] = 14.;
 794
 795   density = 2.33;
 796
 797   AliMaterial(41,"Si",amat[0],zmat[0],density,999.,999.);
 798
 799   // Cu
 800
 801   amat[0] = 63.546;
 802   zmat[0] = 29.;
 803
 804   density = 8.96;
 805
 806   AliMaterial(42,"Cu",amat[0],zmat[0],density,999.,999.);
 807
 808   // Tedlar C2H3F
 809
 810   amat[0] = 12.011;
 811   amat[1] = 1.;
 812   amat[2] = 18.998;
 813
 814   zmat[0] = 6.;
 815   zmat[1] = 1.;
 816   zmat[2] = 9.;
 817
 818   wmat[0] = 2.;
 819   wmat[1] = 3.;
 820   wmat[2] = 1.;
 821
 822   density = 1.71;
 823
 824   AliMixture(43, "Tedlar",amat,zmat,density,-3,wmat);
 825
 826
 827   // Plexiglas  C5H8O2
 828
 829   amat[0]=12.011;
 830   amat[1]=1.;
 831   amat[2]=15.9994;
 832
 833   zmat[0]=6.;
 834   zmat[1]=1.;
 835   zmat[2]=8.;
 836
 837   wmat[0]=5.;
 838   wmat[1]=8.;
 839   wmat[2]=2.;
 840
 841   density=1.18;
 842
 843   AliMixture(44,"Plexiglas",amat,zmat,density,-3,wmat);
 844
 845   // Epoxy - C14 H20 O3
 846
 847
 848   amat[0]=12.011;
 849   amat[1]=1.;
 850   amat[2]=15.9994;
 851
 852   zmat[0]=6.;
 853   zmat[1]=1.;
 854   zmat[2]=8.;
 855
 856   wmat[0]=14.;
 857   wmat[1]=20.;
 858   wmat[2]=3.;
 859
 860   density=1.25;
 861
 862   AliMixture(45,"Epoxy",amat,zmat,density,-3,wmat);
 863
 864   // Carbon
 865
 866   amat[0]=12.011;
 867   zmat[0]=6.;
 868   density= 2.265;
 869
 870   AliMaterial(46,"C",amat[0],zmat[0],density,999.,999.);
 871
 872   // get epoxy
 873
 874   gMC->Gfmate((*fIdmate)[45],namate,amat[1],zmat[1],rho,X0,absl,buf,nbuf);
 875
 876   // Carbon fiber
 877
 878   wmat[0]=0.644; // by weight!
 879   wmat[1]=0.356;
 880
 881   density=0.5*(1.25+2.265);
 882
 883   AliMixture(47,"Cfiber",amat,zmat,density,2,wmat);
 884
 885   // get SiO2
 886
 887   gMC->Gfmate((*fIdmate)[38],namate,amat[0],zmat[0],rho,X0,absl,buf,nbuf);
 888
 889   wmat[0]=0.725; // by weight!
 890   wmat[1]=0.275;
 891
 892   density=1.7;
 893
 894   AliMixture(39,"G10",amat,zmat,density,2,wmat);
 895
 896
 897
 898
 899   //----------------------------------------------------------
 900   // tracking media for gases
 901   //----------------------------------------------------------
 902
 903   AliMedium(0, "Air", 24, 0, iSXFLD, sXMGMX, 10., 999., .1, .01, .1);
 904   AliMedium(1, "Drift gas 1", 31, 0, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
 905   AliMedium(2, "Drift gas 2", 32, 1, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
 906   AliMedium(3,"CO2",10,0, iSXFLD, sXMGMX, 10., 999.,.1, .001, .001);
 907
 908   //-----------------------------------------------------------
 909   // tracking media for solids
 910   //-----------------------------------------------------------
 911
 912   AliMedium(4,"Al",40,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 913   AliMedium(5,"Kevlar",34,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 914   AliMedium(6,"Nomex",35,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 915   AliMedium(7,"Makrolon",36,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 916   AliMedium(8,"Mylar",37,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 917   AliMedium(9,"Tedlar",43,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 918   AliMedium(10,"Cu",42,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 919   AliMedium(11,"Si",41,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 920   AliMedium(12,"G10",39,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 921   AliMedium(13,"Plexiglas",44,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 922   AliMedium(14,"Epoxy",45,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 923   AliMedium(15,"Cfiber",47,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 924
 925 }
 926
 927 void AliTPC::GenerNoise(Int_t tablesize)
 928 {
 929   //
 930   //Generate table with noise
 931   //
 932   if (fTPCParam==0) {
 933     // error message
 934     fNoiseDepth=0;
 935     return;
 936   }
 937   if (fNoiseTable)  delete[] fNoiseTable;
 938   fNoiseTable = new Float_t[tablesize];
 939   fNoiseDepth = tablesize;
 940   fCurrentNoise =0; //!index of the noise in  the noise table
 941
 942   Float_t norm = fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac();
 943   for (Int_t i=0;i<tablesize;i++) fNoiseTable[i]= gRandom->Gaus(0,norm);
 944 }
 945
 946 Float_t AliTPC::GetNoise()
 947 {
 948   // get noise from table
 949   //  if ((fCurrentNoise%10)==0)
 950   //  fCurrentNoise= gRandom->Rndm()*fNoiseDepth;
 951   if (fCurrentNoise>=fNoiseDepth) fCurrentNoise=0;
 952   return fNoiseTable[fCurrentNoise++];
 953   //gRandom->Gaus(0, fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac());
 954 }
 955
 956
 957 Bool_t  AliTPC::IsSectorActive(Int_t sec)
 958 {
 959   //
 960   // check if the sector is active
 961   if (!fActiveSectors) return kTRUE;
 962   else return fActiveSectors[sec];
 963 }
 964
 965 void    AliTPC::SetActiveSectors(Int_t * sectors, Int_t n)
 966 {
 967   // activate interesting sectors
 968   if (fActiveSectors) delete [] fActiveSectors;
 969   fActiveSectors = new Bool_t[fTPCParam->GetNSector()];
 970   for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kFALSE;
 971   for (Int_t i=0;i<n;i++)
 972     if ((sectors[i]>=0) && sectors[i]<fTPCParam->GetNSector())  fActiveSectors[sectors[i]]=kTRUE;
 973
 974 }
 975
 976 void    AliTPC::SetActiveSectors(Int_t flag)
 977 {
 978   //
 979   // activate sectors which were hitted by tracks
 980   //loop over tracks
 981   if (fHitType==0) return;  // if Clones hit - not short volume ID information
 982   if (fActiveSectors) delete [] fActiveSectors;
 983   fActiveSectors = new Bool_t[fTPCParam->GetNSector()];
 984   if (flag) {
 985     for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kTRUE;
 986     return;
 987   }
 988   for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kFALSE;
 989   TBranch * branch=0;
 990   if (fHitType>1) branch = gAlice->TreeH()->GetBranch("TPC2");
 991   else branch = gAlice->TreeH()->GetBranch("TPC");
 992   Stat_t ntracks = gAlice->TreeH()->GetEntries();
 993   // loop over all hits
 994   for(Int_t track=0;track<ntracks;track++){
 995     ResetHits();
 996     //
 997     if (fTrackHits && fHitType&4) {
 998       TBranch * br1 = gAlice->TreeH()->GetBranch("fVolumes");
 999       TBranch * br2 = gAlice->TreeH()->GetBranch("fNVolumes");
1000       br1->GetEvent(track);
1001       br2->GetEvent(track);
1002       Int_t *volumes = fTrackHits->GetVolumes();
1003       for (Int_t j=0;j<fTrackHits->GetNVolumes(); j++)
1004         fActiveSectors[volumes[j]]=kTRUE;
1005     }
1006
1007     //
1008     if (fTrackHitsOld && fHitType&2) {
1009       TBranch * br = gAlice->TreeH()->GetBranch("fTrackHitsInfo");
1010       br->GetEvent(track);
1011       AliObjectArray * ar = fTrackHitsOld->fTrackHitsInfo;
1012       for (UInt_t j=0;j<ar->GetSize();j++){
1013         fActiveSectors[((AliTrackHitsInfo*)ar->At(j))->fVolumeID] =kTRUE;
1014       }
1015     }
1016   }
1017
1018 }
1019
1020
1021
1022
1023 void AliTPC::Digits2Clusters(TFile *of, Int_t eventnumber)
1024 {
1025   //-----------------------------------------------------------------
1026   // This is a simple cluster finder.
1027   //-----------------------------------------------------------------
1028   AliTPCclusterer::Digits2Clusters(fTPCParam,of,eventnumber);
1029 }
1030
1031 extern Double_t SigmaY2(Double_t, Double_t, Double_t);
1032 extern Double_t SigmaZ2(Double_t, Double_t);
1033 //_____________________________________________________________________________
1034 void AliTPC::Hits2Clusters(TFile *of, Int_t eventn)
1035 {
1036   //--------------------------------------------------------
1037   // TPC simple cluster generator from hits
1038   // obtained from the TPC Fast Simulator
1039   // The point errors are taken from the parametrization
1040   //--------------------------------------------------------
1041
1042   //-----------------------------------------------------------------
1043   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1044   //-----------------------------------------------------------------
1045   // Adopted to Marian's cluster data structure by I.Belikov, CERN,
1046   // Jouri.Belikov@cern.ch
1047   //----------------------------------------------------------------
1048
1049   /////////////////////////////////////////////////////////////////////////////
1050   //
1051   //---------------------------------------------------------------------
1052   //   ALICE TPC Cluster Parameters
1053   //--------------------------------------------------------------------
1054
1055
1056
1057   // Cluster width in rphi
1058   const Float_t kACrphi=0.18322;
1059   const Float_t kBCrphi=0.59551e-3;
1060   const Float_t kCCrphi=0.60952e-1;
1061   // Cluster width in z
1062   const Float_t kACz=0.19081;
1063   const Float_t kBCz=0.55938e-3;
1064   const Float_t kCCz=0.30428;
1065
1066   TDirectory *savedir=gDirectory;
1067
1068   if (!of->IsOpen()) {
1069      cerr<<"AliTPC::Hits2Clusters(): output file not open !\n";
1070      return;
1071   }
1072
1073   //if(fDefaults == 0) SetDefaults();
1074
1075   Float_t sigmaRphi,sigmaZ,clRphi,clZ;
1076   //
1077   TParticle *particle; // pointer to a given particle
1078   AliTPChit *tpcHit; // pointer to a sigle TPC hit
1079   Int_t sector;
1080   Int_t ipart;
1081   Float_t xyz[5];
1082   Float_t pl,pt,tanth,rpad,ratio;
1083   Float_t cph,sph;
1084
1085   //---------------------------------------------------------------
1086   //  Get the access to the tracks
1087   //---------------------------------------------------------------
1088
1089   TTree *tH = gAlice->TreeH();
1090   Stat_t ntracks = tH->GetEntries();
1091
1092   //Switch to the output file
1093   of->cd();
1094
1095   char   cname[100];
1096
1097   sprintf(cname,"TreeC_TPC_%d",eventn);
1098
1099   fTPCParam->Write(fTPCParam->GetTitle());
1100   AliTPCClustersArray carray;
1101   carray.Setup(fTPCParam);
1102   carray.SetClusterType("AliTPCcluster");
1103   carray.MakeTree();
1104
1105   Int_t nclusters=0; //cluster counter
1106
1107   //------------------------------------------------------------
1108   // Loop over all sectors (72 sectors for 20 deg
1109   // segmentation for both lower and upper sectors)
1110   // Sectors 0-35 are lower sectors, 0-17 z>0, 17-35 z<0
1111   // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
1112   //
1113   // First cluster for sector 0 starts at "0"
1114   //------------------------------------------------------------
1115
1116   for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++){
1117     //MI change
1118     fTPCParam->AdjustCosSin(isec,cph,sph);
1119
1120     //------------------------------------------------------------
1121     // Loop over tracks
1122     //------------------------------------------------------------
1123
1124     for(Int_t track=0;track<ntracks;track++){
1125       ResetHits();
1126       tH->GetEvent(track);
1127       //
1128       //  Get number of the TPC hits
1129       //
1130        tpcHit = (AliTPChit*)FirstHit(-1);
1131
1132       // Loop over hits
1133       //
1134        while(tpcHit){
1135
1136          if (tpcHit->fQ == 0.) {
1137            tpcHit = (AliTPChit*) NextHit();
1138            continue; //information about track (I.Belikov)
1139          }
1140         sector=tpcHit->fSector; // sector number
1141
1142        if(sector != isec){
1143          tpcHit = (AliTPChit*) NextHit();
1144          continue;
1145        }
1146         ipart=tpcHit->Track();
1147         particle=gAlice->Particle(ipart);
1148         pl=particle->Pz();
1149         pt=particle->Pt();
1150         if(pt < 1.e-9) pt=1.e-9;
1151         tanth=pl/pt;
1152         tanth = TMath::Abs(tanth);
1153         rpad=TMath::Sqrt(tpcHit->X()*tpcHit->X() + tpcHit->Y()*tpcHit->Y());
1154         ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
1155
1156         //   space-point resolutions
1157
1158         sigmaRphi=SigmaY2(rpad,tanth,pt);
1159         sigmaZ   =SigmaZ2(rpad,tanth   );
1160
1161         //   cluster widths
1162
1163         clRphi=kACrphi-kBCrphi*rpad*tanth+kCCrphi*ratio*ratio;
1164         clZ=kACz-kBCz*rpad*tanth+kCCz*tanth*tanth;
1165
1166         // temporary protection
1167
1168         if(sigmaRphi < 0.) sigmaRphi=0.4e-3;
1169         if(sigmaZ < 0.) sigmaZ=0.4e-3;
1170         if(clRphi < 0.) clRphi=2.5e-3;
1171         if(clZ < 0.) clZ=2.5e-5;
1172
1173         //
1174
1175         //
1176         // smearing --> rotate to the 1 (13) or to the 25 (49) sector,
1177         // then the inaccuracy in a X-Y plane is only along Y (pad row)!
1178         //
1179         Float_t xprim= tpcHit->X()*cph + tpcHit->Y()*sph;
1180         Float_t yprim=-tpcHit->X()*sph + tpcHit->Y()*cph;
1181         xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigmaRphi));   // y
1182           Float_t alpha=(isec < fTPCParam->GetNInnerSector()) ?
1183           fTPCParam->GetInnerAngle() : fTPCParam->GetOuterAngle();
1184           Float_t ymax=xprim*TMath::Tan(0.5*alpha);
1185           if (TMath::Abs(xyz[0])>ymax) xyz[0]=yprim;
1186         xyz[1]=gRandom->Gaus(tpcHit->Z(),TMath::Sqrt(sigmaZ)); // z
1187           if (TMath::Abs(xyz[1])>fTPCParam->GetZLength()) xyz[1]=tpcHit->Z();
1188         xyz[2]=sigmaRphi;                                     // fSigmaY2
1189         xyz[3]=sigmaZ;                                        // fSigmaZ2
1190         xyz[4]=tpcHit->fQ;                                    // q
1191
1192         AliTPCClustersRow *clrow=carray.GetRow(sector,tpcHit->fPadRow);
1193         if (!clrow) clrow=carray.CreateRow(sector,tpcHit->fPadRow);
1194
1195         Int_t tracks[3]={tpcHit->Track(), -1, -1};
1196         AliTPCcluster cluster(tracks,xyz);
1197
1198         clrow->InsertCluster(&cluster); nclusters++;
1199
1200         tpcHit = (AliTPChit*)NextHit();
1201
1202
1203       } // end of loop over hits
1204
1205     }   // end of loop over tracks
1206
1207     Int_t nrows=fTPCParam->GetNRow(isec);
1208     for (Int_t irow=0; irow<nrows; irow++) {
1209         AliTPCClustersRow *clrow=carray.GetRow(isec,irow);
1210         if (!clrow) continue;
1211         carray.StoreRow(isec,irow);
1212         carray.ClearRow(isec,irow);
1213     }
1214
1215   } // end of loop over sectors
1216
1217   cerr<<"Number of made clusters : "<<nclusters<<"                        \n";
1218   carray.GetTree()->SetName(cname);
1219   carray.GetTree()->Write();
1220
1221   savedir->cd(); //switch back to the input file
1222
1223 } // end of function
1224
1225 //_________________________________________________________________
1226 void AliTPC::Hits2ExactClustersSector(Int_t isec)
1227 {
1228   //--------------------------------------------------------
1229   //calculate exact cross point of track and given pad row
1230   //resulting values are expressed in "digit" coordinata
1231   //--------------------------------------------------------
1232
1233   //-----------------------------------------------------------------
1234   // Origin: Marian Ivanov  GSI Darmstadt, m.ivanov@gsi.de
1235   //-----------------------------------------------------------------
1236   //
1237   if (fClustersArray==0){
1238     return;
1239   }
1240   //
1241   TParticle *particle; // pointer to a given particle
1242   AliTPChit *tpcHit; // pointer to a sigle TPC hit
1243   //  Int_t sector,nhits;
1244   Int_t ipart;
1245   const Int_t kcmaxhits=30000;
1246   AliTPCFastVector * xxxx = new AliTPCFastVector(kcmaxhits*4);
1247   AliTPCFastVector & xxx = *xxxx;
1248   Int_t maxhits = kcmaxhits;
1249   //construct array for each padrow
1250   for (Int_t i=0; i<fTPCParam->GetNRow(isec);i++)
1251     fClustersArray->CreateRow(isec,i);
1252
1253   //---------------------------------------------------------------
1254   //  Get the access to the tracks
1255   //---------------------------------------------------------------
1256
1257   TTree *tH = gAlice->TreeH();
1258   Stat_t ntracks = tH->GetEntries();
1259   Int_t npart = gAlice->GetNtrack();
1260   //MI change
1261   TBranch * branch=0;
1262   if (fHitType>1) branch = tH->GetBranch("TPC2");
1263   else branch = tH->GetBranch("TPC");
1264
1265   //------------------------------------------------------------
1266   // Loop over tracks
1267   //------------------------------------------------------------
1268
1269   for(Int_t track=0;track<ntracks;track++){
1270     Bool_t isInSector=kTRUE;
1271     ResetHits();
1272      isInSector = TrackInVolume(isec,track);
1273     if (!isInSector) continue;
1274     //MI change
1275     branch->GetEntry(track); // get next track
1276     //
1277     //  Get number of the TPC hits and a pointer
1278     //  to the particles
1279     // Loop over hits
1280     //
1281     Int_t currentIndex=0;
1282     Int_t lastrow=-1;  //last writen row
1283
1284     //M.I. changes
1285
1286     tpcHit = (AliTPChit*)FirstHit(-1);
1287     while(tpcHit){
1288
1289       Int_t sector=tpcHit->fSector; // sector number
1290       if(sector != isec){
1291         tpcHit = (AliTPChit*) NextHit();
1292         continue;
1293       }
1294
1295       ipart=tpcHit->Track();
1296       if (ipart<npart) particle=gAlice->Particle(ipart);
1297
1298       //find row number
1299
1300       Float_t  x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
1301       Int_t    index[3]={1,isec,0};
1302       Int_t    currentrow = fTPCParam->GetPadRow(x,index) ;
1303       if (currentrow<0) {tpcHit = (AliTPChit*)NextHit(); continue;}
1304       if (lastrow<0) lastrow=currentrow;
1305       if (currentrow==lastrow){
1306         if ( currentIndex>=maxhits){
1307           maxhits+=kcmaxhits;
1308           xxx.ResizeTo(4*maxhits);
1309         }
1310         xxx(currentIndex*4)=x[0];
1311         xxx(currentIndex*4+1)=x[1];
1312         xxx(currentIndex*4+2)=x[2];
1313         xxx(currentIndex*4+3)=tpcHit->fQ;
1314         currentIndex++;
1315       }
1316       else
1317         if (currentIndex>2){
1318           Float_t sumx=0;
1319           Float_t sumx2=0;
1320           Float_t sumx3=0;
1321           Float_t sumx4=0;
1322           Float_t sumy=0;
1323           Float_t sumxy=0;
1324           Float_t sumx2y=0;
1325           Float_t sumz=0;
1326           Float_t sumxz=0;
1327           Float_t sumx2z=0;
1328           Float_t sumq=0;
1329           for (Int_t index=0;index<currentIndex;index++){
1330             Float_t x,x2,x3,x4;
1331             x=x2=x3=x4=xxx(index*4);
1332             x2*=x;
1333             x3*=x2;
1334             x4*=x3;
1335             sumx+=x;
1336             sumx2+=x2;
1337             sumx3+=x3;
1338             sumx4+=x4;
1339             sumy+=xxx(index*4+1);
1340             sumxy+=xxx(index*4+1)*x;
1341             sumx2y+=xxx(index*4+1)*x2;
1342             sumz+=xxx(index*4+2);
1343             sumxz+=xxx(index*4+2)*x;
1344             sumx2z+=xxx(index*4+2)*x2;
1345             sumq+=xxx(index*4+3);
1346           }
1347           Float_t centralPad = (fTPCParam->GetNPads(isec,lastrow)-1)/2;
1348           Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
1349             sumx2*(sumx*sumx3-sumx2*sumx2);
1350
1351           Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
1352             sumx2*(sumxy*sumx3-sumx2y*sumx2);
1353           Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
1354             sumx2*(sumxz*sumx3-sumx2z*sumx2);
1355
1356           Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
1357             sumx2*(sumx*sumx2y-sumx2*sumxy);
1358           Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
1359             sumx2*(sumx*sumx2z-sumx2*sumxz);
1360
1361           if (TMath::Abs(det)<0.00001){
1362              tpcHit = (AliTPChit*)NextHit();
1363             continue;
1364           }
1365
1366           Float_t y=detay/det+centralPad;
1367           Float_t z=detaz/det;
1368           Float_t by=detby/det; //y angle
1369           Float_t bz=detbz/det; //z angle
1370           sumy/=Float_t(currentIndex);
1371           sumz/=Float_t(currentIndex);
1372
1373           AliTPCClustersRow * row = (fClustersArray->GetRow(isec,lastrow));
1374           if (row!=0) {
1375             AliComplexCluster* cl = new((AliComplexCluster*)row->Append()) AliComplexCluster ;
1376             //    AliComplexCluster cl;
1377             cl->fX=z;
1378             cl->fY=y;
1379             cl->fQ=sumq;
1380             cl->fSigmaX2=bz;
1381             cl->fSigmaY2=by;
1382             cl->fTracks[0]=ipart;
1383           }
1384           currentIndex=0;
1385           lastrow=currentrow;
1386         } //end of calculating cluster for given row
1387
1388
1389       tpcHit = (AliTPChit*)NextHit();
1390     } // end of loop over hits
1391   }   // end of loop over tracks
1392   //write padrows to tree
1393   for (Int_t ii=0; ii<fTPCParam->GetNRow(isec);ii++) {
1394     fClustersArray->StoreRow(isec,ii);
1395     fClustersArray->ClearRow(isec,ii);
1396   }
1397   xxxx->Delete();
1398
1399 }
1400
1401
1402
1403 //__
1404 void AliTPC::SDigits2Digits2(Int_t eventnumber)
1405 {
1406   //create digits from summable digits
1407   GenerNoise(500000); //create teble with noise
1408   char  sname[100];
1409   char  dname[100];
1410   sprintf(sname,"TreeS_%s_%d",fTPCParam->GetTitle(),eventnumber);
1411   sprintf(dname,"TreeD_%s_%d",fTPCParam->GetTitle(),eventnumber);
1412
1413   //conect tree with sSDigits
1414   TTree *t;
1415   if (gAlice->GetTreeDFile()) {
1416     t = (TTree *)gAlice->GetTreeDFile()->Get(sname);
1417   } else {
1418     t = (TTree *)gDirectory->Get(sname);
1419   }
1420   if (!t) {
1421     cerr<<"TPC tree with sdigits not found"<<endl;
1422     return;
1423   }
1424   AliSimDigits digarr, *dummy=&digarr;
1425   t->GetBranch("Segment")->SetAddress(&dummy);
1426   Stat_t nentries = t->GetEntries();
1427
1428   // set zero suppression
1429
1430   fTPCParam->SetZeroSup(2);
1431
1432   // get zero suppression
1433
1434   Int_t zerosup = fTPCParam->GetZeroSup();
1435
1436   //make tree with digits
1437
1438   AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1439   arr->SetClass("AliSimDigits");
1440   arr->Setup(fTPCParam);
1441 // Note that methods arr->MakeTree have different signatures
1442   if (gAlice->GetTreeDFile()) {
1443     arr->MakeTree(gAlice->GetTreeDFile());
1444   } else {
1445     arr->MakeTree(fDigitsFile);
1446   }
1447
1448   AliTPCParam * par =fTPCParam;
1449
1450   //Loop over segments of the TPC
1451
1452   for (Int_t n=0; n<nentries; n++) {
1453     t->GetEvent(n);
1454     Int_t sec, row;
1455     if (!par->AdjustSectorRow(digarr.GetID(),sec,row)) {
1456       cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
1457       continue;
1458     }
1459     if (!IsSectorActive(sec)) continue;
1460     AliSimDigits * digrow =(AliSimDigits*) arr->CreateRow(sec,row);
1461     Int_t nrows = digrow->GetNRows();
1462     Int_t ncols = digrow->GetNCols();
1463
1464     digrow->ExpandBuffer();
1465     digarr.ExpandBuffer();
1466     digrow->ExpandTrackBuffer();
1467     digarr.ExpandTrackBuffer();
1468
1469
1470     Short_t * pamp0 = digarr.GetDigits();
1471     Int_t   * ptracks0 = digarr.GetTracks();
1472     Short_t * pamp1 = digrow->GetDigits();
1473     Int_t   * ptracks1 = digrow->GetTracks();
1474     Int_t  nelems =nrows*ncols;
1475     Int_t saturation = fTPCParam->GetADCSat();
1476     //use internal structure of the AliDigits - for speed reason
1477     //if you cahnge implementation
1478     //of the Alidigits - it must be rewriten -
1479     for (Int_t i= 0; i<nelems; i++){
1480       //      Float_t q = *pamp0;
1481       //q/=16.;  //conversion faktor
1482       //Float_t noise= GetNoise();
1483       //q+=noise;
1484       //q= TMath::Nint(q);
1485       Float_t q = TMath::Nint(Float_t(*pamp0)/16.+GetNoise());
1486       if (q>zerosup){
1487         if (q>saturation) q=saturation;
1488         *pamp1=(Short_t)q;
1489         //if (ptracks0[0]==0)
1490         //  ptracks1[0]=1;
1491         //else
1492         ptracks1[0]=ptracks0[0];
1493         ptracks1[nelems]=ptracks0[nelems];
1494         ptracks1[2*nelems]=ptracks0[2*nelems];
1495       }
1496       pamp0++;
1497       pamp1++;
1498       ptracks0++;
1499       ptracks1++;
1500     }
1501
1502     arr->StoreRow(sec,row);
1503     arr->ClearRow(sec,row);
1504     // cerr<<sec<<"\t"<<row<<"\n";
1505   }
1506
1507
1508   //write results
1509
1510
1511   arr->GetTree()->SetName(dname);
1512   arr->GetTree()->AutoSave();
1513   delete arr;
1514 }
1515 //_________________________________________
1516 void AliTPC::Merge(TTree * intree, Int_t *mask, Int_t nin, Int_t outid)
1517 {
1518
1519   //intree - pointer to array of trees with s digits
1520   //mask   - mask for each
1521   //nin    - number of inputs
1522   //outid  - event  number of the output event
1523
1524   //create digits from summable digits
1525   //conect tree with sSDigits
1526
1527
1528   AliSimDigits ** digarr = new AliSimDigits*[nin];
1529   for (Int_t i1=0;i1<nin; i1++){
1530     digarr[i1]=0;
1531     intree[i1].GetBranch("Segment")->SetAddress(&digarr[i1]);
1532   }
1533   Stat_t nentries = intree[0].GetEntries();
1534
1535   //make tree with digits
1536   char  dname[100];
1537   sprintf(dname,"TreeD_%s_%d",fTPCParam->GetTitle(),outid);
1538   AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1539   arr->SetClass("AliSimDigits");
1540   arr->Setup(fTPCParam);
1541   arr->MakeTree(fDigitsFile);
1542
1543   // set zero suppression
1544
1545   fTPCParam->SetZeroSup(2);
1546
1547   // get zero suppression
1548
1549   Int_t zerosup = fTPCParam->GetZeroSup();
1550
1551
1552   AliTPCParam * par =fTPCParam;
1553
1554   //Loop over segments of the TPC
1555   for (Int_t n=0; n<nentries; n++) {
1556
1557     for (Int_t i=0;i<nin; i++){
1558       //connect proper digits
1559       intree[i].GetEvent(n);
1560       digarr[i]->ExpandBuffer();
1561       digarr[i]->ExpandTrackBuffer();
1562     }
1563     Int_t sec, row;
1564     if (!par->AdjustSectorRow(digarr[0]->GetID(),sec,row)) {
1565       cerr<<"AliTPC warning: invalid segment ID ! "<<digarr[0]->GetID()<<endl;
1566       continue;
1567     }
1568
1569     AliSimDigits * digrow =(AliSimDigits*) arr->CreateRow(sec,row);
1570     Int_t nrows = digrow->GetNRows();
1571     Int_t ncols = digrow->GetNCols();
1572
1573     digrow->ExpandBuffer();
1574     digrow->ExpandTrackBuffer();
1575
1576     for (Int_t rows=0;rows<nrows; rows++){
1577       for (Int_t col=0;col<ncols; col++){
1578         Float_t q=0;
1579         Int_t label[1000]; // i hope no more than 300 events merged
1580         Int_t labptr = 0;
1581         // looop over digits
1582         for (Int_t i=0;i<nin; i++){
1583           q  += digarr[i]->GetDigitFast(rows,col);
1584           for (Int_t tr=0;tr<3;tr++) {
1585             Int_t lab = digarr[i]->GetTrackIDFast(rows,col,tr);
1586             if ( lab > 1) {
1587               label[labptr]=lab+mask[i]-2;
1588               labptr++;
1589             }
1590           }
1591         }
1592         //add noise
1593         q = gRandom->Gaus(q,fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac()*16);
1594
1595         q/=16;  //conversion faktor
1596         q=(Short_t)q;
1597
1598         if (q> zerosup){
1599
1600           if(q > fTPCParam->GetADCSat()) q = fTPCParam->GetADCSat();
1601           digrow->SetDigitFast((Short_t)q,rows,col);
1602           for (Int_t tr=0;tr<3;tr++){
1603             if (tr<labptr)
1604               ((AliSimDigits*)digrow)->SetTrackIDFast(label[tr],rows,col,tr);
1605             else
1606               ((AliSimDigits*)digrow)->SetTrackIDFast(-1,rows,col,tr);
1607           }
1608         }
1609       }
1610      }
1611       arr->StoreRow(sec,row);
1612
1613       arr->ClearRow(sec,row);
1614
1615   }
1616
1617   delete digarr;
1618   arr->GetTree()->SetName(dname);
1619   arr->GetTree()->Write();
1620
1621   delete arr;
1622
1623 }
1624
1625 /*_________________________________________
1626 void AliTPC::SDigits2Digits(Int_t eventnumber)
1627 {
1628
1629
1630   cerr<<"Digitizing TPC...\n";
1631
1632   Hits2Digits(eventnumber);
1633
1634
1635   //write results
1636
1637   //  char treeName[100];
1638
1639   //  sprintf(treeName,"TreeD_%s_%d",fTPCParam->GetTitle(),eventnumber);
1640
1641   //  GetDigitsArray()->GetTree()->Write(treeName);
1642 }
1643 */
1644 //__________________________________________________________________
1645 void AliTPC::SetDefaults(){
1646
1647
1648    cerr<<"Setting default parameters...\n";
1649
1650   // Set response functions
1651
1652   AliTPCParamSR *param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60");
1653   if(param){
1654     printf("You are using 2 pad-length geom hits with 3 pad-lenght geom digits...\n");
1655     delete param;
1656     param = new AliTPCParamSR();
1657   }
1658   else {
1659     param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60_150x60");
1660   }
1661   if(!param){
1662     printf("No TPC parameters found\n");
1663     exit(4);
1664   }
1665
1666
1667   AliTPCPRF2D    * prfinner   = new AliTPCPRF2D;
1668   AliTPCPRF2D    * prfouter1   = new AliTPCPRF2D;
1669   AliTPCPRF2D    * prfouter2   = new AliTPCPRF2D;
1670   AliTPCRF1D     * rf    = new AliTPCRF1D(kTRUE);
1671   rf->SetGauss(param->GetZSigma(),param->GetZWidth(),1.);
1672   rf->SetOffset(3*param->GetZSigma());
1673   rf->Update();
1674
1675   TDirectory *savedir=gDirectory;
1676   TFile *f=TFile::Open("$ALICE_ROOT/TPC/AliTPCprf2d.root");
1677   if (!f->IsOpen()) {
1678     cerr<<"Can't open $ALICE_ROOT/TPC/AliTPCprf2d.root !\n" ;
1679      exit(3);
1680   }
1681   prfinner->Read("prf_07504_Gati_056068_d02");
1682   prfouter1->Read("prf_10006_Gati_047051_d03");
1683   prfouter2->Read("prf_15006_Gati_047051_d03");
1684   f->Close();
1685   savedir->cd();
1686
1687   param->SetInnerPRF(prfinner);
1688   param->SetOuter1PRF(prfouter1);
1689   param->SetOuter2PRF(prfouter2);
1690   param->SetTimeRF(rf);
1691
1692   // set fTPCParam
1693
1694   SetParam(param);
1695
1696
1697   fDefaults = 1;
1698
1699 }
1700 //__________________________________________________________________
1701 void AliTPC::Hits2Digits(Int_t eventnumber)
1702 {
1703  //----------------------------------------------------
1704  // Loop over all sectors for a single event
1705  //----------------------------------------------------
1706
1707
1708   if(fDefaults == 0) SetDefaults();  // check if the parameters are set
1709   GenerNoise(500000); //create teble with noise
1710
1711   //setup TPCDigitsArray
1712
1713   if(GetDigitsArray()) delete GetDigitsArray();
1714
1715   AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1716   arr->SetClass("AliSimDigits");
1717   arr->Setup(fTPCParam);
1718 // Note that methods arr->MakeTree have different signatures
1719   if (gAlice->GetTreeDFile()) {
1720     arr->MakeTree(gAlice->GetTreeDFile());
1721   } else {
1722     arr->MakeTree(fDigitsFile);
1723   }
1724   SetDigitsArray(arr);
1725
1726   fDigitsSwitch=0; // standard digits
1727
1728   cerr<<"Digitizing TPC -- normal digits...\n";
1729
1730   for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) if (IsSectorActive(isec)) Hits2DigitsSector(isec);
1731
1732   // write results
1733
1734   char treeName[100];
1735
1736   sprintf(treeName,"TreeD_%s_%d",fTPCParam->GetTitle(),eventnumber);
1737
1738   GetDigitsArray()->GetTree()->SetName(treeName);
1739   GetDigitsArray()->GetTree()->AutoSave();
1740
1741
1742 }
1743
1744
1745
1746 //__________________________________________________________________
1747 void AliTPC::Hits2SDigits2(Int_t eventnumber)
1748 {
1749
1750   //-----------------------------------------------------------
1751   //   summable digits - 16 bit "ADC", no noise, no saturation
1752   //-----------------------------------------------------------
1753
1754  //----------------------------------------------------
1755  // Loop over all sectors for a single event
1756  //----------------------------------------------------
1757
1758
1759   if(fDefaults == 0) SetDefaults();
1760   GenerNoise(500000); //create table with noise
1761   //setup TPCDigitsArray
1762
1763   if(GetDigitsArray()) delete GetDigitsArray();
1764
1765   AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1766   arr->SetClass("AliSimDigits");
1767   arr->Setup(fTPCParam);
1768 // Note that methods arr->MakeTree have different signatures
1769   if (gAlice->GetTreeSFile()) {
1770     arr->MakeTree(gAlice->GetTreeSFile());
1771   } else {
1772     arr->MakeTree(fDigitsFile);
1773   }
1774   SetDigitsArray(arr);
1775
1776   cerr<<"Digitizing TPC -- summable digits...\n";
1777
1778   fDigitsSwitch=1; // summable digits
1779
1780     // set zero suppression to "0"
1781
1782   fTPCParam->SetZeroSup(0);
1783
1784  for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) if (IsSectorActive(isec)) Hits2DigitsSector(isec);
1785
1786
1787   // write results
1788
1789   char treeName[100];
1790
1791   sprintf(treeName,"TreeS_%s_%d",fTPCParam->GetTitle(),eventnumber);
1792
1793   GetDigitsArray()->GetTree()->SetName(treeName);
1794   GetDigitsArray()->GetTree()->AutoSave();
1795
1796 }
1797
1798
1799
1800 //__________________________________________________________________
1801 void AliTPC::Hits2SDigits()
1802 {
1803
1804   //-----------------------------------------------------------
1805   //   summable digits - 16 bit "ADC", no noise, no saturation
1806   //-----------------------------------------------------------
1807
1808  //----------------------------------------------------
1809  // Loop over all sectors for a single event
1810  //----------------------------------------------------
1811   //MI change - for pp run
1812   Int_t eventnumber = gAlice->GetEvNumber();
1813
1814   if(fDefaults == 0) SetDefaults();
1815   GenerNoise(500000); //create table with noise
1816
1817   //setup TPCDigitsArray
1818
1819   if(GetDigitsArray()) delete GetDigitsArray();
1820
1821   AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1822   arr->SetClass("AliSimDigits");
1823   arr->Setup(fTPCParam);
1824 // Note that methods arr->MakeTree have different signatures
1825   if (gAlice->GetTreeSFile()) {
1826     arr->MakeTree(gAlice->GetTreeSFile());
1827   } else {
1828     arr->MakeTree(fDigitsFile);
1829   }
1830   SetDigitsArray(arr);
1831
1832   cerr<<"Digitizing TPC -- summable digits...\n";
1833
1834   //  fDigitsSwitch=1; // summable digits  -for the moment direct
1835
1836   for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) if (IsSectorActive(isec)) Hits2DigitsSector(isec);
1837
1838
1839   // write results
1840   char treeName[100];
1841
1842   sprintf(treeName,"TreeD_%s_%d",fTPCParam->GetTitle(),eventnumber);
1843
1844   GetDigitsArray()->GetTree()->SetName(treeName);
1845   GetDigitsArray()->GetTree()->AutoSave();
1846
1847 }
1848
1849
1850 //_____________________________________________________________________________
1851 void AliTPC::Hits2DigitsSector(Int_t isec)
1852 {
1853   //-------------------------------------------------------------------
1854   // TPC conversion from hits to digits.
1855   //-------------------------------------------------------------------
1856
1857   //-----------------------------------------------------------------
1858   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1859   //-----------------------------------------------------------------
1860
1861   //-------------------------------------------------------
1862   //  Get the access to the track hits
1863   //-------------------------------------------------------
1864
1865   // check if the parameters are set - important if one calls this method
1866   // directly, not from the Hits2Digits
1867
1868   if(fDefaults == 0) SetDefaults();
1869
1870   TTree *tH = gAlice->TreeH(); // pointer to the hits tree
1871   Stat_t ntracks = tH->GetEntries();
1872
1873   if( ntracks > 0){
1874
1875   //-------------------------------------------
1876   //  Only if there are any tracks...
1877   //-------------------------------------------
1878
1879     TObjArray **row;
1880
1881     //printf("*** Processing sector number %d ***\n",isec);
1882
1883       Int_t nrows =fTPCParam->GetNRow(isec);
1884
1885       row= new TObjArray* [nrows+2]; // 2 extra rows for cross talk
1886
1887       MakeSector(isec,nrows,tH,ntracks,row);
1888
1889       //--------------------------------------------------------
1890       //   Digitize this sector, row by row
1891       //   row[i] is the pointer to the TObjArray of AliTPCFastVectors,
1892       //   each one containing electrons accepted on this
1893       //   row, assigned into tracks
1894       //--------------------------------------------------------
1895
1896       Int_t i;
1897
1898       if (fDigitsArray->GetTree()==0) fDigitsArray->MakeTree(fDigitsFile);
1899
1900       for (i=0;i<nrows;i++){
1901
1902         AliDigits * dig = fDigitsArray->CreateRow(isec,i);
1903
1904         DigitizeRow(i,isec,row);
1905
1906         fDigitsArray->StoreRow(isec,i);
1907
1908         Int_t ndig = dig->GetDigitSize();
1909         if (gDebug > 10) printf("*** Sector, row, compressed digits %d %d %d ***\n",isec,i,ndig);
1910
1911         fDigitsArray->ClearRow(isec,i);
1912
1913
1914        } // end of the sector digitization
1915
1916       for(i=0;i<nrows+2;i++){
1917         row[i]->Delete();
1918         delete row[i];
1919       }
1920
1921        delete [] row; // delete the array of pointers to TObjArray-s
1922
1923   } // ntracks >0
1924
1925 } // end of Hits2DigitsSector
1926
1927
1928 //_____________________________________________________________________________
1929 void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows)
1930 {
1931   //-----------------------------------------------------------
1932   // Single row digitization, coupling from the neighbouring
1933   // rows taken into account
1934   //-----------------------------------------------------------
1935
1936   //-----------------------------------------------------------------
1937   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1938   // Modified: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
1939   //-----------------------------------------------------------------
1940
1941
1942   Float_t zerosup = fTPCParam->GetZeroSup();
1943   //  Int_t nrows =fTPCParam->GetNRow(isec);
1944   fCurrentIndex[1]= isec;
1945
1946
1947   Int_t nofPads = fTPCParam->GetNPads(isec,irow);
1948   Int_t nofTbins = fTPCParam->GetMaxTBin();
1949   Int_t indexRange[4];
1950   //
1951   //  Integrated signal for this row
1952   //  and a single track signal
1953   //
1954
1955   AliTPCFastMatrix *m1 = new AliTPCFastMatrix(0,nofPads,0,nofTbins); // integrated
1956   AliTPCFastMatrix *m2 = new AliTPCFastMatrix(0,nofPads,0,nofTbins); // single
1957   //
1958   AliTPCFastMatrix &total  = *m1;
1959
1960   //  Array of pointers to the label-signal list
1961
1962   Int_t nofDigits = nofPads*nofTbins; // number of digits for this row
1963   Float_t  **pList = new Float_t* [nofDigits];
1964
1965   Int_t lp;
1966   Int_t i1;
1967   for(lp=0;lp<nofDigits;lp++)pList[lp]=0; // set all pointers to NULL
1968   //
1969   //calculate signal
1970   //
1971   //Int_t row1 = TMath::Max(irow-fTPCParam->GetNCrossRows(),0);
1972   //Int_t row2 = TMath::Min(irow+fTPCParam->GetNCrossRows(),nrows-1);
1973   Int_t row1=irow;
1974   Int_t row2=irow+2;
1975   for (Int_t row= row1;row<=row2;row++){
1976     Int_t nTracks= rows[row]->GetEntries();
1977     for (i1=0;i1<nTracks;i1++){
1978       fCurrentIndex[2]= row;
1979       fCurrentIndex[3]=irow+1;
1980       if (row==irow+1){
1981         m2->Zero();  // clear single track signal matrix
1982         Float_t trackLabel = GetSignal(rows[row],i1,m2,m1,indexRange);
1983         GetList(trackLabel,nofPads,m2,indexRange,pList);
1984       }
1985       else   GetSignal(rows[row],i1,0,m1,indexRange);
1986     }
1987   }
1988
1989   Int_t tracks[3];
1990
1991   AliDigits *dig = fDigitsArray->GetRow(isec,irow);
1992   Int_t gi=-1;
1993   Float_t fzerosup = zerosup+0.5;
1994   for(Int_t it=0;it<nofTbins;it++){
1995     Float_t *pq = &(total.UncheckedAt(0,it));
1996     for(Int_t ip=0;ip<nofPads;ip++){
1997       gi++;
1998       Float_t q=*pq;
1999       pq++;
2000       if(fDigitsSwitch == 0){
2001         q+=GetNoise();
2002         if(q <=fzerosup) continue; // do not fill zeros
2003         q = TMath::Nint(q);
2004         if(q > fTPCParam->GetADCSat()) q = fTPCParam->GetADCSat();  // saturation
2005
2006       }
2007
2008       else {
2009        if(q <= 0.) continue; // do not fill zeros
2010        if(q>2000.) q=2000.;
2011        q *= 16.;
2012        q = TMath::Nint(q);
2013       }
2014
2015       //
2016       //  "real" signal or electronic noise (list = -1)?
2017       //
2018
2019       for(Int_t j1=0;j1<3;j1++){
2020         tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -2;
2021       }
2022
2023 //Begin_Html
2024 /*
2025   <A NAME="AliDigits"></A>
2026   using of AliDigits object
2027 */
2028 //End_Html
2029       dig->SetDigitFast((Short_t)q,it,ip);
2030       if (fDigitsArray->IsSimulated())
2031         {
2032          ((AliSimDigits*)dig)->SetTrackIDFast(tracks[0],it,ip,0);
2033          ((AliSimDigits*)dig)->SetTrackIDFast(tracks[1],it,ip,1);
2034          ((AliSimDigits*)dig)->SetTrackIDFast(tracks[2],it,ip,2);
2035         }
2036
2037
2038     } // end of loop over time buckets
2039   }  // end of lop over pads
2040
2041   //
2042   //  This row has been digitized, delete nonused stuff
2043   //
2044
2045   for(lp=0;lp<nofDigits;lp++){
2046     if(pList[lp]) delete [] pList[lp];
2047   }
2048
2049   delete [] pList;
2050
2051   delete m1;
2052   delete m2;
2053   //  delete m3;
2054
2055 } // end of DigitizeRow
2056
2057 //_____________________________________________________________________________
2058
2059 Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr,
2060              AliTPCFastMatrix *m1, AliTPCFastMatrix *m2,Int_t *indexRange)
2061 {
2062
2063   //---------------------------------------------------------------
2064   //  Calculates 2-D signal (pad,time) for a single track,
2065   //  returns a pointer to the signal matrix and the track label
2066   //  No digitization is performed at this level!!!
2067   //---------------------------------------------------------------
2068
2069   //-----------------------------------------------------------------
2070   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2071   // Modified: Marian Ivanov
2072   //-----------------------------------------------------------------
2073
2074   AliTPCFastVector *tv;
2075
2076   tv = (AliTPCFastVector*)p1->At(ntr); // pointer to a track
2077   AliTPCFastVector &v = *tv;
2078
2079   Float_t label = v(0);
2080   Int_t centralPad = (fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3]-1)-1)/2;
2081
2082   Int_t nElectrons = (tv->GetNrows()-1)/4;
2083   indexRange[0]=9999; // min pad
2084   indexRange[1]=-1; // max pad
2085   indexRange[2]=9999; //min time
2086   indexRange[3]=-1; // max time
2087
2088   AliTPCFastMatrix &signal = *m1;
2089   AliTPCFastMatrix &total = *m2;
2090   //
2091   //  Loop over all electrons
2092   //
2093   for(Int_t nel=0; nel<nElectrons; nel++){
2094     Int_t idx=nel*4;
2095     Float_t aval =  v(idx+4);
2096     Float_t eltoadcfac=aval*fTPCParam->GetTotalNormFac();
2097     Float_t xyz[3]={v(idx+1),v(idx+2),v(idx+3)};
2098     Int_t n = ((AliTPCParamSR*)fTPCParam)->CalcResponseFast(xyz,fCurrentIndex,fCurrentIndex[3]);
2099
2100     Int_t *index = fTPCParam->GetResBin(0);
2101     Float_t *weight = & (fTPCParam->GetResWeight(0));
2102
2103     if (n>0) for (Int_t i =0; i<n; i++){
2104        Int_t pad=index[1]+centralPad;  //in digit coordinates central pad has coordinate 0
2105
2106          if (pad>=0){
2107          Int_t time=index[2];
2108          Float_t qweight = *(weight)*eltoadcfac;
2109
2110          if (m1!=0) signal.UncheckedAt(pad,time)+=qweight;
2111          total.UncheckedAt(pad,time)+=qweight;
2112          if (indexRange[0]>pad) indexRange[0]=pad;
2113          if (indexRange[1]<pad) indexRange[1]=pad;
2114          if (indexRange[2]>time) indexRange[2]=time;
2115          if (indexRange[3]<time) indexRange[3]=time;
2116
2117          index+=3;
2118          weight++;
2119
2120        }
2121     }
2122   } // end of loop over electrons
2123
2124   return label; // returns track label when finished
2125 }
2126
2127 //_____________________________________________________________________________
2128 void AliTPC::GetList(Float_t label,Int_t np,AliTPCFastMatrix *m,
2129                      Int_t *indexRange, Float_t **pList)
2130 {
2131   //----------------------------------------------------------------------
2132   //  Updates the list of tracks contributing to digits for a given row
2133   //----------------------------------------------------------------------
2134
2135   //-----------------------------------------------------------------
2136   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2137   //-----------------------------------------------------------------
2138
2139   AliTPCFastMatrix &signal = *m;
2140
2141   // lop over nonzero digits
2142
2143   for(Int_t it=indexRange[2];it<indexRange[3]+1;it++){
2144     for(Int_t ip=indexRange[0];ip<indexRange[1]+1;ip++){
2145
2146
2147         // accept only the contribution larger than 500 electrons (1/2 s_noise)
2148
2149         if(signal(ip,it)<0.5) continue;
2150
2151
2152         Int_t globalIndex = it*np+ip; // globalIndex starts from 0!
2153
2154         if(!pList[globalIndex]){
2155
2156           //
2157           // Create new list (6 elements - 3 signals and 3 labels),
2158           //
2159
2160           pList[globalIndex] = new Float_t [6];
2161
2162           // set list to -1
2163
2164           *pList[globalIndex] = -1.;
2165           *(pList[globalIndex]+1) = -1.;
2166           *(pList[globalIndex]+2) = -1.;
2167           *(pList[globalIndex]+3) = -1.;
2168           *(pList[globalIndex]+4) = -1.;
2169           *(pList[globalIndex]+5) = -1.;
2170
2171
2172           *pList[globalIndex] = label;
2173           *(pList[globalIndex]+3) = signal(ip,it);
2174         }
2175         else{
2176
2177           // check the signal magnitude
2178
2179           Float_t highest = *(pList[globalIndex]+3);
2180           Float_t middle = *(pList[globalIndex]+4);
2181           Float_t lowest = *(pList[globalIndex]+5);
2182
2183           //
2184           //  compare the new signal with already existing list
2185           //
2186
2187           if(signal(ip,it)<lowest) continue; // neglect this track
2188
2189           //
2190
2191           if (signal(ip,it)>highest){
2192             *(pList[globalIndex]+5) = middle;
2193             *(pList[globalIndex]+4) = highest;
2194             *(pList[globalIndex]+3) = signal(ip,it);
2195
2196             *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
2197             *(pList[globalIndex]+1) = *pList[globalIndex];
2198             *pList[globalIndex] = label;
2199           }
2200           else if (signal(ip,it)>middle){
2201             *(pList[globalIndex]+5) = middle;
2202             *(pList[globalIndex]+4) = signal(ip,it);
2203
2204             *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
2205             *(pList[globalIndex]+1) = label;
2206           }
2207           else{
2208             *(pList[globalIndex]+5) = signal(ip,it);
2209             *(pList[globalIndex]+2) = label;
2210           }
2211         }
2212
2213     } // end of loop over pads
2214   } // end of loop over time bins
2215
2216
2217
2218 }//end of GetList
2219 //___________________________________________________________________
2220 void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
2221                         Stat_t ntracks,TObjArray **row)
2222 {
2223
2224   //-----------------------------------------------------------------
2225   // Prepares the sector digitization, creates the vectors of
2226   // tracks for each row of this sector. The track vector
2227   // contains the track label and the position of electrons.
2228   //-----------------------------------------------------------------
2229
2230   //-----------------------------------------------------------------
2231   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2232   //-----------------------------------------------------------------
2233
2234   Float_t gasgain = fTPCParam->GetGasGain();
2235   Int_t i;
2236   Float_t xyz[4];
2237
2238   AliTPChit *tpcHit; // pointer to a sigle TPC hit
2239   //MI change
2240   TBranch * branch=0;
2241   if (fHitType>1) branch = TH->GetBranch("TPC2");
2242   else branch = TH->GetBranch("TPC");
2243
2244
2245   //----------------------------------------------
2246   // Create TObjArray-s, one for each row,
2247   // each TObjArray will store the AliTPCFastVectors
2248   // of electrons, one AliTPCFastVectors per each track.
2249   //----------------------------------------------
2250
2251   Int_t *nofElectrons = new Int_t [nrows+2]; // electron counter for each row
2252   AliTPCFastVector **tracks = new AliTPCFastVector* [nrows+2]; //pointers to the track vectors
2253
2254   for(i=0; i<nrows+2; i++){
2255     row[i] = new TObjArray;
2256     nofElectrons[i]=0;
2257     tracks[i]=0;
2258   }
2259
2260
2261
2262   //--------------------------------------------------------------------
2263   //  Loop over tracks, the "track" contains the full history
2264   //--------------------------------------------------------------------
2265
2266   Int_t previousTrack,currentTrack;
2267   previousTrack = -1; // nothing to store so far!
2268
2269   for(Int_t track=0;track<ntracks;track++){
2270     Bool_t isInSector=kTRUE;
2271     ResetHits();
2272     isInSector = TrackInVolume(isec,track);
2273     if (!isInSector) continue;
2274     //MI change
2275     branch->GetEntry(track); // get next track
2276
2277     //M.I. changes
2278
2279     tpcHit = (AliTPChit*)FirstHit(-1);
2280
2281     //--------------------------------------------------------------
2282     //  Loop over hits
2283     //--------------------------------------------------------------
2284
2285
2286     while(tpcHit){
2287
2288       Int_t sector=tpcHit->fSector; // sector number
2289       if(sector != isec){
2290         tpcHit = (AliTPChit*) NextHit();
2291         continue;
2292       }
2293
2294         currentTrack = tpcHit->Track(); // track number
2295
2296
2297         if(currentTrack != previousTrack){
2298
2299            // store already filled fTrack
2300
2301            for(i=0;i<nrows+2;i++){
2302              if(previousTrack != -1){
2303                if(nofElectrons[i]>0){
2304                  AliTPCFastVector &v = *tracks[i];
2305                  v(0) = previousTrack;
2306                  tracks[i]->ResizeTo(4*nofElectrons[i]+1); // shrink if necessary
2307                  row[i]->Add(tracks[i]);
2308                }
2309                else{
2310                  delete tracks[i]; // delete empty AliTPCFastVector
2311                  tracks[i]=0;
2312                }
2313              }
2314
2315              nofElectrons[i]=0;
2316              tracks[i] = new AliTPCFastVector(481); // AliTPCFastVectors for the next fTrack
2317
2318            } // end of loop over rows
2319
2320            previousTrack=currentTrack; // update track label
2321         }
2322
2323         Int_t qI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons)
2324
2325        //---------------------------------------------------
2326        //  Calculate the electron attachment probability
2327        //---------------------------------------------------
2328
2329
2330         Float_t time = 1.e6*(fTPCParam->GetZLength()-TMath::Abs(tpcHit->Z()))
2331                                                         /fTPCParam->GetDriftV();
2332         // in microseconds!
2333         Float_t attProb = fTPCParam->GetAttCoef()*
2334           fTPCParam->GetOxyCont()*time; //  fraction!
2335
2336         //-----------------------------------------------
2337         //  Loop over electrons
2338         //-----------------------------------------------
2339         Int_t index[3];
2340         index[1]=isec;
2341         for(Int_t nel=0;nel<qI;nel++){
2342           // skip if electron lost due to the attachment
2343           if((gRandom->Rndm(0)) < attProb) continue; // electron lost!
2344           xyz[0]=tpcHit->X();
2345           xyz[1]=tpcHit->Y();
2346           xyz[2]=tpcHit->Z();
2347           //
2348           // protection for the nonphysical avalanche size (10**6 maximum)
2349           //
2350           Double_t rn=TMath::Max(gRandom->Rndm(0),1.93e-22);
2351           xyz[3]= (Float_t) (-gasgain*TMath::Log(rn));
2352           index[0]=1;
2353
2354           TransportElectron(xyz,index);
2355           Int_t rowNumber;
2356           fTPCParam->GetPadRow(xyz,index);
2357           // row 0 - cross talk from the innermost row
2358           // row fNRow+1 cross talk from the outermost row
2359           rowNumber = index[2]+1;
2360           //transform position to local digit coordinates
2361           //relative to nearest pad row
2362           if ((rowNumber<0)||rowNumber>fTPCParam->GetNRow(isec)+1) continue;
2363           Float_t x1,y1;
2364           if (isec <fTPCParam->GetNInnerSector()) {
2365             x1 = xyz[1]*fTPCParam->GetInnerPadPitchWidth();
2366             y1 = fTPCParam->GetYInner(rowNumber);
2367           }
2368           else{
2369             x1=xyz[1]*fTPCParam->GetOuterPadPitchWidth();
2370             y1 = fTPCParam->GetYOuter(rowNumber);
2371           }
2372           // gain inefficiency at the wires edges - linear
2373           x1=TMath::Abs(x1);
2374           y1-=1.;
2375           if(x1>y1) xyz[3]*=TMath::Max(1.e-6,(y1-x1+1.));
2376
2377           nofElectrons[rowNumber]++;
2378           //----------------------------------
2379           // Expand vector if necessary
2380           //----------------------------------
2381           if(nofElectrons[rowNumber]>120){
2382             Int_t range = tracks[rowNumber]->GetNrows();
2383             if((nofElectrons[rowNumber])>(range-1)/4){
2384
2385               tracks[rowNumber]->ResizeTo(range+400); // Add 100 electrons
2386             }
2387           }
2388
2389           AliTPCFastVector &v = *tracks[rowNumber];
2390           Int_t idx = 4*nofElectrons[rowNumber]-3;
2391           Real_t * position = &(((AliTPCFastVector&)v).UncheckedAt(idx)); //make code faster
2392           memcpy(position,xyz,4*sizeof(Float_t));
2393
2394         } // end of loop over electrons
2395
2396         tpcHit = (AliTPChit*)NextHit();
2397
2398       } // end of loop over hits
2399     } // end of loop over tracks
2400
2401     //
2402     //   store remaining track (the last one) if not empty
2403     //
2404
2405      for(i=0;i<nrows+2;i++){
2406        if(nofElectrons[i]>0){
2407           AliTPCFastVector &v = *tracks[i];
2408           v(0) = previousTrack;
2409           tracks[i]->ResizeTo(4*nofElectrons[i]+1); // shrink if necessary
2410           row[i]->Add(tracks[i]);
2411         }
2412         else{
2413           delete tracks[i];
2414           tracks[i]=0;
2415         }
2416       }
2417
2418           delete [] tracks;
2419           delete [] nofElectrons;
2420
2421
2422 } // end of MakeSector
2423
2424
2425 //_____________________________________________________________________________
2426 void AliTPC::Init()
2427 {
2428   //
2429   // Initialise TPC detector after definition of geometry
2430   //
2431   Int_t i;
2432   //
2433   if(fDebug) {
2434     printf("\n%s: ",ClassName());
2435     for(i=0;i<35;i++) printf("*");
2436     printf(" TPC_INIT ");
2437     for(i=0;i<35;i++) printf("*");
2438     printf("\n%s: ",ClassName());
2439     //
2440     for(i=0;i<80;i++) printf("*");
2441     printf("\n");
2442   }
2443 }
2444
2445 //_____________________________________________________________________________
2446 void AliTPC::MakeBranch(Option_t* option, const char *file)
2447 {
2448   //
2449   // Create Tree branches for the TPC.
2450   //
2451   Int_t buffersize = 4000;
2452   char branchname[10];
2453   sprintf(branchname,"%s",GetName());
2454
2455   AliDetector::MakeBranch(option,file);
2456
2457   const char *d = strstr(option,"D");
2458
2459   if (fDigits   && gAlice->TreeD() && d) {
2460       MakeBranchInTree(gAlice->TreeD(),
2461                        branchname, &fDigits, buffersize, file);
2462   }
2463
2464   if (fHitType>1) MakeBranch2(option,file); // MI change 14.09.2000
2465 }
2466
2467 //_____________________________________________________________________________
2468 void AliTPC::ResetDigits()
2469 {
2470   //
2471   // Reset number of digits and the digits array for this detector
2472   //
2473   fNdigits   = 0;
2474   if (fDigits)   fDigits->Clear();
2475 }
2476
2477 //_____________________________________________________________________________
2478 void AliTPC::SetSecAL(Int_t sec)
2479 {
2480   //---------------------------------------------------
2481   // Activate/deactivate selection for lower sectors
2482   //---------------------------------------------------
2483
2484   //-----------------------------------------------------------------
2485   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2486   //-----------------------------------------------------------------
2487
2488   fSecAL = sec;
2489 }
2490
2491 //_____________________________________________________________________________
2492 void AliTPC::SetSecAU(Int_t sec)
2493 {
2494   //----------------------------------------------------
2495   // Activate/deactivate selection for upper sectors
2496   //---------------------------------------------------
2497
2498   //-----------------------------------------------------------------
2499   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2500   //-----------------------------------------------------------------
2501
2502   fSecAU = sec;
2503 }
2504
2505 //_____________________________________________________________________________
2506 void AliTPC::SetSecLows(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6)
2507 {
2508   //----------------------------------------
2509   // Select active lower sectors
2510   //----------------------------------------
2511
2512   //-----------------------------------------------------------------
2513   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2514   //-----------------------------------------------------------------
2515
2516   fSecLows[0] = s1;
2517   fSecLows[1] = s2;
2518   fSecLows[2] = s3;
2519   fSecLows[3] = s4;
2520   fSecLows[4] = s5;
2521   fSecLows[5] = s6;
2522 }
2523
2524 //_____________________________________________________________________________
2525 void AliTPC::SetSecUps(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6,
2526                        Int_t s7, Int_t s8 ,Int_t s9 ,Int_t s10,
2527                        Int_t s11 , Int_t s12)
2528 {
2529   //--------------------------------
2530   // Select active upper sectors
2531   //--------------------------------
2532
2533   //-----------------------------------------------------------------
2534   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2535   //-----------------------------------------------------------------
2536
2537   fSecUps[0] = s1;
2538   fSecUps[1] = s2;
2539   fSecUps[2] = s3;
2540   fSecUps[3] = s4;
2541   fSecUps[4] = s5;
2542   fSecUps[5] = s6;
2543   fSecUps[6] = s7;
2544   fSecUps[7] = s8;
2545   fSecUps[8] = s9;
2546   fSecUps[9] = s10;
2547   fSecUps[10] = s11;
2548   fSecUps[11] = s12;
2549 }
2550
2551 //_____________________________________________________________________________
2552 void AliTPC::SetSens(Int_t sens)
2553 {
2554
2555   //-------------------------------------------------------------
2556   // Activates/deactivates the sensitive strips at the center of
2557   // the pad row -- this is for the space-point resolution calculations
2558   //-------------------------------------------------------------
2559
2560   //-----------------------------------------------------------------
2561   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2562   //-----------------------------------------------------------------
2563
2564   fSens = sens;
2565 }
2566
2567
2568 void AliTPC::SetSide(Float_t side=0.)
2569 {
2570   // choice of the TPC side
2571
2572   fSide = side;
2573
2574 }
2575 //____________________________________________________________________________
2576 void AliTPC::SetGasMixt(Int_t nc,Int_t c1,Int_t c2,Int_t c3,Float_t p1,
2577                            Float_t p2,Float_t p3)
2578 {
2579
2580   // gax mixture definition
2581
2582  fNoComp = nc;
2583
2584  fMixtComp[0]=c1;
2585  fMixtComp[1]=c2;
2586  fMixtComp[2]=c3;
2587
2588  fMixtProp[0]=p1;
2589  fMixtProp[1]=p2;
2590  fMixtProp[2]=p3;
2591
2592
2593 }
2594 //_____________________________________________________________________________
2595
2596 void AliTPC::TransportElectron(Float_t *xyz, Int_t *index)
2597 {
2598   //
2599   // electron transport taking into account:
2600   // 1. diffusion,
2601   // 2.ExB at the wires
2602   // 3. nonisochronity
2603   //
2604   // xyz and index must be already transformed to system 1
2605   //
2606
2607   fTPCParam->Transform1to2(xyz,index);
2608
2609   //add diffusion
2610   Float_t driftl=xyz[2];
2611   if(driftl<0.01) driftl=0.01;
2612   driftl=TMath::Sqrt(driftl);
2613   Float_t sigT = driftl*(fTPCParam->GetDiffT());
2614   Float_t sigL = driftl*(fTPCParam->GetDiffL());
2615   xyz[0]=gRandom->Gaus(xyz[0],sigT);
2616   xyz[1]=gRandom->Gaus(xyz[1],sigT);
2617   xyz[2]=gRandom->Gaus(xyz[2],sigL);
2618
2619   // ExB
2620
2621   if (fTPCParam->GetMWPCReadout()==kTRUE){
2622     Float_t dx = fTPCParam->Transform2to2NearestWire(xyz,index);
2623     xyz[1]+=dx*(fTPCParam->GetOmegaTau());
2624   }
2625   //add nonisochronity (not implemented yet)
2626 }
2627
2628 ClassImp(AliTPCdigit)
2629
2630 //_____________________________________________________________________________
2631 AliTPCdigit::AliTPCdigit(Int_t *tracks, Int_t *digits):
2632   AliDigit(tracks)
2633 {
2634   //
2635   // Creates a TPC digit object
2636   //
2637   fSector     = digits[0];
2638   fPadRow     = digits[1];
2639   fPad        = digits[2];
2640   fTime       = digits[3];
2641   fSignal     = digits[4];
2642 }
2643
2644
2645 ClassImp(AliTPChit)
2646
2647 //_____________________________________________________________________________
2648 AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
2649 AliHit(shunt,track)
2650 {
2651   //
2652   // Creates a TPC hit object
2653   //
2654   fSector     = vol[0];
2655   fPadRow     = vol[1];
2656   fX          = hits[0];
2657   fY          = hits[1];
2658   fZ          = hits[2];
2659   fQ          = hits[3];
2660 }
2661
2662
2663 //________________________________________________________________________
2664 // Additional code because of the AliTPCTrackHitsV2
2665
2666 void AliTPC::MakeBranch2(Option_t *option,const char *file)
2667 {
2668   //
2669   // Create a new branch in the current Root Tree
2670   // The branch of fHits is automatically split
2671   // MI change 14.09.2000
2672   if (fHitType<2) return;
2673   char branchname[10];
2674   sprintf(branchname,"%s2",GetName());
2675   //
2676   // Get the pointer to the header
2677   const char *cH = strstr(option,"H");
2678   //
2679   if (fTrackHits   && gAlice->TreeH() && cH && fHitType&4) {
2680     //    AliObjectBranch * branch = new AliObjectBranch(branchname,"AliTPCTrackHitsV2",&fTrackHits,
2681     //                                             gAlice->TreeH(),fBufferSize,99);
2682     //TBranch * branch =
2683     gAlice->TreeH()->Branch(branchname,"AliTPCTrackHitsV2",&fTrackHits,
2684                                                    fBufferSize,99);
2685
2686     // gAlice->TreeH()->GetListOfBranches()->Add(branch);
2687     if (GetDebug()>1)
2688       printf("* AliDetector::MakeBranch * Making Branch %s for trackhits\n",branchname);
2689     const char folder [] = "RunMC/Event/Data";
2690     if (GetDebug())
2691       printf("%15s: Publishing %s to %s\n",ClassName(),branchname,folder);
2692     Publish(folder,&fTrackHits,branchname);
2693     if (file) {
2694         TBranch *b = gAlice->TreeH()->GetBranch(branchname);
2695         TDirectory *wd = gDirectory;
2696         b->SetFile(file);
2697         TIter next( b->GetListOfBranches());
2698         while ((b=(TBranch*)next())) {
2699           b->SetFile(file);
2700         }
2701         wd->cd();
2702         if (GetDebug()>1)
2703               cout << "Diverting branch " << branchname << " to file " << file << endl;
2704     }
2705   }
2706
2707   if (fTrackHitsOld   && gAlice->TreeH() && cH && fHitType&2) {
2708     AliObjectBranch * branch = new AliObjectBranch(branchname,"AliTPCTrackHits",&fTrackHitsOld,
2709                                                    gAlice->TreeH(),fBufferSize,99);
2710     //TBranch * branch = gAlice->TreeH()->Branch(branchname,"AliTPCTrackHitsV2",&fTrackHits,
2711     //                                                     fBufferSize,99);
2712
2713     gAlice->TreeH()->GetListOfBranches()->Add(branch);
2714     if (GetDebug()>1)
2715       printf("* AliDetector::MakeBranch * Making Branch %s for trackhits\n",branchname);
2716     const char folder [] = "RunMC/Event/Data";
2717     if (GetDebug())
2718       printf("%15s: Publishing %s to %s\n",ClassName(),branchname,folder);
2719     Publish(folder,&fTrackHitsOld,branchname);
2720     if (file) {
2721         TBranch *b = gAlice->TreeH()->GetBranch(branchname);
2722         TDirectory *wd = gDirectory;
2723         b->SetFile(file);
2724         TIter next( b->GetListOfBranches());
2725         while ((b=(TBranch*)next())) {
2726           b->SetFile(file);
2727         }
2728         wd->cd();
2729         if (GetDebug()>1)
2730               cout << "Diverting branch " << branchname << " to file " << file << endl;
2731     }
2732   }
2733 }
2734
2735 void AliTPC::SetTreeAddress()
2736 {
2737   if (fHitType<=1) AliDetector::SetTreeAddress();
2738   if (fHitType>1) SetTreeAddress2();
2739 }
2740
2741 void AliTPC::SetTreeAddress2()
2742 {
2743   //
2744   // Set branch address for the TrackHits Tree
2745   //
2746   TBranch *branch;
2747   char branchname[20];
2748   sprintf(branchname,"%s2",GetName());
2749   //
2750   // Branch address for hit tree
2751   TTree *treeH = gAlice->TreeH();
2752   if ((treeH)&&(fHitType&4)) {
2753     branch = treeH->GetBranch(branchname);
2754     if (branch) branch->SetAddress(&fTrackHits);
2755   }
2756   if ((treeH)&&(fHitType&2)) {
2757     branch = treeH->GetBranch(branchname);
2758     if (branch) branch->SetAddress(&fTrackHitsOld);
2759   }
2760   //set address to TREETR
2761   TTree *treeTR = gAlice->TreeTR();
2762   if (treeTR && fTrackReferences) {
2763     branch = treeTR->GetBranch(GetName());
2764     if (branch) branch->SetAddress(&fTrackReferences);
2765   }
2766
2767
2768 }
2769
2770 void AliTPC::FinishPrimary()
2771 {
2772   if (fTrackHits &&fHitType&4)      fTrackHits->FlushHitStack();
2773   if (fTrackHitsOld && fHitType&2)  fTrackHitsOld->FlushHitStack();
2774 }
2775
2776
2777 void AliTPC::AddHit2(Int_t track, Int_t *vol, Float_t *hits)
2778 {
2779   //
2780   // add hit to the list
2781   Int_t rtrack;
2782   if (fIshunt) {
2783     int primary = gAlice->GetPrimary(track);
2784     gAlice->Particle(primary)->SetBit(kKeepBit);
2785     rtrack=primary;
2786   } else {
2787     rtrack=track;
2788     gAlice->FlagTrack(track);
2789   }
2790   //AliTPChit *hit = (AliTPChit*)fHits->UncheckedAt(fNhits-1);
2791   //if (hit->fTrack!=rtrack)
2792   //  cout<<"bad track number\n";
2793   if (fTrackHits && fHitType&4)
2794     fTrackHits->AddHitKartez(vol[0],rtrack, hits[0],
2795                              hits[1],hits[2],(Int_t)hits[3]);
2796   if (fTrackHitsOld &&fHitType&2 )
2797     fTrackHitsOld->AddHitKartez(vol[0],rtrack, hits[0],
2798                              hits[1],hits[2],(Int_t)hits[3]);
2799
2800 }
2801
2802 void AliTPC::ResetHits()
2803 {
2804   if (fHitType&1) AliDetector::ResetHits();
2805   if (fHitType>1) ResetHits2();
2806 }
2807
2808 void AliTPC::ResetHits2()
2809 {
2810   //
2811   //reset hits
2812   if (fTrackHits && fHitType&4) fTrackHits->Clear();
2813   if (fTrackHitsOld && fHitType&2) fTrackHitsOld->Clear();
2814
2815 }
2816
2817 AliHit* AliTPC::FirstHit(Int_t track)
2818 {
2819   if (fHitType>1) return FirstHit2(track);
2820   return AliDetector::FirstHit(track);
2821 }
2822 AliHit* AliTPC::NextHit()
2823 {
2824   if (fHitType>1) return NextHit2();
2825
2826   return AliDetector::NextHit();
2827 }
2828
2829 AliHit* AliTPC::FirstHit2(Int_t track)
2830 {
2831   //
2832   // Initialise the hit iterator
2833   // Return the address of the first hit for track
2834   // If track>=0 the track is read from disk
2835   // while if track<0 the first hit of the current
2836   // track is returned
2837   //
2838   if(track>=0) {
2839     gAlice->ResetHits();
2840     gAlice->TreeH()->GetEvent(track);
2841   }
2842   //
2843   if (fTrackHits && fHitType&4) {
2844     fTrackHits->First();
2845     return fTrackHits->GetHit();
2846   }
2847   if (fTrackHitsOld && fHitType&2) {
2848     fTrackHitsOld->First();
2849     return fTrackHitsOld->GetHit();
2850   }
2851
2852   else return 0;
2853 }
2854
2855 AliHit* AliTPC::NextHit2()
2856 {
2857   //
2858   //Return the next hit for the current track
2859
2860
2861   if (fTrackHitsOld && fHitType&2) {
2862     fTrackHitsOld->Next();
2863     return fTrackHitsOld->GetHit();
2864   }
2865   if (fTrackHits) {
2866     fTrackHits->Next();
2867     return fTrackHits->GetHit();
2868   }
2869   else
2870     return 0;
2871 }
2872
2873 void AliTPC::LoadPoints(Int_t)
2874 {
2875   //
2876   Int_t a = 0;
2877   /*  if(fHitType==1) return AliDetector::LoadPoints(a);
2878   LoadPoints2(a);
2879   */
2880   if(fHitType==1) AliDetector::LoadPoints(a);
2881   else LoadPoints2(a);
2882
2883   // LoadPoints3(a);
2884
2885 }
2886
2887
2888 void AliTPC::RemapTrackHitIDs(Int_t *map)
2889 {
2890   if (!fTrackHits) return;
2891
2892   if (fTrackHitsOld && fHitType&2){
2893     AliObjectArray * arr = fTrackHitsOld->fTrackHitsInfo;
2894     for (UInt_t i=0;i<arr->GetSize();i++){
2895       AliTrackHitsInfo * info = (AliTrackHitsInfo *)(arr->At(i));
2896       info->fTrackID = map[info->fTrackID];
2897     }
2898   }
2899   if (fTrackHitsOld && fHitType&4){
2900     TClonesArray * arr = fTrackHits->GetArray();;
2901     for (Int_t i=0;i<arr->GetEntriesFast();i++){
2902       AliTrackHitsParamV2 * info = (AliTrackHitsParamV2 *)(arr->At(i));
2903       info->fTrackID = map[info->fTrackID];
2904     }
2905   }
2906 }
2907
2908 Bool_t   AliTPC::TrackInVolume(Int_t id,Int_t track)
2909 {
2910   //return bool information - is track in given volume
2911   //load only part of the track information
2912   //return true if current track is in volume
2913   //
2914   //  return kTRUE;
2915   if (fTrackHitsOld && fHitType&2) {
2916     TBranch * br = gAlice->TreeH()->GetBranch("fTrackHitsInfo");
2917     br->GetEvent(track);
2918     AliObjectArray * ar = fTrackHitsOld->fTrackHitsInfo;
2919     for (UInt_t j=0;j<ar->GetSize();j++){
2920       if (  ((AliTrackHitsInfo*)ar->At(j))->fVolumeID==id) return kTRUE;
2921     }
2922   }
2923
2924   if (fTrackHits && fHitType&4) {
2925     TBranch * br1 = gAlice->TreeH()->GetBranch("fVolumes");
2926     TBranch * br2 = gAlice->TreeH()->GetBranch("fNVolumes");
2927     br2->GetEvent(track);
2928     br1->GetEvent(track);
2929     Int_t *volumes = fTrackHits->GetVolumes();
2930     Int_t nvolumes = fTrackHits->GetNVolumes();
2931     if (!volumes && nvolumes>0) {
2932       printf("Problematic track\t%d\t%d",track,nvolumes);
2933       return kFALSE;
2934     }
2935     for (Int_t j=0;j<nvolumes; j++)
2936       if (volumes[j]==id) return kTRUE;
2937   }
2938
2939   if (fHitType&1) {
2940     TBranch * br = gAlice->TreeH()->GetBranch("fSector");
2941     br->GetEvent(track);
2942     for (Int_t j=0;j<fHits->GetEntriesFast();j++){
2943       if (  ((AliTPChit*)fHits->At(j))->fSector==id) return kTRUE;
2944     }
2945   }
2946   return kFALSE;
2947
2948 }
2949
2950 //_____________________________________________________________________________
2951 void AliTPC::LoadPoints2(Int_t)
2952 {
2953   //
2954   // Store x, y, z of all hits in memory
2955   //
2956   if (fTrackHits == 0 && fTrackHitsOld==0) return;
2957   //
2958   Int_t nhits =0;
2959   if (fHitType&4) nhits = fTrackHits->GetEntriesFast();
2960   if (fHitType&2) nhits = fTrackHitsOld->GetEntriesFast();
2961
2962   if (nhits == 0) return;
2963   Int_t tracks = gAlice->GetNtrack();
2964   if (fPoints == 0) fPoints = new TObjArray(tracks);
2965   AliHit *ahit;
2966   //
2967   Int_t *ntrk=new Int_t[tracks];
2968   Int_t *limi=new Int_t[tracks];
2969   Float_t **coor=new Float_t*[tracks];
2970   for(Int_t i=0;i<tracks;i++) {
2971     ntrk[i]=0;
2972     coor[i]=0;
2973     limi[i]=0;
2974   }
2975   //
2976   AliPoints *points = 0;
2977   Float_t *fp=0;
2978   Int_t trk;
2979   Int_t chunk=nhits/4+1;
2980   //
2981   // Loop over all the hits and store their position
2982   //
2983   ahit = FirstHit2(-1);
2984   while (ahit){
2985     trk=ahit->GetTrack();
2986     if(ntrk[trk]==limi[trk]) {
2987       //
2988       // Initialise a new track
2989       fp=new Float_t[3*(limi[trk]+chunk)];
2990       if(coor[trk]) {
2991         memcpy(fp,coor[trk],sizeof(Float_t)*3*limi[trk]);
2992         delete [] coor[trk];
2993       }
2994       limi[trk]+=chunk;
2995       coor[trk] = fp;
2996     } else {
2997       fp = coor[trk];
2998     }
2999     fp[3*ntrk[trk]  ] = ahit->X();
3000     fp[3*ntrk[trk]+1] = ahit->Y();
3001     fp[3*ntrk[trk]+2] = ahit->Z();
3002     ntrk[trk]++;
3003     ahit = NextHit2();
3004   }
3005
3006
3007
3008   //
3009   for(trk=0; trk<tracks; ++trk) {
3010     if(ntrk[trk]) {
3011       points = new AliPoints();
3012       points->SetMarkerColor(GetMarkerColor());
3013       points->SetMarkerSize(GetMarkerSize());
3014       points->SetDetector(this);
3015       points->SetParticle(trk);
3016       points->SetPolyMarker(ntrk[trk],coor[trk],GetMarkerStyle());
3017       fPoints->AddAt(points,trk);
3018       delete [] coor[trk];
3019       coor[trk]=0;
3020     }
3021   }
3022   delete [] coor;
3023   delete [] ntrk;
3024   delete [] limi;
3025 }
3026
3027
3028 //_____________________________________________________________________________
3029 void AliTPC::LoadPoints3(Int_t)
3030 {
3031   //
3032   // Store x, y, z of all hits in memory
3033   // - only intersection point with pad row
3034   if (fTrackHits == 0) return;
3035   //
3036   Int_t nhits = fTrackHits->GetEntriesFast();
3037   if (nhits == 0) return;
3038   Int_t tracks = gAlice->GetNtrack();
3039   if (fPoints == 0) fPoints = new TObjArray(2*tracks);
3040   fPoints->Expand(2*tracks);
3041   AliHit *ahit;
3042   //
3043   Int_t *ntrk=new Int_t[tracks];
3044   Int_t *limi=new Int_t[tracks];
3045   Float_t **coor=new Float_t*[tracks];
3046   for(Int_t i=0;i<tracks;i++) {
3047     ntrk[i]=0;
3048     coor[i]=0;
3049     limi[i]=0;
3050   }
3051   //
3052   AliPoints *points = 0;
3053   Float_t *fp=0;
3054   Int_t trk;
3055   Int_t chunk=nhits/4+1;
3056   //
3057   // Loop over all the hits and store their position
3058   //
3059   ahit = FirstHit2(-1);
3060   //for (Int_t hit=0;hit<nhits;hit++) {
3061
3062   Int_t lastrow = -1;
3063   while (ahit){
3064     //    ahit = (AliHit*)fHits->UncheckedAt(hit);
3065     trk=ahit->GetTrack();
3066     Float_t  x[3]={ahit->X(),ahit->Y(),ahit->Z()};
3067     Int_t    index[3]={1,((AliTPChit*)ahit)->fSector,0};
3068     Int_t    currentrow = fTPCParam->GetPadRow(x,index) ;
3069     if (currentrow!=lastrow){
3070       lastrow = currentrow;
3071       //later calculate intersection point
3072       if(ntrk[trk]==limi[trk]) {
3073         //
3074         // Initialise a new track
3075         fp=new Float_t[3*(limi[trk]+chunk)];
3076         if(coor[trk]) {
3077           memcpy(fp,coor[trk],sizeof(Float_t)*3*limi[trk]);
3078           delete [] coor[trk];
3079         }
3080         limi[trk]+=chunk;
3081         coor[trk] = fp;
3082       } else {
3083         fp = coor[trk];
3084       }
3085       fp[3*ntrk[trk]  ] = ahit->X();
3086       fp[3*ntrk[trk]+1] = ahit->Y();
3087       fp[3*ntrk[trk]+2] = ahit->Z();
3088       ntrk[trk]++;
3089     }
3090     ahit = NextHit2();
3091   }
3092
3093   //
3094   for(trk=0; trk<tracks; ++trk) {
3095     if(ntrk[trk]) {
3096       points = new AliPoints();
3097       points->SetMarkerColor(GetMarkerColor()+1);
3098       points->SetMarkerStyle(5);
3099       points->SetMarkerSize(0.2);
3100       points->SetDetector(this);
3101       points->SetParticle(trk);
3102       //      points->SetPolyMarker(ntrk[trk],coor[trk],GetMarkerStyle()20);
3103       points->SetPolyMarker(ntrk[trk],coor[trk],30);
3104       fPoints->AddAt(points,tracks+trk);
3105       delete [] coor[trk];
3106       coor[trk]=0;
3107     }
3108   }
3109   delete [] coor;
3110   delete [] ntrk;
3111   delete [] limi;
3112 }
3113
3114
3115
3116 void AliTPC::FindTrackHitsIntersection(TClonesArray * arr)
3117 {
3118
3119   //
3120   //fill clones array with intersection of current point with the
3121   //middle of the row
3122   Int_t sector;
3123   Int_t ipart;
3124
3125   const Int_t kcmaxhits=30000;
3126   AliTPCFastVector * xxxx = new AliTPCFastVector(kcmaxhits*4);
3127   AliTPCFastVector & xxx = *xxxx;
3128   Int_t maxhits = kcmaxhits;
3129
3130   //
3131   AliTPChit * tpcHit=0;
3132   tpcHit = (AliTPChit*)FirstHit2(-1);
3133   Int_t currentIndex=0;
3134   Int_t lastrow=-1;  //last writen row
3135
3136   while (tpcHit){
3137     if (tpcHit==0) continue;
3138     sector=tpcHit->fSector; // sector number
3139     ipart=tpcHit->Track();
3140
3141     //find row number
3142
3143     Float_t  x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
3144     Int_t    index[3]={1,sector,0};
3145     Int_t    currentrow = fTPCParam->GetPadRow(x,index) ;
3146     if (currentrow<0) continue;
3147     if (lastrow<0) lastrow=currentrow;
3148     if (currentrow==lastrow){
3149       if ( currentIndex>=maxhits){
3150         maxhits+=kcmaxhits;
3151         xxx.ResizeTo(4*maxhits);
3152       }
3153       xxx(currentIndex*4)=x[0];
3154       xxx(currentIndex*4+1)=x[1];
3155       xxx(currentIndex*4+2)=x[2];
3156       xxx(currentIndex*4+3)=tpcHit->fQ;
3157       currentIndex++;
3158     }
3159     else
3160       if (currentIndex>2){
3161         Float_t sumx=0;
3162         Float_t sumx2=0;
3163         Float_t sumx3=0;
3164         Float_t sumx4=0;
3165         Float_t sumy=0;
3166         Float_t sumxy=0;
3167         Float_t sumx2y=0;
3168         Float_t sumz=0;
3169         Float_t sumxz=0;
3170         Float_t sumx2z=0;
3171         Float_t sumq=0;
3172         for (Int_t index=0;index<currentIndex;index++){
3173           Float_t x,x2,x3,x4;
3174           x=x2=x3=x4=xxx(index*4);
3175           x2*=x;
3176           x3*=x2;
3177           x4*=x3;
3178           sumx+=x;
3179           sumx2+=x2;
3180           sumx3+=x3;
3181           sumx4+=x4;
3182           sumy+=xxx(index*4+1);
3183           sumxy+=xxx(index*4+1)*x;
3184           sumx2y+=xxx(index*4+1)*x2;
3185           sumz+=xxx(index*4+2);
3186           sumxz+=xxx(index*4+2)*x;
3187           sumx2z+=xxx(index*4+2)*x2;
3188           sumq+=xxx(index*4+3);
3189         }
3190         Float_t centralPad = (fTPCParam->GetNPads(sector,lastrow)-1)/2;
3191         Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
3192           sumx2*(sumx*sumx3-sumx2*sumx2);
3193
3194         Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
3195           sumx2*(sumxy*sumx3-sumx2y*sumx2);
3196         Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
3197           sumx2*(sumxz*sumx3-sumx2z*sumx2);
3198
3199         Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
3200           sumx2*(sumx*sumx2y-sumx2*sumxy);
3201         Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
3202           sumx2*(sumx*sumx2z-sumx2*sumxz);
3203
3204         Float_t y=detay/det+centralPad;
3205         Float_t z=detaz/det;
3206         Float_t by=detby/det; //y angle
3207         Float_t bz=detbz/det; //z angle
3208         sumy/=Float_t(currentIndex);
3209         sumz/=Float_t(currentIndex);
3210
3211         AliComplexCluster cl;
3212         cl.fX=z;
3213         cl.fY=y;
3214         cl.fQ=sumq;
3215         cl.fSigmaX2=bz;
3216         cl.fSigmaY2=by;
3217         cl.fTracks[0]=ipart;
3218
3219         AliTPCClustersRow * row = (fClustersArray->GetRow(sector,lastrow));
3220         if (row!=0) row->InsertCluster(&cl);
3221         currentIndex=0;
3222         lastrow=currentrow;
3223       } //end of calculating cluster for given row
3224
3225   } // end of loop over hits
3226   xxxx->Delete();
3227
3228 }
3229 //_______________________________________________________________________________
3230 void AliTPC::Digits2Reco(Int_t firstevent,Int_t lastevent)
3231 {
3232   // produces rec points from digits and writes them on the root file
3233   // AliTPCclusters.root
3234
3235   TDirectory *cwd = gDirectory;
3236
3237
3238   AliTPCParamSR *dig=(AliTPCParamSR *)gDirectory->Get("75x40_100x60");
3239   if(dig){
3240     printf("You are running 2 pad-length geom hits with 3 pad-length geom digits\n");
3241     delete dig;
3242     dig = new AliTPCParamSR();
3243   }
3244   else
3245   {
3246    dig=(AliTPCParamSR *)gDirectory->Get("75x40_100x60_150x60");
3247   }
3248   if(!dig){
3249    printf("No TPC parameters found\n");
3250    exit(3);
3251   }
3252
3253   SetParam(dig);
3254   cout<<"AliTPC::Digits2Reco: TPC parameteres have been set"<<endl;
3255   TFile *out;
3256   if(!gSystem->Getenv("CONFIG_FILE")){
3257     out=TFile::Open("AliTPCclusters.root","recreate");
3258   }
3259   else {
3260     const char *tmp1;
3261     const char *tmp2;
3262     char tmp3[80];
3263     tmp1=gSystem->Getenv("CONFIG_FILE_PREFIX");
3264     tmp2=gSystem->Getenv("CONFIG_OUTDIR");
3265     sprintf(tmp3,"%s%s/AliTPCclusters.root",tmp1,tmp2);
3266     out=TFile::Open(tmp3,"recreate");
3267   }
3268
3269   TStopwatch timer;
3270   cout<<"AliTPC::Digits2Reco - determination of rec points begins"<<endl;
3271   timer.Start();
3272   cwd->cd();
3273   for(Int_t iev=firstevent;iev<lastevent+1;iev++){
3274
3275     printf("Processing event %d\n",iev);
3276     Digits2Clusters(out,iev);
3277   }
3278   cout<<"AliTPC::Digits2Reco - determination of rec points ended"<<endl;
3279   timer.Stop();
3280   timer.Print();
3281   out->Close();
3282   cwd->cd();
3283
3284
3285 }