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.26  2001/01/10 07:59:43  kowal2
  19 Corrections to load points from the noncompressed hits.
  20
  21 Revision 1.25  2000/11/02 07:25:31  kowal2
  22 Changes due to the new hit structure.
  23 Memory leak removed.
  24
  25 Revision 1.24  2000/10/05 16:06:09  kowal2
  26 Forward declarations. Changes due to a new class AliComplexCluster.
  27
  28 Revision 1.23  2000/10/02 21:28:18  fca
  29 Removal of useless dependecies via forward declarations
  30
  31 Revision 1.22  2000/07/10 20:57:39  hristov
  32 Update of TPC code and macros by M.Kowalski
  33
  34 Revision 1.19.2.4  2000/06/26 07:39:42  kowal2
  35 Changes to obey the coding rules
  36
  37 Revision 1.19.2.3  2000/06/25 08:38:41  kowal2
  38 Splitted from AliTPCtracking
  39
  40 Revision 1.19.2.2  2000/06/16 12:59:28  kowal2
  41 Changed parameter settings
  42
  43 Revision 1.19.2.1  2000/06/09 07:15:07  kowal2
  44
  45 Defaults loaded automatically (hard-wired)
  46 Optional parameters can be set via macro called in the constructor
  47
  48 Revision 1.19  2000/04/18 19:00:59  fca
  49 Small bug fixes to TPC files
  50
  51 Revision 1.18  2000/04/17 09:37:33  kowal2
  52 removed obsolete AliTPCDigitsDisplay.C
  53
  54 Revision 1.17.2.2  2000/04/10 08:15:12  kowal2
  55
  56 New, experimental data structure from M. Ivanov
  57 New tracking algorithm
  58 Different pad geometry for different sectors
  59 Digitization rewritten
  60
  61 Revision 1.17.2.1  2000/04/10 07:56:53  kowal2
  62 Not used anymore - removed
  63
  64 Revision 1.17  2000/01/19 17:17:30  fca
  65 Introducing a list of lists of hits -- more hits allowed for detector now
  66
  67 Revision 1.16  1999/11/05 09:29:23  fca
  68 Accept only signals > 0
  69
  70 Revision 1.15  1999/10/08 06:26:53  fca
  71 Removed ClustersIndex - not used anymore
  72
  73 Revision 1.14  1999/09/29 09:24:33  fca
  74 Introduction of the Copyright and cvs Log
  75
  76 */
  77
  78 ///////////////////////////////////////////////////////////////////////////////
  79 //                                                                           //
  80 //  Time Projection Chamber                                                  //
  81 //  This class contains the basic functions for the Time Projection Chamber  //
  82 //  detector. Functions specific to one particular geometry are              //
  83 //  contained in the derived classes                                         //
  84 //                                                                           //
  85 //Begin_Html
  86 /*
  87 <img src="picts/AliTPCClass.gif">
  88 */
  89 //End_Html
  90 //                                                                           //
  91 //                                                                          //
  92 ///////////////////////////////////////////////////////////////////////////////
  93
  94 //
  95
  96 #include <TMath.h>
  97 #include <TRandom.h>
  98 #include <TVector.h>
  99 #include <TMatrix.h>
 100 #include <TGeometry.h>
 101 #include <TNode.h>
 102 #include <TTUBS.h>
 103 #include <TObjectTable.h>
 104 #include "TParticle.h"
 105 #include "AliTPC.h"
 106 #include <TFile.h>
 107 #include "AliRun.h"
 108 #include <iostream.h>
 109 #include <fstream.h>
 110 #include "AliMC.h"
 111 #include "AliMagF.h"
 112
 113
 114 #include "AliTPCParamSR.h"
 115 #include "AliTPCPRF2D.h"
 116 #include "AliTPCRF1D.h"
 117 #include "AliDigits.h"
 118 #include "AliSimDigits.h"
 119 #include "AliTPCTrackHits.h"
 120 #include "AliPoints.h"
 121 #include "AliArrayBranch.h"
 122
 123
 124 #include "AliTPCDigitsArray.h"
 125 #include "AliComplexCluster.h"
 126 #include "AliClusters.h"
 127 #include "AliTPCClustersRow.h"
 128 #include "AliTPCClustersArray.h"
 129
 130 #include "AliTPCcluster.h"
 131 #include "AliTPCclusterer.h"
 132 #include "AliTPCtracker.h"
 133
 134 #include <TInterpreter.h>
 135 #include <TTree.h>
 136
 137
 138
 139 ClassImp(AliTPC)
 140
 141 //_____________________________________________________________________________
 142 AliTPC::AliTPC()
 143 {
 144   //
 145   // Default constructor
 146   //
 147   fIshunt   = 0;
 148   fHits     = 0;
 149   fDigits   = 0;
 150   fNsectors = 0;
 151   //MI changes
 152   fDigitsArray = 0;
 153   fClustersArray = 0;
 154   fTPCParam=0;
 155   fTrackHits = 0;
 156   fHitType = 2;
 157   fTPCParam = 0;
 158 }
 159
 160 //_____________________________________________________________________________
 161 AliTPC::AliTPC(const char *name, const char *title)
 162       : AliDetector(name,title)
 163 {
 164   //
 165   // Standard constructor
 166   //
 167
 168   //
 169   // Initialise arrays of hits and digits
 170   fHits     = new TClonesArray("AliTPChit",  176);
 171   gAlice->AddHitList(fHits);
 172   //MI change
 173   fDigitsArray = 0;
 174   fClustersArray= 0;
 175   //
 176   fTrackHits = new AliTPCTrackHits;  //MI - 13.09.2000
 177   fTrackHits->SetHitPrecision(0.002);
 178   fTrackHits->SetStepPrecision(0.003);
 179   fTrackHits->SetMaxDistance(100);
 180
 181   fHitType = 2;
 182   //
 183   // Initialise counters
 184   fNsectors = 0;
 185
 186   //
 187   fIshunt     =  0;
 188   //
 189   // Initialise color attributes
 190   SetMarkerColor(kYellow);
 191
 192   //
 193   //  Set TPC parameters
 194   //
 195
 196   if (!strcmp(title,"Default")) {
 197      fTPCParam = new AliTPCParamSR;
 198   } else {
 199     cerr<<"AliTPC warning: in Config.C you must set non-default parameters\n";
 200     fTPCParam=0;
 201   }
 202
 203 }
 204
 205 //_____________________________________________________________________________
 206 AliTPC::~AliTPC()
 207 {
 208   //
 209   // TPC destructor
 210   //
 211
 212   fIshunt   = 0;
 213   delete fHits;
 214   delete fDigits;
 215   delete fTPCParam;
 216   delete fTrackHits; //MI 15.09.2000
 217 }
 218
 219 //_____________________________________________________________________________
 220 void AliTPC::AddHit(Int_t track, Int_t *vol, Float_t *hits)
 221 {
 222   //
 223   // Add a hit to the list
 224   //
 225   //  TClonesArray &lhits = *fHits;
 226   //  new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits);
 227   if (fHitType&1){
 228     TClonesArray &lhits = *fHits;
 229     new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits);
 230   }
 231   if (fHitType&2)
 232    AddHit2(track,vol,hits);
 233 }
 234
 235 //_____________________________________________________________________________
 236 void AliTPC::BuildGeometry()
 237 {
 238
 239   //
 240   // Build TPC ROOT TNode geometry for the event display
 241   //
 242   TNode *nNode, *nTop;
 243   TTUBS *tubs;
 244   Int_t i;
 245   const int kColorTPC=19;
 246   char name[5], title[25];
 247   const Double_t kDegrad=TMath::Pi()/180;
 248   const Double_t kRaddeg=180./TMath::Pi();
 249
 250
 251   Float_t innerOpenAngle = fTPCParam->GetInnerAngle();
 252   Float_t outerOpenAngle = fTPCParam->GetOuterAngle();
 253
 254   Float_t innerAngleShift = fTPCParam->GetInnerAngleShift();
 255   Float_t outerAngleShift = fTPCParam->GetOuterAngleShift();
 256
 257   Int_t nLo = fTPCParam->GetNInnerSector()/2;
 258   Int_t nHi = fTPCParam->GetNOuterSector()/2;
 259
 260   const Double_t kloAng = (Double_t)TMath::Nint(innerOpenAngle*kRaddeg);
 261   const Double_t khiAng = (Double_t)TMath::Nint(outerOpenAngle*kRaddeg);
 262   const Double_t kloAngSh = (Double_t)TMath::Nint(innerAngleShift*kRaddeg);
 263   const Double_t khiAngSh = (Double_t)TMath::Nint(outerAngleShift*kRaddeg);
 264
 265
 266   const Double_t kloCorr = 1/TMath::Cos(0.5*kloAng*kDegrad);
 267   const Double_t khiCorr = 1/TMath::Cos(0.5*khiAng*kDegrad);
 268
 269   Double_t rl,ru;
 270
 271
 272   //
 273   // Get ALICE top node
 274   //
 275
 276   nTop=gAlice->GetGeometry()->GetNode("alice");
 277
 278   //  inner sectors
 279
 280   rl = fTPCParam->GetInnerRadiusLow();
 281   ru = fTPCParam->GetInnerRadiusUp();
 282
 283
 284   for(i=0;i<nLo;i++) {
 285     sprintf(name,"LS%2.2d",i);
 286     name[4]='\0';
 287     sprintf(title,"TPC low sector %3d",i);
 288     title[24]='\0';
 289
 290     tubs = new TTUBS(name,title,"void",rl*kloCorr,ru*kloCorr,250.,
 291                      kloAng*(i-0.5)+kloAngSh,kloAng*(i+0.5)+kloAngSh);
 292     tubs->SetNumberOfDivisions(1);
 293     nTop->cd();
 294     nNode = new TNode(name,title,name,0,0,0,"");
 295     nNode->SetLineColor(kColorTPC);
 296     fNodes->Add(nNode);
 297   }
 298
 299   // Outer sectors
 300
 301   rl = fTPCParam->GetOuterRadiusLow();
 302   ru = fTPCParam->GetOuterRadiusUp();
 303
 304   for(i=0;i<nHi;i++) {
 305     sprintf(name,"US%2.2d",i);
 306     name[4]='\0';
 307     sprintf(title,"TPC upper sector %d",i);
 308     title[24]='\0';
 309     tubs = new TTUBS(name,title,"void",rl*khiCorr,ru*khiCorr,250,
 310                      khiAng*(i-0.5)+khiAngSh,khiAng*(i+0.5)+khiAngSh);
 311     tubs->SetNumberOfDivisions(1);
 312     nTop->cd();
 313     nNode = new TNode(name,title,name,0,0,0,"");
 314     nNode->SetLineColor(kColorTPC);
 315     fNodes->Add(nNode);
 316   }
 317
 318 }
 319
 320 //_____________________________________________________________________________
 321 Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t )
 322 {
 323   //
 324   // Calculate distance from TPC to mouse on the display
 325   // Dummy procedure
 326   //
 327   return 9999;
 328 }
 329
 330 void AliTPC::Clusters2Tracks(TFile *of) {
 331   //-----------------------------------------------------------------
 332   // This is a track finder.
 333   //-----------------------------------------------------------------
 334   AliTPCtracker::Clusters2Tracks(fTPCParam,of);
 335 }
 336
 337 //_____________________________________________________________________________
 338 void AliTPC::CreateMaterials()
 339 {
 340   //-----------------------------------------------
 341   // Create Materials for for TPC simulations
 342   //-----------------------------------------------
 343
 344   //-----------------------------------------------------------------
 345   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
 346   //-----------------------------------------------------------------
 347
 348   Int_t iSXFLD=gAlice->Field()->Integ();
 349   Float_t sXMGMX=gAlice->Field()->Max();
 350
 351   Float_t amat[5]; // atomic numbers
 352   Float_t zmat[5]; // z
 353   Float_t wmat[5]; // proportions
 354
 355   Float_t density;
 356   Float_t apure[2];
 357
 358
 359   //***************** Gases *************************
 360
 361   //-------------------------------------------------
 362   // pure gases
 363   //-------------------------------------------------
 364
 365   // Neon
 366
 367
 368   amat[0]= 20.18;
 369   zmat[0]= 10.;
 370   density = 0.0009;
 371
 372   apure[0]=amat[0];
 373
 374   AliMaterial(20,"Ne",amat[0],zmat[0],density,999.,999.);
 375
 376   // Argon
 377
 378   amat[0]= 39.948;
 379   zmat[0]= 18.;
 380   density = 0.001782;
 381
 382   apure[1]=amat[0];
 383
 384   AliMaterial(21,"Ar",amat[0],zmat[0],density,999.,999.);
 385
 386
 387   //--------------------------------------------------------------
 388   // gases - compounds
 389   //--------------------------------------------------------------
 390
 391   Float_t amol[3];
 392
 393   // CO2
 394
 395   amat[0]=12.011;
 396   amat[1]=15.9994;
 397
 398   zmat[0]=6.;
 399   zmat[1]=8.;
 400
 401   wmat[0]=1.;
 402   wmat[1]=2.;
 403
 404   density=0.001977;
 405
 406   amol[0] = amat[0]*wmat[0]+amat[1]*wmat[1];
 407
 408   AliMixture(10,"CO2",amat,zmat,density,-2,wmat);
 409
 410   // CF4
 411
 412   amat[0]=12.011;
 413   amat[1]=18.998;
 414
 415   zmat[0]=6.;
 416   zmat[1]=9.;
 417
 418   wmat[0]=1.;
 419   wmat[1]=4.;
 420
 421   density=0.003034;
 422
 423   amol[1] = amat[0]*wmat[0]+amat[1]*wmat[1];
 424
 425   AliMixture(11,"CF4",amat,zmat,density,-2,wmat);
 426
 427
 428   // CH4
 429
 430   amat[0]=12.011;
 431   amat[1]=1.;
 432
 433   zmat[0]=6.;
 434   zmat[1]=1.;
 435
 436   wmat[0]=1.;
 437   wmat[1]=4.;
 438
 439   density=0.000717;
 440
 441   amol[2] = amat[0]*wmat[0]+amat[1]*wmat[1];
 442
 443   AliMixture(12,"CH4",amat,zmat,density,-2,wmat);
 444
 445   //----------------------------------------------------------------
 446   // gases - mixtures, ID >= 20 pure gases, <= 10 ID < 20 -compounds
 447   //----------------------------------------------------------------
 448
 449   char namate[21];
 450   density = 0.;
 451   Float_t am=0;
 452   Int_t nc;
 453   Float_t rho,absl,X0,buf[1];
 454   Int_t nbuf;
 455   Float_t a,z;
 456
 457   for(nc = 0;nc<fNoComp;nc++)
 458     {
 459
 460       // retrive material constants
 461
 462       gMC->Gfmate((*fIdmate)[fMixtComp[nc]],namate,a,z,rho,X0,absl,buf,nbuf);
 463
 464       amat[nc] = a;
 465       zmat[nc] = z;
 466
 467       Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10;
 468
 469       am += fMixtProp[nc]*((fMixtComp[nc]>=20) ? apure[nnc] : amol[nnc]);
 470       density += fMixtProp[nc]*rho;  // density of the mixture
 471
 472     }
 473
 474   // mixture proportions by weight!
 475
 476   for(nc = 0;nc<fNoComp;nc++)
 477     {
 478
 479       Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10;
 480
 481       wmat[nc] = fMixtProp[nc]*((fMixtComp[nc]>=20) ?
 482                  apure[nnc] : amol[nnc])/am;
 483
 484     }
 485
 486   // Drift gases 1 - nonsensitive, 2 - sensitive
 487
 488   AliMixture(31,"Drift gas 1",amat,zmat,density,fNoComp,wmat);
 489   AliMixture(32,"Drift gas 2",amat,zmat,density,fNoComp,wmat);
 490
 491
 492   // Air
 493
 494   amat[0] = 14.61;
 495   zmat[0] = 7.3;
 496   density = 0.001205;
 497
 498   AliMaterial(24,"Air",amat[0],zmat[0],density,999.,999.);
 499
 500
 501   //----------------------------------------------------------------------
 502   //               solid materials
 503   //----------------------------------------------------------------------
 504
 505
 506   // Kevlar C14H22O2N2
 507
 508   amat[0] = 12.011;
 509   amat[1] = 1.;
 510   amat[2] = 15.999;
 511   amat[3] = 14.006;
 512
 513   zmat[0] = 6.;
 514   zmat[1] = 1.;
 515   zmat[2] = 8.;
 516   zmat[3] = 7.;
 517
 518   wmat[0] = 14.;
 519   wmat[1] = 22.;
 520   wmat[2] = 2.;
 521   wmat[3] = 2.;
 522
 523   density = 1.45;
 524
 525   AliMixture(34,"Kevlar",amat,zmat,density,-4,wmat);
 526
 527   // NOMEX
 528
 529   amat[0] = 12.011;
 530   amat[1] = 1.;
 531   amat[2] = 15.999;
 532   amat[3] = 14.006;
 533
 534   zmat[0] = 6.;
 535   zmat[1] = 1.;
 536   zmat[2] = 8.;
 537   zmat[3] = 7.;
 538
 539   wmat[0] = 14.;
 540   wmat[1] = 22.;
 541   wmat[2] = 2.;
 542   wmat[3] = 2.;
 543
 544   density = 0.03;
 545
 546
 547   AliMixture(35,"NOMEX",amat,zmat,density,-4,wmat);
 548
 549   // Makrolon C16H18O3
 550
 551   amat[0] = 12.011;
 552   amat[1] = 1.;
 553   amat[2] = 15.999;
 554
 555   zmat[0] = 6.;
 556   zmat[1] = 1.;
 557   zmat[2] = 8.;
 558
 559   wmat[0] = 16.;
 560   wmat[1] = 18.;
 561   wmat[2] = 3.;
 562
 563   density = 1.2;
 564
 565   AliMixture(36,"Makrolon",amat,zmat,density,-3,wmat);
 566
 567   // Mylar C5H4O2
 568
 569   amat[0]=12.011;
 570   amat[1]=1.;
 571   amat[2]=15.9994;
 572
 573   zmat[0]=6.;
 574   zmat[1]=1.;
 575   zmat[2]=8.;
 576
 577   wmat[0]=5.;
 578   wmat[1]=4.;
 579   wmat[2]=2.;
 580
 581   density = 1.39;
 582
 583   AliMixture(37, "Mylar",amat,zmat,density,-3,wmat);
 584
 585   // G10 60% SiO2 + 40% epoxy, I use A and Z for SiO2
 586
 587   amat[0]=28.086;
 588   amat[1]=15.9994;
 589
 590   zmat[0]=14.;
 591   zmat[1]=8.;
 592
 593   wmat[0]=1.;
 594   wmat[1]=2.;
 595
 596   density = 1.7;
 597
 598   AliMixture(38,"SiO2",amat,zmat,2.2,-2,wmat); //SiO2 - quartz
 599
 600   gMC->Gfmate((*fIdmate)[38],namate,amat[0],zmat[0],rho,X0,absl,buf,nbuf);
 601
 602   AliMaterial(39,"G10",amat[0],zmat[0],density,999.,999.);
 603
 604   // Al
 605
 606   amat[0] = 26.98;
 607   zmat[0] = 13.;
 608
 609   density = 2.7;
 610
 611   AliMaterial(40,"Al",amat[0],zmat[0],density,999.,999.);
 612
 613   // Si
 614
 615   amat[0] = 28.086;
 616   zmat[0] = 14.,
 617
 618   density = 2.33;
 619
 620   AliMaterial(41,"Si",amat[0],zmat[0],density,999.,999.);
 621
 622   // Cu
 623
 624   amat[0] = 63.546;
 625   zmat[0] = 29.;
 626
 627   density = 8.96;
 628
 629   AliMaterial(42,"Cu",amat[0],zmat[0],density,999.,999.);
 630
 631   // Tedlar C2H3F
 632
 633   amat[0] = 12.011;
 634   amat[1] = 1.;
 635   amat[2] = 18.998;
 636
 637   zmat[0] = 6.;
 638   zmat[1] = 1.;
 639   zmat[2] = 9.;
 640
 641   wmat[0] = 2.;
 642   wmat[1] = 3.;
 643   wmat[2] = 1.;
 644
 645   density = 1.71;
 646
 647   AliMixture(43, "Tedlar",amat,zmat,density,-3,wmat);
 648
 649
 650   // Plexiglas  C5H8O2
 651
 652   amat[0]=12.011;
 653   amat[1]=1.;
 654   amat[2]=15.9994;
 655
 656   zmat[0]=6.;
 657   zmat[1]=1.;
 658   zmat[2]=8.;
 659
 660   wmat[0]=5.;
 661   wmat[1]=8.;
 662   wmat[2]=2.;
 663
 664   density=1.18;
 665
 666   AliMixture(44,"Plexiglas",amat,zmat,density,-3,wmat);
 667
 668
 669
 670   //----------------------------------------------------------
 671   // tracking media for gases
 672   //----------------------------------------------------------
 673
 674   AliMedium(0, "Air", 24, 0, iSXFLD, sXMGMX, 10., 999., .1, .01, .1);
 675   AliMedium(1, "Drift gas 1", 31, 0, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
 676   AliMedium(2, "Drift gas 2", 32, 1, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
 677   AliMedium(3,"CO2",10,0, iSXFLD, sXMGMX, 10., 999.,.1, .001, .001);
 678
 679   //-----------------------------------------------------------
 680   // tracking media for solids
 681   //-----------------------------------------------------------
 682
 683   AliMedium(4,"Al",40,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 684   AliMedium(5,"Kevlar",34,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 685   AliMedium(6,"Nomex",35,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 686   AliMedium(7,"Makrolon",36,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 687   AliMedium(8,"Mylar",37,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 688   AliMedium(9,"Tedlar",43,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 689   AliMedium(10,"Cu",42,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 690   AliMedium(11,"Si",41,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 691   AliMedium(12,"G10",39,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 692   AliMedium(13,"Plexiglas",44,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 693
 694 }
 695
 696
 697 void AliTPC::Digits2Clusters(TFile *of)
 698 {
 699   //-----------------------------------------------------------------
 700   // This is a simple cluster finder.
 701   //-----------------------------------------------------------------
 702   AliTPCclusterer::Digits2Clusters(fTPCParam,of);
 703 }
 704
 705 extern Double_t SigmaY2(Double_t, Double_t, Double_t);
 706 extern Double_t SigmaZ2(Double_t, Double_t);
 707 //_____________________________________________________________________________
 708 void AliTPC::Hits2Clusters(TFile *of)
 709 {
 710   //--------------------------------------------------------
 711   // TPC simple cluster generator from hits
 712   // obtained from the TPC Fast Simulator
 713   // The point errors are taken from the parametrization
 714   //--------------------------------------------------------
 715
 716   //-----------------------------------------------------------------
 717   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
 718   //-----------------------------------------------------------------
 719   // Adopted to Marian's cluster data structure by I.Belikov, CERN,
 720   // Jouri.Belikov@cern.ch
 721   //----------------------------------------------------------------
 722
 723   /////////////////////////////////////////////////////////////////////////////
 724   //
 725   //---------------------------------------------------------------------
 726   //   ALICE TPC Cluster Parameters
 727   //--------------------------------------------------------------------
 728
 729
 730
 731   // Cluster width in rphi
 732   const Float_t kACrphi=0.18322;
 733   const Float_t kBCrphi=0.59551e-3;
 734   const Float_t kCCrphi=0.60952e-1;
 735   // Cluster width in z
 736   const Float_t kACz=0.19081;
 737   const Float_t kBCz=0.55938e-3;
 738   const Float_t kCCz=0.30428;
 739
 740   TDirectory *savedir=gDirectory;
 741
 742   if (!of->IsOpen()) {
 743      cerr<<"AliTPC::Hits2Clusters(): output file not open !\n";
 744      return;
 745   }
 746
 747    if(fTPCParam == 0){
 748      printf("AliTPCParam MUST be created firstly\n");
 749      return;
 750    }
 751
 752   Float_t sigmaRphi,sigmaZ,clRphi,clZ;
 753   //
 754   TParticle *particle; // pointer to a given particle
 755   AliTPChit *tpcHit; // pointer to a sigle TPC hit
 756   TClonesArray *particles; //pointer to the particle list
 757   Int_t sector;
 758   Int_t ipart;
 759   Float_t xyz[5];
 760   Float_t pl,pt,tanth,rpad,ratio;
 761   Float_t cph,sph;
 762
 763   //---------------------------------------------------------------
 764   //  Get the access to the tracks
 765   //---------------------------------------------------------------
 766
 767   TTree *tH = gAlice->TreeH();
 768   Stat_t ntracks = tH->GetEntries();
 769   particles=gAlice->Particles();
 770
 771   //Switch to the output file
 772   of->cd();
 773
 774   fTPCParam->Write(fTPCParam->GetTitle());
 775   AliTPCClustersArray carray;
 776   carray.Setup(fTPCParam);
 777   carray.SetClusterType("AliTPCcluster");
 778   carray.MakeTree();
 779
 780   Int_t nclusters=0; //cluster counter
 781
 782   //------------------------------------------------------------
 783   // Loop over all sectors (72 sectors for 20 deg
 784   // segmentation for both lower and upper sectors)
 785   // Sectors 0-35 are lower sectors, 0-17 z>0, 17-35 z<0
 786   // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
 787   //
 788   // First cluster for sector 0 starts at "0"
 789   //------------------------------------------------------------
 790
 791   for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++){
 792     //MI change
 793     fTPCParam->AdjustCosSin(isec,cph,sph);
 794
 795     //------------------------------------------------------------
 796     // Loop over tracks
 797     //------------------------------------------------------------
 798
 799     for(Int_t track=0;track<ntracks;track++){
 800       ResetHits();
 801       tH->GetEvent(track);
 802       //
 803       //  Get number of the TPC hits
 804       //
 805       // nhits=fHits->GetEntriesFast();
 806       //
 807
 808        tpcHit = (AliTPChit*)FirstHit(-1);
 809
 810       // Loop over hits
 811       //
 812        //   for(Int_t hit=0;hit<nhits;hit++){
 813        //tpcHit=(AliTPChit*)fHits->UncheckedAt(hit);
 814
 815        while(tpcHit){
 816
 817          if (tpcHit->fQ == 0.) {
 818            tpcHit = (AliTPChit*) NextHit();
 819            continue; //information about track (I.Belikov)
 820          }
 821         sector=tpcHit->fSector; // sector number
 822
 823
 824         //      if(sector != isec) continue; //terminate iteration
 825
 826        if(sector != isec){
 827          tpcHit = (AliTPChit*) NextHit();
 828          continue;
 829        }
 830         ipart=tpcHit->Track();
 831         particle=(TParticle*)particles->UncheckedAt(ipart);
 832         pl=particle->Pz();
 833         pt=particle->Pt();
 834         if(pt < 1.e-9) pt=1.e-9;
 835         tanth=pl/pt;
 836         tanth = TMath::Abs(tanth);
 837         rpad=TMath::Sqrt(tpcHit->X()*tpcHit->X() + tpcHit->Y()*tpcHit->Y());
 838         ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
 839
 840         //   space-point resolutions
 841
 842         sigmaRphi=SigmaY2(rpad,tanth,pt);
 843         sigmaZ   =SigmaZ2(rpad,tanth   );
 844
 845         //   cluster widths
 846
 847         clRphi=kACrphi-kBCrphi*rpad*tanth+kCCrphi*ratio*ratio;
 848         clZ=kACz-kBCz*rpad*tanth+kCCz*tanth*tanth;
 849
 850         // temporary protection
 851
 852         if(sigmaRphi < 0.) sigmaRphi=0.4e-3;
 853         if(sigmaZ < 0.) sigmaZ=0.4e-3;
 854         if(clRphi < 0.) clRphi=2.5e-3;
 855         if(clZ < 0.) clZ=2.5e-5;
 856
 857         //
 858
 859         //
 860         // smearing --> rotate to the 1 (13) or to the 25 (49) sector,
 861         // then the inaccuracy in a X-Y plane is only along Y (pad row)!
 862         //
 863         Float_t xprim= tpcHit->X()*cph + tpcHit->Y()*sph;
 864         Float_t yprim=-tpcHit->X()*sph + tpcHit->Y()*cph;
 865         xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigmaRphi));   // y
 866           Float_t alpha=(isec < fTPCParam->GetNInnerSector()) ?
 867           fTPCParam->GetInnerAngle() : fTPCParam->GetOuterAngle();
 868           Float_t ymax=xprim*TMath::Tan(0.5*alpha);
 869           if (TMath::Abs(xyz[0])>ymax) xyz[0]=yprim;
 870         xyz[1]=gRandom->Gaus(tpcHit->Z(),TMath::Sqrt(sigmaZ)); // z
 871           if (TMath::Abs(xyz[1])>fTPCParam->GetZLength()) xyz[1]=tpcHit->Z();
 872         xyz[2]=sigmaRphi;                                     // fSigmaY2
 873         xyz[3]=sigmaZ;                                        // fSigmaZ2
 874         xyz[4]=tpcHit->fQ;                                    // q
 875
 876         AliTPCClustersRow *clrow=carray.GetRow(sector,tpcHit->fPadRow);
 877         if (!clrow) clrow=carray.CreateRow(sector,tpcHit->fPadRow);
 878
 879         Int_t tracks[3]={tpcHit->Track(), -1, -1};
 880         AliTPCcluster cluster(tracks,xyz);
 881
 882         clrow->InsertCluster(&cluster); nclusters++;
 883
 884         tpcHit = (AliTPChit*)NextHit();
 885
 886
 887       } // end of loop over hits
 888
 889     }   // end of loop over tracks
 890
 891     Int_t nrows=fTPCParam->GetNRow(isec);
 892     for (Int_t irow=0; irow<nrows; irow++) {
 893         AliTPCClustersRow *clrow=carray.GetRow(isec,irow);
 894         if (!clrow) continue;
 895         carray.StoreRow(isec,irow);
 896         carray.ClearRow(isec,irow);
 897     }
 898
 899   } // end of loop over sectors
 900
 901   cerr<<"Number of made clusters : "<<nclusters<<"                        \n";
 902
 903   carray.GetTree()->Write();
 904
 905   savedir->cd(); //switch back to the input file
 906
 907 } // end of function
 908
 909 //_________________________________________________________________
 910 void AliTPC::Hits2ExactClustersSector(Int_t isec)
 911 {
 912   //--------------------------------------------------------
 913   //calculate exact cross point of track and given pad row
 914   //resulting values are expressed in "digit" coordinata
 915   //--------------------------------------------------------
 916
 917   //-----------------------------------------------------------------
 918   // Origin: Marian Ivanov  GSI Darmstadt, m.ivanov@gsi.de
 919   //-----------------------------------------------------------------
 920   //
 921   if (fClustersArray==0){
 922     return;
 923   }
 924   //
 925   TParticle *particle; // pointer to a given particle
 926   AliTPChit *tpcHit; // pointer to a sigle TPC hit
 927   TClonesArray *particles; //pointer to the particle list
 928   Int_t sector,nhits;
 929   Int_t ipart;
 930   const Int_t kcmaxhits=30000;
 931   TVector * xxxx = new TVector(kcmaxhits*4);
 932   TVector & xxx = *xxxx;
 933   Int_t maxhits = kcmaxhits;
 934   //construct array for each padrow
 935   for (Int_t i=0; i<fTPCParam->GetNRow(isec);i++)
 936     fClustersArray->CreateRow(isec,i);
 937
 938   //---------------------------------------------------------------
 939   //  Get the access to the tracks
 940   //---------------------------------------------------------------
 941
 942   TTree *tH = gAlice->TreeH();
 943   Stat_t ntracks = tH->GetEntries();
 944   particles=gAlice->Particles();
 945   Int_t npart = particles->GetEntriesFast();
 946
 947   //------------------------------------------------------------
 948   // Loop over tracks
 949   //------------------------------------------------------------
 950
 951   for(Int_t track=0;track<ntracks;track++){
 952     ResetHits();
 953     tH->GetEvent(track);
 954     //
 955     //  Get number of the TPC hits and a pointer
 956     //  to the particles
 957     //
 958     nhits=fHits->GetEntriesFast();
 959     //
 960     // Loop over hits
 961     //
 962     Int_t currentIndex=0;
 963     Int_t lastrow=-1;  //last writen row
 964     for(Int_t hit=0;hit<nhits;hit++){
 965       tpcHit=(AliTPChit*)fHits->UncheckedAt(hit);
 966       if (tpcHit==0) continue;
 967       sector=tpcHit->fSector; // sector number
 968       if(sector != isec) continue;
 969       ipart=tpcHit->Track();
 970       if (ipart<npart) particle=(TParticle*)particles->UncheckedAt(ipart);
 971
 972       //find row number
 973
 974       Float_t  x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
 975       Int_t    index[3]={1,isec,0};
 976       Int_t    currentrow = fTPCParam->GetPadRow(x,index) ;
 977       if (currentrow<0) continue;
 978       if (lastrow<0) lastrow=currentrow;
 979       if (currentrow==lastrow){
 980         if ( currentIndex>=maxhits){
 981           maxhits+=kcmaxhits;
 982           xxx.ResizeTo(4*maxhits);
 983         }
 984         xxx(currentIndex*4)=x[0];
 985         xxx(currentIndex*4+1)=x[1];
 986         xxx(currentIndex*4+2)=x[2];
 987         xxx(currentIndex*4+3)=tpcHit->fQ;
 988         currentIndex++;
 989       }
 990       else
 991         if (currentIndex>2){
 992           Float_t sumx=0;
 993           Float_t sumx2=0;
 994           Float_t sumx3=0;
 995           Float_t sumx4=0;
 996           Float_t sumy=0;
 997           Float_t sumxy=0;
 998           Float_t sumx2y=0;
 999           Float_t sumz=0;
1000           Float_t sumxz=0;
1001           Float_t sumx2z=0;
1002           Float_t sumq=0;
1003           for (Int_t index=0;index<currentIndex;index++){
1004             Float_t x,x2,x3,x4;
1005             x=x2=x3=x4=xxx(index*4);
1006             x2*=x;
1007             x3*=x2;
1008             x4*=x3;
1009             sumx+=x;
1010             sumx2+=x2;
1011             sumx3+=x3;
1012             sumx4+=x4;
1013             sumy+=xxx(index*4+1);
1014             sumxy+=xxx(index*4+1)*x;
1015             sumx2y+=xxx(index*4+1)*x2;
1016             sumz+=xxx(index*4+2);
1017             sumxz+=xxx(index*4+2)*x;
1018             sumx2z+=xxx(index*4+2)*x2;
1019             sumq+=xxx(index*4+3);
1020           }
1021           Float_t centralPad = (fTPCParam->GetNPads(isec,lastrow)-1)/2;
1022           Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
1023             sumx2*(sumx*sumx3-sumx2*sumx2);
1024
1025           Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
1026             sumx2*(sumxy*sumx3-sumx2y*sumx2);
1027           Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
1028             sumx2*(sumxz*sumx3-sumx2z*sumx2);
1029
1030           Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
1031             sumx2*(sumx*sumx2y-sumx2*sumxy);
1032           Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
1033             sumx2*(sumx*sumx2z-sumx2*sumxz);
1034
1035           Float_t y=detay/det+centralPad;
1036           Float_t z=detaz/det;
1037           Float_t by=detby/det; //y angle
1038           Float_t bz=detbz/det; //z angle
1039           sumy/=Float_t(currentIndex);
1040           sumz/=Float_t(currentIndex);
1041           AliComplexCluster cl;
1042           cl.fX=z;
1043           cl.fY=y;
1044           cl.fQ=sumq;
1045           cl.fSigmaX2=bz;
1046           cl.fSigmaY2=by;
1047           cl.fTracks[0]=ipart;
1048
1049           AliTPCClustersRow * row = (fClustersArray->GetRow(isec,lastrow));
1050           if (row!=0) row->InsertCluster(&cl);
1051           currentIndex=0;
1052           lastrow=currentrow;
1053         } //end of calculating cluster for given row
1054
1055
1056
1057     } // end of loop over hits
1058   }   // end of loop over tracks
1059   //write padrows to tree
1060   for (Int_t ii=0; ii<fTPCParam->GetNRow(isec);ii++) {
1061     fClustersArray->StoreRow(isec,ii);
1062     fClustersArray->ClearRow(isec,ii);
1063   }
1064   xxxx->Delete();
1065
1066 }
1067
1068 //__________________________________________________________________
1069 void AliTPC::Hits2Digits()
1070 {
1071  //----------------------------------------------------
1072  // Loop over all sectors
1073  //----------------------------------------------------
1074
1075   if(fTPCParam == 0){
1076     printf("AliTPCParam MUST be created firstly\n");
1077     return;
1078   }
1079
1080  for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++) Hits2DigitsSector(isec);
1081
1082 }
1083
1084
1085 //_____________________________________________________________________________
1086 void AliTPC::Hits2DigitsSector(Int_t isec)
1087 {
1088   //-------------------------------------------------------------------
1089   // TPC conversion from hits to digits.
1090   //-------------------------------------------------------------------
1091
1092   //-----------------------------------------------------------------
1093   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1094   //-----------------------------------------------------------------
1095
1096   //-------------------------------------------------------
1097   //  Get the access to the track hits
1098   //-------------------------------------------------------
1099
1100
1101   TTree *tH = gAlice->TreeH(); // pointer to the hits tree
1102   Stat_t ntracks = tH->GetEntries();
1103
1104   if( ntracks > 0){
1105
1106   //-------------------------------------------
1107   //  Only if there are any tracks...
1108   //-------------------------------------------
1109
1110     TObjArray **row;
1111
1112       printf("*** Processing sector number %d ***\n",isec);
1113
1114       Int_t nrows =fTPCParam->GetNRow(isec);
1115
1116       row= new TObjArray* [nrows];
1117
1118       MakeSector(isec,nrows,tH,ntracks,row);
1119
1120       //--------------------------------------------------------
1121       //   Digitize this sector, row by row
1122       //   row[i] is the pointer to the TObjArray of TVectors,
1123       //   each one containing electrons accepted on this
1124       //   row, assigned into tracks
1125       //--------------------------------------------------------
1126
1127       Int_t i;
1128
1129       if (fDigitsArray->GetTree()==0) fDigitsArray->MakeTree();
1130
1131       for (i=0;i<nrows;i++){
1132
1133         AliDigits * dig = fDigitsArray->CreateRow(isec,i);
1134
1135         DigitizeRow(i,isec,row);
1136
1137         fDigitsArray->StoreRow(isec,i);
1138
1139         Int_t ndig = dig->GetDigitSize();
1140
1141         printf("*** Sector, row, compressed digits %d %d %d ***\n",isec,i,ndig);
1142
1143         fDigitsArray->ClearRow(isec,i);
1144
1145
1146        } // end of the sector digitization
1147
1148       for(i=0;i<nrows;i++){
1149         row[i]->Delete();
1150         delete row[i];
1151       }
1152
1153        delete [] row; // delete the array of pointers to TObjArray-s
1154
1155   } // ntracks >0
1156
1157 } // end of Hits2DigitsSector
1158
1159
1160 //_____________________________________________________________________________
1161 void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows)
1162 {
1163   //-----------------------------------------------------------
1164   // Single row digitization, coupling from the neighbouring
1165   // rows taken into account
1166   //-----------------------------------------------------------
1167
1168   //-----------------------------------------------------------------
1169   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1170   // Modified: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
1171   //-----------------------------------------------------------------
1172
1173
1174   Float_t zerosup = fTPCParam->GetZeroSup();
1175   Int_t nrows =fTPCParam->GetNRow(isec);
1176   fCurrentIndex[1]= isec;
1177
1178
1179   Int_t nofPads = fTPCParam->GetNPads(isec,irow);
1180   Int_t nofTbins = fTPCParam->GetMaxTBin();
1181   Int_t indexRange[4];
1182   //
1183   //  Integrated signal for this row
1184   //  and a single track signal
1185   //
1186   TMatrix *m1   = new TMatrix(0,nofPads,0,nofTbins); // integrated
1187   TMatrix *m2   = new TMatrix(0,nofPads,0,nofTbins); // single
1188   //
1189   TMatrix &total  = *m1;
1190
1191   //  Array of pointers to the label-signal list
1192
1193   Int_t nofDigits = nofPads*nofTbins; // number of digits for this row
1194   Float_t  **pList = new Float_t* [nofDigits];
1195
1196   Int_t lp;
1197   Int_t i1;
1198   for(lp=0;lp<nofDigits;lp++)pList[lp]=0; // set all pointers to NULL
1199   //
1200   //calculate signal
1201   //
1202   Int_t row1 = TMath::Max(irow-fTPCParam->GetNCrossRows(),0);
1203   Int_t row2 = TMath::Min(irow+fTPCParam->GetNCrossRows(),nrows-1);
1204   for (Int_t row= row1;row<=row2;row++){
1205     Int_t nTracks= rows[row]->GetEntries();
1206     for (i1=0;i1<nTracks;i1++){
1207       fCurrentIndex[2]= row;
1208       fCurrentIndex[3]=irow;
1209       if (row==irow){
1210         m2->Zero();  // clear single track signal matrix
1211         Float_t trackLabel = GetSignal(rows[row],i1,m2,m1,indexRange);
1212         GetList(trackLabel,nofPads,m2,indexRange,pList);
1213       }
1214       else   GetSignal(rows[row],i1,0,m1,indexRange);
1215     }
1216   }
1217
1218   Int_t tracks[3];
1219
1220   AliDigits *dig = fDigitsArray->GetRow(isec,irow);
1221   for(Int_t ip=0;ip<nofPads;ip++){
1222     for(Int_t it=0;it<nofTbins;it++){
1223
1224       Float_t q = total(ip,it);
1225
1226       Int_t gi =it*nofPads+ip; // global index
1227
1228       q = gRandom->Gaus(q,fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac());
1229
1230       q = (Int_t)q;
1231
1232       if(q <=zerosup) continue; // do not fill zeros
1233       if(q > fTPCParam->GetADCSat()) q = fTPCParam->GetADCSat();  // saturation
1234
1235       //
1236       //  "real" signal or electronic noise (list = -1)?
1237       //
1238
1239       for(Int_t j1=0;j1<3;j1++){
1240         tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -1;
1241       }
1242
1243 //Begin_Html
1244 /*
1245   <A NAME="AliDigits"></A>
1246   using of AliDigits object
1247 */
1248 //End_Html
1249       dig->SetDigitFast((Short_t)q,it,ip);
1250       if (fDigitsArray->IsSimulated())
1251         {
1252          ((AliSimDigits*)dig)->SetTrackIDFast(tracks[0],it,ip,0);
1253          ((AliSimDigits*)dig)->SetTrackIDFast(tracks[1],it,ip,1);
1254          ((AliSimDigits*)dig)->SetTrackIDFast(tracks[2],it,ip,2);
1255         }
1256
1257
1258     } // end of loop over time buckets
1259   }  // end of lop over pads
1260
1261   //
1262   //  This row has been digitized, delete nonused stuff
1263   //
1264
1265   for(lp=0;lp<nofDigits;lp++){
1266     if(pList[lp]) delete [] pList[lp];
1267   }
1268
1269   delete [] pList;
1270
1271   delete m1;
1272   delete m2;
1273   //  delete m3;
1274
1275 } // end of DigitizeRow
1276
1277 //_____________________________________________________________________________
1278
1279 Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr, TMatrix *m1, TMatrix *m2,
1280                           Int_t *indexRange)
1281 {
1282
1283   //---------------------------------------------------------------
1284   //  Calculates 2-D signal (pad,time) for a single track,
1285   //  returns a pointer to the signal matrix and the track label
1286   //  No digitization is performed at this level!!!
1287   //---------------------------------------------------------------
1288
1289   //-----------------------------------------------------------------
1290   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1291   // Modified: Marian Ivanov
1292   //-----------------------------------------------------------------
1293
1294   TVector *tv;
1295
1296   tv = (TVector*)p1->At(ntr); // pointer to a track
1297   TVector &v = *tv;
1298
1299   Float_t label = v(0);
1300   Int_t centralPad = (fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3])-1)/2;
1301
1302   Int_t nElectrons = (tv->GetNrows()-1)/4;
1303   indexRange[0]=9999; // min pad
1304   indexRange[1]=-1; // max pad
1305   indexRange[2]=9999; //min time
1306   indexRange[3]=-1; // max time
1307
1308   //  Float_t IneffFactor = 0.5; // inefficiency in the gain close to the edge, as above
1309
1310   TMatrix &signal = *m1;
1311   TMatrix &total = *m2;
1312   //
1313   //  Loop over all electrons
1314   //
1315   for(Int_t nel=0; nel<nElectrons; nel++){
1316     Int_t idx=nel*4;
1317     Float_t aval =  v(idx+4);
1318     Float_t eltoadcfac=aval*fTPCParam->GetTotalNormFac();
1319     Float_t xyz[3]={v(idx+1),v(idx+2),v(idx+3)};
1320     Int_t n = fTPCParam->CalcResponse(xyz,fCurrentIndex,fCurrentIndex[3]);
1321
1322     if (n>0) for (Int_t i =0; i<n; i++){
1323        Int_t *index = fTPCParam->GetResBin(i);
1324        Int_t pad=index[1]+centralPad;  //in digit coordinates central pad has coordinate 0
1325        if ( ( pad<(fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3]))) && (pad>0)) {
1326          Int_t time=index[2];
1327          Float_t weight = fTPCParam->GetResWeight(i); //we normalise response to ADC channel
1328          weight *= eltoadcfac;
1329
1330          if (m1!=0) signal(pad,time)+=weight;
1331          total(pad,time)+=weight;
1332          indexRange[0]=TMath::Min(indexRange[0],pad);
1333          indexRange[1]=TMath::Max(indexRange[1],pad);
1334          indexRange[2]=TMath::Min(indexRange[2],time);
1335          indexRange[3]=TMath::Max(indexRange[3],time);
1336        }
1337     }
1338   } // end of loop over electrons
1339
1340   return label; // returns track label when finished
1341 }
1342
1343 //_____________________________________________________________________________
1344 void AliTPC::GetList(Float_t label,Int_t np,TMatrix *m,Int_t *indexRange,
1345                      Float_t **pList)
1346 {
1347   //----------------------------------------------------------------------
1348   //  Updates the list of tracks contributing to digits for a given row
1349   //----------------------------------------------------------------------
1350
1351   //-----------------------------------------------------------------
1352   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1353   //-----------------------------------------------------------------
1354
1355   TMatrix &signal = *m;
1356
1357   // lop over nonzero digits
1358
1359   for(Int_t it=indexRange[2];it<indexRange[3]+1;it++){
1360     for(Int_t ip=indexRange[0];ip<indexRange[1]+1;ip++){
1361
1362
1363         // accept only the contribution larger than 500 electrons (1/2 s_noise)
1364
1365         if(signal(ip,it)<0.5) continue;
1366
1367
1368         Int_t globalIndex = it*np+ip; // globalIndex starts from 0!
1369
1370         if(!pList[globalIndex]){
1371
1372           //
1373           // Create new list (6 elements - 3 signals and 3 labels),
1374           //
1375
1376           pList[globalIndex] = new Float_t [6];
1377
1378           // set list to -1
1379
1380           *pList[globalIndex] = -1.;
1381           *(pList[globalIndex]+1) = -1.;
1382           *(pList[globalIndex]+2) = -1.;
1383           *(pList[globalIndex]+3) = -1.;
1384           *(pList[globalIndex]+4) = -1.;
1385           *(pList[globalIndex]+5) = -1.;
1386
1387
1388           *pList[globalIndex] = label;
1389           *(pList[globalIndex]+3) = signal(ip,it);
1390         }
1391         else{
1392
1393           // check the signal magnitude
1394
1395           Float_t highest = *(pList[globalIndex]+3);
1396           Float_t middle = *(pList[globalIndex]+4);
1397           Float_t lowest = *(pList[globalIndex]+5);
1398
1399           //
1400           //  compare the new signal with already existing list
1401           //
1402
1403           if(signal(ip,it)<lowest) continue; // neglect this track
1404
1405           //
1406
1407           if (signal(ip,it)>highest){
1408             *(pList[globalIndex]+5) = middle;
1409             *(pList[globalIndex]+4) = highest;
1410             *(pList[globalIndex]+3) = signal(ip,it);
1411
1412             *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
1413             *(pList[globalIndex]+1) = *pList[globalIndex];
1414             *pList[globalIndex] = label;
1415           }
1416           else if (signal(ip,it)>middle){
1417             *(pList[globalIndex]+5) = middle;
1418             *(pList[globalIndex]+4) = signal(ip,it);
1419
1420             *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
1421             *(pList[globalIndex]+1) = label;
1422           }
1423           else{
1424             *(pList[globalIndex]+5) = signal(ip,it);
1425             *(pList[globalIndex]+2) = label;
1426           }
1427         }
1428
1429     } // end of loop over pads
1430   } // end of loop over time bins
1431
1432
1433
1434 }//end of GetList
1435 //___________________________________________________________________
1436 void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
1437                         Stat_t ntracks,TObjArray **row)
1438 {
1439
1440   //-----------------------------------------------------------------
1441   // Prepares the sector digitization, creates the vectors of
1442   // tracks for each row of this sector. The track vector
1443   // contains the track label and the position of electrons.
1444   //-----------------------------------------------------------------
1445
1446   //-----------------------------------------------------------------
1447   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1448   //-----------------------------------------------------------------
1449
1450   Float_t gasgain = fTPCParam->GetGasGain();
1451   Int_t i;
1452   Float_t xyz[4];
1453
1454   AliTPChit *tpcHit; // pointer to a sigle TPC hit
1455   //MI change
1456   TBranch * branch=0;
1457   if (fHitType&2) branch = TH->GetBranch("TPC2");
1458   else branch = TH->GetBranch("TPC");
1459
1460
1461   //----------------------------------------------
1462   // Create TObjArray-s, one for each row,
1463   // each TObjArray will store the TVectors
1464   // of electrons, one TVector per each track.
1465   //----------------------------------------------
1466
1467   Int_t *nofElectrons = new Int_t [nrows]; // electron counter for each row
1468   TVector **tracks = new TVector* [nrows]; //pointers to the track vectors
1469   for(i=0; i<nrows; i++){
1470     row[i] = new TObjArray;
1471     nofElectrons[i]=0;
1472     tracks[i]=0;
1473   }
1474
1475
1476
1477   //--------------------------------------------------------------------
1478   //  Loop over tracks, the "track" contains the full history
1479   //--------------------------------------------------------------------
1480
1481   Int_t previousTrack,currentTrack;
1482   previousTrack = -1; // nothing to store so far!
1483
1484   for(Int_t track=0;track<ntracks;track++){
1485     Bool_t isInSector=kTRUE;
1486     ResetHits();
1487
1488     if (fHitType&2) {
1489       isInSector=kFALSE;
1490       TBranch * br = TH->GetBranch("fTrackHitsInfo");
1491       br->GetEvent(track);
1492       AliObjectArray * ar = fTrackHits->fTrackHitsInfo;
1493       for (UInt_t j=0;j<ar->GetSize();j++){
1494         if (  ((AliTrackHitsInfo*)ar->At(j))->fVolumeID==isec) isInSector=kTRUE;
1495       }
1496     }
1497     if (!isInSector) continue;
1498     //MI change
1499     branch->GetEntry(track); // get next track
1500
1501     //M.I. changes
1502
1503     tpcHit = (AliTPChit*)FirstHit(-1);
1504
1505     //--------------------------------------------------------------
1506     //  Loop over hits
1507     //--------------------------------------------------------------
1508
1509
1510     while(tpcHit){
1511
1512       Int_t sector=tpcHit->fSector; // sector number
1513       //      if(sector != isec) continue;
1514       if(sector != isec){
1515         tpcHit = (AliTPChit*) NextHit();
1516         continue;
1517       }
1518
1519         currentTrack = tpcHit->Track(); // track number
1520
1521
1522         if(currentTrack != previousTrack){
1523
1524            // store already filled fTrack
1525
1526            for(i=0;i<nrows;i++){
1527              if(previousTrack != -1){
1528                if(nofElectrons[i]>0){
1529                  TVector &v = *tracks[i];
1530                  v(0) = previousTrack;
1531                  tracks[i]->ResizeTo(4*nofElectrons[i]+1); // shrink if necessary
1532                  row[i]->Add(tracks[i]);
1533                }
1534                else{
1535                  delete tracks[i]; // delete empty TVector
1536                  tracks[i]=0;
1537                }
1538              }
1539
1540              nofElectrons[i]=0;
1541              tracks[i] = new TVector(481); // TVectors for the next fTrack
1542
1543            } // end of loop over rows
1544
1545            previousTrack=currentTrack; // update track label
1546         }
1547
1548         Int_t qI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons)
1549
1550        //---------------------------------------------------
1551        //  Calculate the electron attachment probability
1552        //---------------------------------------------------
1553
1554
1555         Float_t time = 1.e6*(fTPCParam->GetZLength()-TMath::Abs(tpcHit->Z()))
1556                                                         /fTPCParam->GetDriftV();
1557         // in microseconds!
1558         Float_t attProb = fTPCParam->GetAttCoef()*
1559           fTPCParam->GetOxyCont()*time; //  fraction!
1560
1561         //-----------------------------------------------
1562         //  Loop over electrons
1563         //-----------------------------------------------
1564         Int_t index[3];
1565         index[1]=isec;
1566         for(Int_t nel=0;nel<qI;nel++){
1567           // skip if electron lost due to the attachment
1568           if((gRandom->Rndm(0)) < attProb) continue; // electron lost!
1569           xyz[0]=tpcHit->X();
1570           xyz[1]=tpcHit->Y();
1571           xyz[2]=tpcHit->Z();
1572           xyz[3]= (Float_t) (-gasgain*TMath::Log(gRandom->Rndm()));
1573           index[0]=1;
1574
1575           TransportElectron(xyz,index); //MI change -august
1576           Int_t rowNumber;
1577           fTPCParam->GetPadRow(xyz,index); //MI change august
1578           rowNumber = index[2];
1579           //transform position to local digit coordinates
1580           //relative to nearest pad row
1581           if ((rowNumber<0)||rowNumber>=fTPCParam->GetNRow(isec)) continue;
1582           nofElectrons[rowNumber]++;
1583           //----------------------------------
1584           // Expand vector if necessary
1585           //----------------------------------
1586           if(nofElectrons[rowNumber]>120){
1587             Int_t range = tracks[rowNumber]->GetNrows();
1588             if((nofElectrons[rowNumber])>(range-1)/4){
1589
1590               tracks[rowNumber]->ResizeTo(range+400); // Add 100 electrons
1591             }
1592           }
1593
1594           TVector &v = *tracks[rowNumber];
1595           Int_t idx = 4*nofElectrons[rowNumber]-3;
1596
1597           v(idx)=  xyz[0];   // X - pad row coordinate
1598           v(idx+1)=xyz[1];   // Y - pad coordinate (along the pad-row)
1599           v(idx+2)=xyz[2];   // Z - time bin coordinate
1600           v(idx+3)=xyz[3];   // avalanche size
1601         } // end of loop over electrons
1602
1603         tpcHit = (AliTPChit*)NextHit();
1604
1605       } // end of loop over hits
1606     } // end of loop over tracks
1607
1608     //
1609     //   store remaining track (the last one) if not empty
1610     //
1611
1612      for(i=0;i<nrows;i++){
1613        if(nofElectrons[i]>0){
1614           TVector &v = *tracks[i];
1615           v(0) = previousTrack;
1616           tracks[i]->ResizeTo(4*nofElectrons[i]+1); // shrink if necessary
1617           row[i]->Add(tracks[i]);
1618         }
1619         else{
1620           delete tracks[i];
1621           tracks[i]=0;
1622         }
1623       }
1624
1625           delete [] tracks;
1626           delete [] nofElectrons;
1627
1628
1629 } // end of MakeSector
1630
1631
1632 //_____________________________________________________________________________
1633 void AliTPC::Init()
1634 {
1635   //
1636   // Initialise TPC detector after definition of geometry
1637   //
1638   Int_t i;
1639   //
1640   printf("\n");
1641   for(i=0;i<35;i++) printf("*");
1642   printf(" TPC_INIT ");
1643   for(i=0;i<35;i++) printf("*");
1644   printf("\n");
1645   //
1646   for(i=0;i<80;i++) printf("*");
1647   printf("\n");
1648 }
1649
1650 //_____________________________________________________________________________
1651 void AliTPC::MakeBranch(Option_t* option)
1652 {
1653   //
1654   // Create Tree branches for the TPC.
1655   //
1656   Int_t buffersize = 4000;
1657   char branchname[10];
1658   sprintf(branchname,"%s",GetName());
1659
1660   AliDetector::MakeBranch(option);
1661
1662   char *d = strstr(option,"D");
1663
1664   if (fDigits   && gAlice->TreeD() && d) {
1665     gAlice->TreeD()->Branch(branchname,&fDigits, buffersize);
1666     printf("Making Branch %s for digits\n",branchname);
1667   }
1668
1669   if (fHitType&2) MakeBranch2(option); // MI change 14.09.2000
1670 }
1671
1672 //_____________________________________________________________________________
1673 void AliTPC::ResetDigits()
1674 {
1675   //
1676   // Reset number of digits and the digits array for this detector
1677   //
1678   fNdigits   = 0;
1679   if (fDigits)   fDigits->Clear();
1680 }
1681
1682 //_____________________________________________________________________________
1683 void AliTPC::SetSecAL(Int_t sec)
1684 {
1685   //---------------------------------------------------
1686   // Activate/deactivate selection for lower sectors
1687   //---------------------------------------------------
1688
1689   //-----------------------------------------------------------------
1690   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1691   //-----------------------------------------------------------------
1692
1693   fSecAL = sec;
1694 }
1695
1696 //_____________________________________________________________________________
1697 void AliTPC::SetSecAU(Int_t sec)
1698 {
1699   //----------------------------------------------------
1700   // Activate/deactivate selection for upper sectors
1701   //---------------------------------------------------
1702
1703   //-----------------------------------------------------------------
1704   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1705   //-----------------------------------------------------------------
1706
1707   fSecAU = sec;
1708 }
1709
1710 //_____________________________________________________________________________
1711 void AliTPC::SetSecLows(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6)
1712 {
1713   //----------------------------------------
1714   // Select active lower sectors
1715   //----------------------------------------
1716
1717   //-----------------------------------------------------------------
1718   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1719   //-----------------------------------------------------------------
1720
1721   fSecLows[0] = s1;
1722   fSecLows[1] = s2;
1723   fSecLows[2] = s3;
1724   fSecLows[3] = s4;
1725   fSecLows[4] = s5;
1726   fSecLows[5] = s6;
1727 }
1728
1729 //_____________________________________________________________________________
1730 void AliTPC::SetSecUps(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6,
1731                        Int_t s7, Int_t s8 ,Int_t s9 ,Int_t s10,
1732                        Int_t s11 , Int_t s12)
1733 {
1734   //--------------------------------
1735   // Select active upper sectors
1736   //--------------------------------
1737
1738   //-----------------------------------------------------------------
1739   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1740   //-----------------------------------------------------------------
1741
1742   fSecUps[0] = s1;
1743   fSecUps[1] = s2;
1744   fSecUps[2] = s3;
1745   fSecUps[3] = s4;
1746   fSecUps[4] = s5;
1747   fSecUps[5] = s6;
1748   fSecUps[6] = s7;
1749   fSecUps[7] = s8;
1750   fSecUps[8] = s9;
1751   fSecUps[9] = s10;
1752   fSecUps[10] = s11;
1753   fSecUps[11] = s12;
1754 }
1755
1756 //_____________________________________________________________________________
1757 void AliTPC::SetSens(Int_t sens)
1758 {
1759
1760   //-------------------------------------------------------------
1761   // Activates/deactivates the sensitive strips at the center of
1762   // the pad row -- this is for the space-point resolution calculations
1763   //-------------------------------------------------------------
1764
1765   //-----------------------------------------------------------------
1766   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1767   //-----------------------------------------------------------------
1768
1769   fSens = sens;
1770 }
1771
1772
1773 void AliTPC::SetSide(Float_t side=0.)
1774 {
1775   // choice of the TPC side
1776
1777   fSide = side;
1778
1779 }
1780 //____________________________________________________________________________
1781 void AliTPC::SetGasMixt(Int_t nc,Int_t c1,Int_t c2,Int_t c3,Float_t p1,
1782                            Float_t p2,Float_t p3)
1783 {
1784
1785   // gax mixture definition
1786
1787  fNoComp = nc;
1788
1789  fMixtComp[0]=c1;
1790  fMixtComp[1]=c2;
1791  fMixtComp[2]=c3;
1792
1793  fMixtProp[0]=p1;
1794  fMixtProp[1]=p2;
1795  fMixtProp[2]=p3;
1796
1797
1798 }
1799 //_____________________________________________________________________________
1800
1801 void AliTPC::TransportElectron(Float_t *xyz, Int_t *index)
1802 {
1803   //
1804   // electron transport taking into account:
1805   // 1. diffusion,
1806   // 2.ExB at the wires
1807   // 3. nonisochronity
1808   //
1809   // xyz and index must be already transformed to system 1
1810   //
1811
1812   fTPCParam->Transform1to2(xyz,index);
1813
1814   //add diffusion
1815   Float_t driftl=xyz[2];
1816   if(driftl<0.01) driftl=0.01;
1817   driftl=TMath::Sqrt(driftl);
1818   Float_t sigT = driftl*(fTPCParam->GetDiffT());
1819   Float_t sigL = driftl*(fTPCParam->GetDiffL());
1820   xyz[0]=gRandom->Gaus(xyz[0],sigT);
1821   xyz[1]=gRandom->Gaus(xyz[1],sigT);
1822   xyz[2]=gRandom->Gaus(xyz[2],sigL);
1823
1824   // ExB
1825
1826   if (fTPCParam->GetMWPCReadout()==kTRUE){
1827     Float_t x1=xyz[0];
1828     fTPCParam->Transform2to2NearestWire(xyz,index);
1829     Float_t dx=xyz[0]-x1;
1830     xyz[1]+=dx*(fTPCParam->GetOmegaTau());
1831   }
1832   //add nonisochronity (not implemented yet)
1833
1834 }
1835 //_____________________________________________________________________________
1836 void AliTPC::Streamer(TBuffer &R__b)
1837 {
1838   //
1839   // Stream an object of class AliTPC.
1840   //
1841    if (R__b.IsReading()) {
1842       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1843       AliDetector::Streamer(R__b);
1844       R__b >> fTPCParam;
1845       if (R__v < 2) return;
1846       R__b >> fNsectors;
1847       R__b >> fHitType;
1848
1849    } else {
1850       R__b.WriteVersion(AliTPC::IsA());
1851       AliDetector::Streamer(R__b);
1852       R__b << fTPCParam;
1853       R__b << fNsectors;
1854       R__b << fHitType;
1855    }
1856 }
1857
1858 ClassImp(AliTPCdigit)
1859
1860 //_____________________________________________________________________________
1861 AliTPCdigit::AliTPCdigit(Int_t *tracks, Int_t *digits):
1862   AliDigit(tracks)
1863 {
1864   //
1865   // Creates a TPC digit object
1866   //
1867   fSector     = digits[0];
1868   fPadRow     = digits[1];
1869   fPad        = digits[2];
1870   fTime       = digits[3];
1871   fSignal     = digits[4];
1872 }
1873
1874
1875 ClassImp(AliTPChit)
1876
1877 //_____________________________________________________________________________
1878 AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
1879 AliHit(shunt,track)
1880 {
1881   //
1882   // Creates a TPC hit object
1883   //
1884   fSector     = vol[0];
1885   fPadRow     = vol[1];
1886   fX          = hits[0];
1887   fY          = hits[1];
1888   fZ          = hits[2];
1889   fQ          = hits[3];
1890 }
1891
1892
1893 //________________________________________________________________________
1894 // Additional code because of the AliTPCTrackHits
1895
1896 void AliTPC::MakeBranch2(Option_t *option)
1897 {
1898   //
1899   // Create a new branch in the current Root Tree
1900   // The branch of fHits is automatically split
1901   // MI change 14.09.2000
1902   if (fHitType&2==0) return;
1903   char branchname[10];
1904   sprintf(branchname,"%s2",GetName());
1905   //
1906   // Get the pointer to the header
1907   char *cH = strstr(option,"H");
1908   //
1909   if (fTrackHits   && gAlice->TreeH() && cH) {
1910     //    gAlice->TreeH()->Branch(branchname,&fTrackHits, fBufferSize);
1911     AliObjectBranch * branch = new AliObjectBranch(branchname,"AliTPCTrackHits",&fTrackHits,
1912                                                    gAlice->TreeH(),fBufferSize,1);
1913     gAlice->TreeH()->GetListOfBranches()->Add(branch);
1914     printf("* AliDetector::MakeBranch * Making Branch %s for trackhits\n",branchname);
1915   }
1916 }
1917
1918 void AliTPC::SetTreeAddress()
1919 {
1920   if (fHitType&1) AliDetector::SetTreeAddress();
1921   if (fHitType&2) SetTreeAddress2();
1922 }
1923
1924 void AliTPC::SetTreeAddress2()
1925 {
1926   //
1927   // Set branch address for the TrackHits Tree
1928   //
1929   TBranch *branch;
1930   char branchname[20];
1931   sprintf(branchname,"%s2",GetName());
1932   //
1933   // Branch address for hit tree
1934   TTree *treeH = gAlice->TreeH();
1935   if (treeH) {
1936     branch = treeH->GetBranch(branchname);
1937     if (branch) branch->SetAddress(&fTrackHits);
1938   }
1939 }
1940
1941 void AliTPC::FinishPrimary()
1942 {
1943   if (fTrackHits) fTrackHits->FlushHitStack();
1944 }
1945
1946
1947 void AliTPC::AddHit2(Int_t track, Int_t *vol, Float_t *hits)
1948 {
1949   //
1950   // add hit to the list
1951   TClonesArray &particles = *(gAlice->Particles());
1952   Int_t rtrack;
1953   if (fIshunt) {
1954     int primary = gAlice->GetPrimary(track);
1955     ((TParticle *)particles[primary])->SetBit(kKeepBit);
1956     rtrack=primary;
1957   } else {
1958     rtrack=track;
1959     gAlice->FlagTrack(track);
1960   }
1961   //AliTPChit *hit = (AliTPChit*)fHits->UncheckedAt(fNhits-1);
1962   //if (hit->fTrack!=rtrack)
1963   //  cout<<"bad track number\n";
1964   if (fTrackHits)
1965     fTrackHits->AddHitKartez(vol[0],rtrack, hits[0],
1966                              hits[1],hits[2],(Int_t)hits[3]);
1967 }
1968
1969 void AliTPC::ResetHits()
1970 {
1971   if (fHitType&1) AliDetector::ResetHits();
1972   if (fHitType&2) ResetHits2();
1973 }
1974
1975 void AliTPC::ResetHits2()
1976 {
1977   //
1978   //reset hits
1979   if (fTrackHits) fTrackHits->Clear();
1980 }
1981
1982 AliHit* AliTPC::FirstHit(Int_t track)
1983 {
1984   if (fHitType&2) return FirstHit2(track);
1985   return AliDetector::FirstHit(track);
1986 }
1987 AliHit* AliTPC::NextHit()
1988 {
1989   if (fHitType&2) return NextHit2();
1990   return AliDetector::NextHit();
1991 }
1992
1993 AliHit* AliTPC::FirstHit2(Int_t track)
1994 {
1995   //
1996   // Initialise the hit iterator
1997   // Return the address of the first hit for track
1998   // If track>=0 the track is read from disk
1999   // while if track<0 the first hit of the current
2000   // track is returned
2001   //
2002   if(track>=0) {
2003     gAlice->ResetHits();
2004     gAlice->TreeH()->GetEvent(track);
2005   }
2006   //
2007   if (fTrackHits) {
2008     fTrackHits->First();
2009     return fTrackHits->GetHit();
2010   }
2011   else return 0;
2012 }
2013
2014 AliHit* AliTPC::NextHit2()
2015 {
2016   //
2017   //Return the next hit for the current track
2018
2019   if (fTrackHits) {
2020     fTrackHits->Next();
2021     return fTrackHits->GetHit();
2022   }
2023   else
2024     return 0;
2025 }
2026
2027 void AliTPC::LoadPoints(Int_t)
2028 {
2029   //
2030   Int_t a = 0;
2031   /*  if(fHitType==1) return AliDetector::LoadPoints(a);
2032   LoadPoints2(a);
2033   */
2034   if(fHitType==1) AliDetector::LoadPoints(a);
2035   else LoadPoints2(a);
2036
2037   // LoadPoints3(a);
2038
2039 }
2040
2041
2042 void AliTPC::RemapTrackHitIDs(Int_t *map)
2043 {
2044   if (!fTrackHits) return;
2045   AliObjectArray * arr = fTrackHits->fTrackHitsInfo;
2046   for (UInt_t i=0;i<arr->GetSize();i++){
2047     AliTrackHitsInfo * info = (AliTrackHitsInfo *)(arr->At(i));
2048     info->fTrackID = map[info->fTrackID];
2049   }
2050
2051 }
2052
2053
2054 //_____________________________________________________________________________
2055 void AliTPC::LoadPoints2(Int_t)
2056 {
2057   //
2058   // Store x, y, z of all hits in memory
2059   //
2060   if (fTrackHits == 0) return;
2061   //
2062   Int_t nhits = fTrackHits->GetEntriesFast();
2063   if (nhits == 0) return;
2064   Int_t tracks = gAlice->GetNtrack();
2065   if (fPoints == 0) fPoints = new TObjArray(tracks);
2066   AliHit *ahit;
2067   //
2068   Int_t *ntrk=new Int_t[tracks];
2069   Int_t *limi=new Int_t[tracks];
2070   Float_t **coor=new Float_t*[tracks];
2071   for(Int_t i=0;i<tracks;i++) {
2072     ntrk[i]=0;
2073     coor[i]=0;
2074     limi[i]=0;
2075   }
2076   //
2077   AliPoints *points = 0;
2078   Float_t *fp=0;
2079   Int_t trk;
2080   Int_t chunk=nhits/4+1;
2081   //
2082   // Loop over all the hits and store their position
2083   //
2084   ahit = FirstHit2(-1);
2085   //for (Int_t hit=0;hit<nhits;hit++) {
2086   while (ahit){
2087     //    ahit = (AliHit*)fHits->UncheckedAt(hit);
2088     trk=ahit->GetTrack();
2089     if(ntrk[trk]==limi[trk]) {
2090       //
2091       // Initialise a new track
2092       fp=new Float_t[3*(limi[trk]+chunk)];
2093       if(coor[trk]) {
2094         memcpy(fp,coor[trk],sizeof(Float_t)*3*limi[trk]);
2095         delete [] coor[trk];
2096       }
2097       limi[trk]+=chunk;
2098       coor[trk] = fp;
2099     } else {
2100       fp = coor[trk];
2101     }
2102     fp[3*ntrk[trk]  ] = ahit->X();
2103     fp[3*ntrk[trk]+1] = ahit->Y();
2104     fp[3*ntrk[trk]+2] = ahit->Z();
2105     ntrk[trk]++;
2106     ahit = NextHit2();
2107   }
2108   //
2109   for(trk=0; trk<tracks; ++trk) {
2110     if(ntrk[trk]) {
2111       points = new AliPoints();
2112       points->SetMarkerColor(GetMarkerColor());
2113       points->SetMarkerSize(GetMarkerSize());
2114       points->SetDetector(this);
2115       points->SetParticle(trk);
2116       points->SetPolyMarker(ntrk[trk],coor[trk],GetMarkerStyle());
2117       fPoints->AddAt(points,trk);
2118       delete [] coor[trk];
2119       coor[trk]=0;
2120     }
2121   }
2122   delete [] coor;
2123   delete [] ntrk;
2124   delete [] limi;
2125 }
2126
2127
2128 //_____________________________________________________________________________
2129 void AliTPC::LoadPoints3(Int_t)
2130 {
2131   //
2132   // Store x, y, z of all hits in memory
2133   // - only intersection point with pad row
2134   if (fTrackHits == 0) return;
2135   //
2136   Int_t nhits = fTrackHits->GetEntriesFast();
2137   if (nhits == 0) return;
2138   Int_t tracks = gAlice->GetNtrack();
2139   if (fPoints == 0) fPoints = new TObjArray(2*tracks);
2140   fPoints->Expand(2*tracks);
2141   AliHit *ahit;
2142   //
2143   Int_t *ntrk=new Int_t[tracks];
2144   Int_t *limi=new Int_t[tracks];
2145   Float_t **coor=new Float_t*[tracks];
2146   for(Int_t i=0;i<tracks;i++) {
2147     ntrk[i]=0;
2148     coor[i]=0;
2149     limi[i]=0;
2150   }
2151   //
2152   AliPoints *points = 0;
2153   Float_t *fp=0;
2154   Int_t trk;
2155   Int_t chunk=nhits/4+1;
2156   //
2157   // Loop over all the hits and store their position
2158   //
2159   ahit = FirstHit2(-1);
2160   //for (Int_t hit=0;hit<nhits;hit++) {
2161
2162   Int_t lastrow = -1;
2163   while (ahit){
2164     //    ahit = (AliHit*)fHits->UncheckedAt(hit);
2165     trk=ahit->GetTrack();
2166     Float_t  x[3]={ahit->X(),ahit->Y(),ahit->Z()};
2167     Int_t    index[3]={1,((AliTPChit*)ahit)->fSector,0};
2168     Int_t    currentrow = fTPCParam->GetPadRow(x,index) ;
2169     if (currentrow!=lastrow){
2170       lastrow = currentrow;
2171       //later calculate intersection point
2172       if(ntrk[trk]==limi[trk]) {
2173         //
2174         // Initialise a new track
2175         fp=new Float_t[3*(limi[trk]+chunk)];
2176         if(coor[trk]) {
2177           memcpy(fp,coor[trk],sizeof(Float_t)*3*limi[trk]);
2178           delete [] coor[trk];
2179         }
2180         limi[trk]+=chunk;
2181         coor[trk] = fp;
2182       } else {
2183         fp = coor[trk];
2184       }
2185       fp[3*ntrk[trk]  ] = ahit->X();
2186       fp[3*ntrk[trk]+1] = ahit->Y();
2187       fp[3*ntrk[trk]+2] = ahit->Z();
2188       ntrk[trk]++;
2189     }
2190     ahit = NextHit2();
2191   }
2192
2193   //
2194   for(trk=0; trk<tracks; ++trk) {
2195     if(ntrk[trk]) {
2196       points = new AliPoints();
2197       points->SetMarkerColor(GetMarkerColor()+1);
2198       points->SetMarkerStyle(5);
2199       points->SetMarkerSize(0.2);
2200       points->SetDetector(this);
2201       points->SetParticle(trk);
2202       //      points->SetPolyMarker(ntrk[trk],coor[trk],GetMarkerStyle()20);
2203       points->SetPolyMarker(ntrk[trk],coor[trk],30);
2204       fPoints->AddAt(points,tracks+trk);
2205       delete [] coor[trk];
2206       coor[trk]=0;
2207     }
2208   }
2209   delete [] coor;
2210   delete [] ntrk;
2211   delete [] limi;
2212 }
2213
2214
2215
2216 void AliTPC::FindTrackHitsIntersection(TClonesArray * arr)
2217 {
2218
2219   //
2220   //fill clones array with intersection of current point with the
2221   //middle of the row
2222   Int_t sector;
2223   Int_t ipart;
2224
2225   const Int_t kcmaxhits=30000;
2226   TVector * xxxx = new TVector(kcmaxhits*4);
2227   TVector & xxx = *xxxx;
2228   Int_t maxhits = kcmaxhits;
2229
2230   //
2231   AliTPChit * tpcHit=0;
2232   tpcHit = (AliTPChit*)FirstHit2(-1);
2233   Int_t currentIndex=0;
2234   Int_t lastrow=-1;  //last writen row
2235
2236   while (tpcHit){
2237     if (tpcHit==0) continue;
2238     sector=tpcHit->fSector; // sector number
2239     ipart=tpcHit->Track();
2240
2241     //find row number
2242
2243     Float_t  x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
2244     Int_t    index[3]={1,sector,0};
2245     Int_t    currentrow = fTPCParam->GetPadRow(x,index) ;
2246     if (currentrow<0) continue;
2247     if (lastrow<0) lastrow=currentrow;
2248     if (currentrow==lastrow){
2249       if ( currentIndex>=maxhits){
2250         maxhits+=kcmaxhits;
2251         xxx.ResizeTo(4*maxhits);
2252       }
2253       xxx(currentIndex*4)=x[0];
2254       xxx(currentIndex*4+1)=x[1];
2255       xxx(currentIndex*4+2)=x[2];
2256       xxx(currentIndex*4+3)=tpcHit->fQ;
2257       currentIndex++;
2258     }
2259     else
2260       if (currentIndex>2){
2261         Float_t sumx=0;
2262         Float_t sumx2=0;
2263         Float_t sumx3=0;
2264         Float_t sumx4=0;
2265         Float_t sumy=0;
2266         Float_t sumxy=0;
2267         Float_t sumx2y=0;
2268         Float_t sumz=0;
2269         Float_t sumxz=0;
2270         Float_t sumx2z=0;
2271         Float_t sumq=0;
2272         for (Int_t index=0;index<currentIndex;index++){
2273           Float_t x,x2,x3,x4;
2274           x=x2=x3=x4=xxx(index*4);
2275           x2*=x;
2276           x3*=x2;
2277           x4*=x3;
2278           sumx+=x;
2279           sumx2+=x2;
2280           sumx3+=x3;
2281           sumx4+=x4;
2282           sumy+=xxx(index*4+1);
2283           sumxy+=xxx(index*4+1)*x;
2284           sumx2y+=xxx(index*4+1)*x2;
2285           sumz+=xxx(index*4+2);
2286           sumxz+=xxx(index*4+2)*x;
2287           sumx2z+=xxx(index*4+2)*x2;
2288           sumq+=xxx(index*4+3);
2289         }
2290         Float_t centralPad = (fTPCParam->GetNPads(sector,lastrow)-1)/2;
2291         Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
2292           sumx2*(sumx*sumx3-sumx2*sumx2);
2293
2294         Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
2295           sumx2*(sumxy*sumx3-sumx2y*sumx2);
2296         Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
2297           sumx2*(sumxz*sumx3-sumx2z*sumx2);
2298
2299         Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
2300           sumx2*(sumx*sumx2y-sumx2*sumxy);
2301         Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
2302           sumx2*(sumx*sumx2z-sumx2*sumxz);
2303
2304         Float_t y=detay/det+centralPad;
2305         Float_t z=detaz/det;
2306         Float_t by=detby/det; //y angle
2307         Float_t bz=detbz/det; //z angle
2308         sumy/=Float_t(currentIndex);
2309         sumz/=Float_t(currentIndex);
2310
2311         AliComplexCluster cl;
2312         cl.fX=z;
2313         cl.fY=y;
2314         cl.fQ=sumq;
2315         cl.fSigmaX2=bz;
2316         cl.fSigmaY2=by;
2317         cl.fTracks[0]=ipart;
2318
2319         AliTPCClustersRow * row = (fClustersArray->GetRow(sector,lastrow));
2320         if (row!=0) row->InsertCluster(&cl);
2321         currentIndex=0;
2322         lastrow=currentrow;
2323       } //end of calculating cluster for given row
2324
2325   } // end of loop over hits
2326   xxxx->Delete();
2327
2328
2329
2330 }