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