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 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 //  Time Projection Chamber                                                  //
  21 //  This class contains the basic functions for the Time Projection Chamber  //
  22 //  detector. Functions specific to one particular geometry are              //
  23 //  contained in the derived classes                                         //
  24 //                                                                           //
  25 //Begin_Html
  26 /*
  27 <img src="picts/AliTPCClass.gif">
  28 */
  29 //End_Html
  30 //                                                                           //
  31 //                                                                          //
  32 ///////////////////////////////////////////////////////////////////////////////
  33
  34 //
  35
  36 #include <Riostream.h>
  37 #include <stdlib.h>
  38
  39 #include <TFile.h>
  40 #include <TGeometry.h>
  41 #include <TInterpreter.h>
  42 #include <TMath.h>
  43 #include <TMatrixF.h>
  44 #include <TVector.h>
  45 #include <TNode.h>
  46 #include <TObjectTable.h>
  47 #include <TParticle.h>
  48 #include <TROOT.h>
  49 #include <TRandom.h>
  50 #include <TSystem.h>
  51 #include <TTUBS.h>
  52 #include <TTree.h>
  53 #include <TVirtualMC.h>
  54 #include <TString.h>
  55 #include <TF2.h>
  56 #include <TStopwatch.h>
  57
  58 #include "AliDigits.h"
  59 #include "AliMagF.h"
  60 #include "AliPoints.h"
  61 #include "AliRun.h"
  62 #include "AliRunLoader.h"
  63 #include "AliSimDigits.h"
  64 #include "AliTPC.h"
  65 #include "AliTPC.h"
  66 #include "AliTPCDigitsArray.h"
  67 #include "AliTPCLoader.h"
  68 #include "AliTPCPRF2D.h"
  69 #include "AliTPCParamSR.h"
  70 #include "AliTPCRF1D.h"
  71 //#include "AliTPCTrackHits.h"
  72 #include "AliTPCTrackHitsV2.h"
  73 #include "AliTrackReference.h"
  74 #include "AliMC.h"
  75 #include "AliStack.h"
  76 #include "AliTPCDigitizer.h"
  77 #include "AliTPCBuffer.h"
  78 #include "AliTPCDDLRawData.h"
  79 #include "AliLog.h"
  80 #include "AliTPCcalibDB.h"
  81 #include "AliTPCCalPad.h"
  82 #include "AliTPCCalROC.h"
  83 #include "AliTPCExB.h"
  84 #include "AliRawReader.h"
  85 #include "AliTPCRawStream.h"
  86
  87 ClassImp(AliTPC)
  88 //_____________________________________________________________________________
  89   AliTPC::AliTPC():AliDetector(),
  90                    fDefaults(0),
  91                    fSens(0),
  92                    fNsectors(0),
  93                    fDigitsArray(0),
  94                    fTPCParam(0),
  95                    fTrackHits(0),
  96                    fHitType(0),
  97                    fDigitsSwitch(0),
  98                    fSide(0),
  99                    fPrimaryIonisation(0),
 100                    fNoiseDepth(0),
 101                    fNoiseTable(0),
 102                    fCurrentNoise(0),
 103                    fActiveSectors(0)
 104
 105
 106 {
 107   //
 108   // Default constructor
 109   //
 110   fIshunt   = 0;
 111
 112   //  fTrackHitsOld = 0;
 113 #if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,1)
 114   fHitType = 4; // ROOT containers
 115 #else
 116   fHitType = 2; //default CONTAINERS - based on ROOT structure
 117 #endif
 118 }
 119
 120 //_____________________________________________________________________________
 121 AliTPC::AliTPC(const char *name, const char *title)
 122   : AliDetector(name,title),
 123                    fDefaults(0),
 124                    fSens(0),
 125                    fNsectors(0),
 126                    fDigitsArray(0),
 127                    fTPCParam(0),
 128                    fTrackHits(0),
 129                    fHitType(0),
 130                    fDigitsSwitch(0),
 131                    fSide(0),
 132                    fPrimaryIonisation(0),
 133                    fNoiseDepth(0),
 134                    fNoiseTable(0),
 135                    fCurrentNoise(0),
 136                    fActiveSectors(0)
 137
 138 {
 139   //
 140   // Standard constructor
 141   //
 142
 143   //
 144   // Initialise arrays of hits and digits
 145   fHits     = new TClonesArray("AliTPChit",  176);
 146   gAlice->GetMCApp()->AddHitList(fHits);
 147   //
 148   fTrackHits = new AliTPCTrackHitsV2;
 149   fTrackHits->SetHitPrecision(0.002);
 150   fTrackHits->SetStepPrecision(0.003);
 151   fTrackHits->SetMaxDistance(100);
 152
 153   //fTrackHitsOld = new AliTPCTrackHits;  //MI - 13.09.2000
 154   //fTrackHitsOld->SetHitPrecision(0.002);
 155   //fTrackHitsOld->SetStepPrecision(0.003);
 156   //fTrackHitsOld->SetMaxDistance(100);
 157
 158
 159 #if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,1)
 160   fHitType = 4; // ROOT containers
 161 #else
 162   fHitType = 2;
 163 #endif
 164
 165
 166
 167   //
 168   fIshunt     =  0;
 169   //
 170   // Initialise color attributes
 171   //PH SetMarkerColor(kYellow);
 172
 173   //
 174   //  Set TPC parameters
 175   //
 176
 177
 178   if (!strcmp(title,"Default")) {
 179     fTPCParam = new AliTPCParamSR;
 180   } else {
 181     AliWarning("In Config.C you must set non-default parameters.");
 182     fTPCParam=0;
 183   }
 184 }
 185
 186 //_____________________________________________________________________________
 187 AliTPC::~AliTPC()
 188 {
 189   //
 190   // TPC destructor
 191   //
 192
 193   fIshunt   = 0;
 194   delete fHits;
 195   delete fDigits;
 196   delete fTPCParam;
 197   delete fTrackHits; //MI 15.09.2000
 198   //  delete fTrackHitsOld; //MI 10.12.2001
 199
 200   fDigitsArray = 0x0;
 201   delete [] fNoiseTable;
 202   delete [] fActiveSectors;
 203
 204 }
 205
 206 //_____________________________________________________________________________
 207 void AliTPC::AddHit(Int_t track, Int_t *vol, Float_t *hits)
 208 {
 209   //
 210   // Add a hit to the list
 211   //
 212   if (fHitType&1){
 213     TClonesArray &lhits = *fHits;
 214     new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits);
 215   }
 216   if (fHitType>1)
 217     AddHit2(track,vol,hits);
 218 }
 219
 220 //_____________________________________________________________________________
 221 void AliTPC::BuildGeometry()
 222 {
 223
 224   //
 225   // Build TPC ROOT TNode geometry for the event display
 226   //
 227   TNode *nNode, *nTop;
 228   TTUBS *tubs;
 229   Int_t i;
 230   const int kColorTPC=19;
 231   char name[5], title[25];
 232   const Double_t kDegrad=TMath::Pi()/180;
 233   const Double_t kRaddeg=180./TMath::Pi();
 234
 235
 236   Float_t innerOpenAngle = fTPCParam->GetInnerAngle();
 237   Float_t outerOpenAngle = fTPCParam->GetOuterAngle();
 238
 239   Float_t innerAngleShift = fTPCParam->GetInnerAngleShift();
 240   Float_t outerAngleShift = fTPCParam->GetOuterAngleShift();
 241
 242   Int_t nLo = fTPCParam->GetNInnerSector()/2;
 243   Int_t nHi = fTPCParam->GetNOuterSector()/2;
 244
 245   const Double_t kloAng = (Double_t)TMath::Nint(innerOpenAngle*kRaddeg);
 246   const Double_t khiAng = (Double_t)TMath::Nint(outerOpenAngle*kRaddeg);
 247   const Double_t kloAngSh = (Double_t)TMath::Nint(innerAngleShift*kRaddeg);
 248   const Double_t khiAngSh = (Double_t)TMath::Nint(outerAngleShift*kRaddeg);
 249
 250
 251   const Double_t kloCorr = 1/TMath::Cos(0.5*kloAng*kDegrad);
 252   const Double_t khiCorr = 1/TMath::Cos(0.5*khiAng*kDegrad);
 253
 254   Double_t rl,ru;
 255
 256
 257   //
 258   // Get ALICE top node
 259   //
 260
 261   nTop=gAlice->GetGeometry()->GetNode("alice");
 262
 263   //  inner sectors
 264
 265   rl = fTPCParam->GetInnerRadiusLow();
 266   ru = fTPCParam->GetInnerRadiusUp();
 267
 268
 269   for(i=0;i<nLo;i++) {
 270     sprintf(name,"LS%2.2d",i);
 271     name[4]='\0';
 272     sprintf(title,"TPC low sector %3d",i);
 273     title[24]='\0';
 274
 275     tubs = new TTUBS(name,title,"void",rl*kloCorr,ru*kloCorr,250.,
 276                      kloAng*(i-0.5)+kloAngSh,kloAng*(i+0.5)+kloAngSh);
 277     tubs->SetNumberOfDivisions(1);
 278     nTop->cd();
 279     nNode = new TNode(name,title,name,0,0,0,"");
 280     nNode->SetLineColor(kColorTPC);
 281     fNodes->Add(nNode);
 282   }
 283
 284   // Outer sectors
 285
 286   rl = fTPCParam->GetOuterRadiusLow();
 287   ru = fTPCParam->GetOuterRadiusUp();
 288
 289   for(i=0;i<nHi;i++) {
 290     sprintf(name,"US%2.2d",i);
 291     name[4]='\0';
 292     sprintf(title,"TPC upper sector %d",i);
 293     title[24]='\0';
 294     tubs = new TTUBS(name,title,"void",rl*khiCorr,ru*khiCorr,250,
 295                      khiAng*(i-0.5)+khiAngSh,khiAng*(i+0.5)+khiAngSh);
 296     tubs->SetNumberOfDivisions(1);
 297     nTop->cd();
 298     nNode = new TNode(name,title,name,0,0,0,"");
 299     nNode->SetLineColor(kColorTPC);
 300     fNodes->Add(nNode);
 301   }
 302
 303 }
 304
 305 //_____________________________________________________________________________
 306 void AliTPC::CreateMaterials()
 307 {
 308   //-----------------------------------------------
 309   // Create Materials for for TPC simulations
 310   //-----------------------------------------------
 311
 312   //-----------------------------------------------------------------
 313   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
 314   //-----------------------------------------------------------------
 315
 316    Int_t iSXFLD=gAlice->Field()->Integ();
 317   Float_t sXMGMX=gAlice->Field()->Max();
 318
 319   Float_t amat[5]; // atomic numbers
 320   Float_t zmat[5]; // z
 321   Float_t wmat[5]; // proportions
 322
 323   Float_t density;
 324
 325
 326
 327   //***************** Gases *************************
 328
 329
 330   //--------------------------------------------------------------
 331   // gases - air and CO2
 332   //--------------------------------------------------------------
 333
 334   // CO2
 335
 336   amat[0]=12.011;
 337   amat[1]=15.9994;
 338
 339   zmat[0]=6.;
 340   zmat[1]=8.;
 341
 342   wmat[0]=0.2729;
 343   wmat[1]=0.7271;
 344
 345   density=0.001977;
 346
 347
 348   AliMixture(10,"CO2",amat,zmat,density,2,wmat);
 349   //
 350   // Air
 351   //
 352   amat[0]=15.9994;
 353   amat[1]=14.007;
 354   //
 355   zmat[0]=8.;
 356   zmat[1]=7.;
 357   //
 358   wmat[0]=0.233;
 359   wmat[1]=0.767;
 360   //
 361   density=0.001205;
 362
 363   AliMixture(11,"Air",amat,zmat,density,2,wmat);
 364
 365   //----------------------------------------------------------------
 366   // drift gases
 367   //----------------------------------------------------------------
 368
 369
 370   // Drift gases 1 - nonsensitive, 2 - sensitive
 371   // Ne-CO2-N (85-10-5)
 372
 373   amat[0]= 20.18;
 374   amat[1]=12.011;
 375   amat[2]=15.9994;
 376   amat[3]=14.007;
 377
 378   zmat[0]= 10.;
 379   zmat[1]=6.;
 380   zmat[2]=8.;
 381   zmat[3]=7.;
 382
 383   wmat[0]=0.7707;
 384   wmat[1]=0.0539;
 385   wmat[2]=0.1438;
 386   wmat[3]=0.0316;
 387
 388   density=0.0010252;
 389
 390   AliMixture(12,"Ne-CO2-N-1",amat,zmat,density,4,wmat);
 391   AliMixture(13,"Ne-CO2-N-2",amat,zmat,density,4,wmat);
 392   AliMixture(30,"Ne-CO2-N-3",amat,zmat,density,4,wmat);
 393   //----------------------------------------------------------------------
 394   //               solid materials
 395   //----------------------------------------------------------------------
 396
 397
 398   // Kevlar C14H22O2N2
 399
 400   amat[0] = 12.011;
 401   amat[1] = 1.;
 402   amat[2] = 15.999;
 403   amat[3] = 14.006;
 404
 405   zmat[0] = 6.;
 406   zmat[1] = 1.;
 407   zmat[2] = 8.;
 408   zmat[3] = 7.;
 409
 410   wmat[0] = 14.;
 411   wmat[1] = 22.;
 412   wmat[2] = 2.;
 413   wmat[3] = 2.;
 414
 415   density = 1.45;
 416
 417   AliMixture(14,"Kevlar",amat,zmat,density,-4,wmat);
 418
 419   // NOMEX
 420
 421   amat[0] = 12.011;
 422   amat[1] = 1.;
 423   amat[2] = 15.999;
 424   amat[3] = 14.006;
 425
 426   zmat[0] = 6.;
 427   zmat[1] = 1.;
 428   zmat[2] = 8.;
 429   zmat[3] = 7.;
 430
 431   wmat[0] = 14.;
 432   wmat[1] = 22.;
 433   wmat[2] = 2.;
 434   wmat[3] = 2.;
 435
 436   density = 0.029;
 437
 438   AliMixture(15,"NOMEX",amat,zmat,density,-4,wmat);
 439
 440   // Makrolon C16H18O3
 441
 442   amat[0] = 12.011;
 443   amat[1] = 1.;
 444   amat[2] = 15.999;
 445
 446   zmat[0] = 6.;
 447   zmat[1] = 1.;
 448   zmat[2] = 8.;
 449
 450   wmat[0] = 16.;
 451   wmat[1] = 18.;
 452   wmat[2] = 3.;
 453
 454   density = 1.2;
 455
 456   AliMixture(16,"Makrolon",amat,zmat,density,-3,wmat);
 457
 458   // Tedlar C2H3F
 459
 460   amat[0] = 12.011;
 461   amat[1] = 1.;
 462   amat[2] = 18.998;
 463
 464   zmat[0] = 6.;
 465   zmat[1] = 1.;
 466   zmat[2] = 9.;
 467
 468   wmat[0] = 2.;
 469   wmat[1] = 3.;
 470   wmat[2] = 1.;
 471
 472   density = 1.71;
 473
 474   AliMixture(17, "Tedlar",amat,zmat,density,-3,wmat);
 475
 476   // Mylar C5H4O2
 477
 478   amat[0]=12.011;
 479   amat[1]=1.;
 480   amat[2]=15.9994;
 481
 482   zmat[0]=6.;
 483   zmat[1]=1.;
 484   zmat[2]=8.;
 485
 486   wmat[0]=5.;
 487   wmat[1]=4.;
 488   wmat[2]=2.;
 489
 490   density = 1.39;
 491
 492   AliMixture(18, "Mylar",amat,zmat,density,-3,wmat);
 493   // material for "prepregs"
 494   // Epoxy - C14 H20 O3
 495   // Quartz SiO2
 496   // Carbon C
 497   // prepreg1 60% C-fiber, 40% epoxy (vol)
 498   amat[0]=12.011;
 499   amat[1]=1.;
 500   amat[2]=15.994;
 501
 502   zmat[0]=6.;
 503   zmat[1]=1.;
 504   zmat[2]=8.;
 505
 506   wmat[0]=0.923;
 507   wmat[1]=0.023;
 508   wmat[2]=0.054;
 509
 510   density=1.859;
 511
 512   AliMixture(19, "Prepreg1",amat,zmat,density,3,wmat);
 513
 514   //prepreg2 60% glass-fiber, 40% epoxy
 515
 516   amat[0]=12.01;
 517   amat[1]=1.;
 518   amat[2]=15.994;
 519   amat[3]=28.086;
 520
 521   zmat[0]=6.;
 522   zmat[1]=1.;
 523   zmat[2]=8.;
 524   zmat[3]=14.;
 525
 526   wmat[0]=0.194;
 527   wmat[1]=0.023;
 528   wmat[2]=0.443;
 529   wmat[3]=0.340;
 530
 531   density=1.82;
 532
 533   AliMixture(20, "Prepreg2",amat,zmat,density,4,wmat);
 534
 535   //prepreg3 50% glass-fiber, 50% epoxy
 536
 537   amat[0]=12.01;
 538   amat[1]=1.;
 539   amat[2]=15.994;
 540   amat[3]=28.086;
 541
 542   zmat[0]=6.;
 543   zmat[1]=1.;
 544   zmat[2]=8.;
 545   zmat[3]=14.;
 546
 547   wmat[0]=0.225;
 548   wmat[1]=0.03;
 549   wmat[2]=0.443;
 550   wmat[3]=0.3;
 551
 552   density=1.163;
 553
 554   AliMixture(21, "Prepreg3",amat,zmat,density,4,wmat);
 555
 556   // G10 60% SiO2 40% epoxy
 557
 558   amat[0]=12.01;
 559   amat[1]=1.;
 560   amat[2]=15.994;
 561   amat[3]=28.086;
 562
 563   zmat[0]=6.;
 564   zmat[1]=1.;
 565   zmat[2]=8.;
 566   zmat[3]=14.;
 567
 568   wmat[0]=0.194;
 569   wmat[1]=0.023;
 570   wmat[2]=0.443;
 571   wmat[3]=0.340;
 572
 573   density=1.7;
 574
 575   AliMixture(22, "G10",amat,zmat,density,4,wmat);
 576
 577   // Al
 578
 579   amat[0] = 26.98;
 580   zmat[0] = 13.;
 581
 582   density = 2.7;
 583
 584   AliMaterial(23,"Al",amat[0],zmat[0],density,999.,999.);
 585
 586   // Si (for electronics
 587
 588   amat[0] = 28.086;
 589   zmat[0] = 14.;
 590
 591   density = 2.33;
 592
 593   AliMaterial(24,"Si",amat[0],zmat[0],density,999.,999.);
 594
 595   // Cu
 596
 597   amat[0] = 63.546;
 598   zmat[0] = 29.;
 599
 600   density = 8.96;
 601
 602   AliMaterial(25,"Cu",amat[0],zmat[0],density,999.,999.);
 603
 604   // Epoxy - C14 H20 O3
 605
 606   amat[0]=12.011;
 607   amat[1]=1.;
 608   amat[2]=15.9994;
 609
 610   zmat[0]=6.;
 611   zmat[1]=1.;
 612   zmat[2]=8.;
 613
 614   wmat[0]=14.;
 615   wmat[1]=20.;
 616   wmat[2]=3.;
 617
 618   density=1.25;
 619
 620   AliMixture(26,"Epoxy",amat,zmat,density,-3,wmat);
 621
 622   // Plexiglas  C5H8O2
 623
 624   amat[0]=12.011;
 625   amat[1]=1.;
 626   amat[2]=15.9994;
 627
 628   zmat[0]=6.;
 629   zmat[1]=1.;
 630   zmat[2]=8.;
 631
 632   wmat[0]=5.;
 633   wmat[1]=8.;
 634   wmat[2]=2.;
 635
 636   density=1.18;
 637
 638   AliMixture(27,"Plexiglas",amat,zmat,density,-3,wmat);
 639
 640   // Carbon
 641
 642   amat[0]=12.011;
 643   zmat[0]=6.;
 644   density= 2.265;
 645
 646   AliMaterial(28,"C",amat[0],zmat[0],density,999.,999.);
 647
 648   // Fe (steel for the inner heat screen)
 649
 650   amat[0]=55.845;
 651
 652   zmat[0]=26.;
 653
 654   density=7.87;
 655
 656   AliMaterial(29,"Fe",amat[0],zmat[0],density,999.,999.);
 657
 658   //----------------------------------------------------------
 659   // tracking media for gases
 660   //----------------------------------------------------------
 661
 662   AliMedium(0, "Air", 11, 0, iSXFLD, sXMGMX, 10., 999., .1, .01, .1);
 663   AliMedium(1, "Ne-CO2-N-1", 12, 0, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
 664   AliMedium(2, "Ne-CO2-N-2", 13, 1, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
 665   AliMedium(3,"CO2",10,0, iSXFLD, sXMGMX, 10., 999.,.1, .001, .001);
 666   AliMedium(20, "Ne-CO2-N-3", 30, 1, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
 667   //-----------------------------------------------------------
 668   // tracking media for solids
 669   //-----------------------------------------------------------
 670
 671   AliMedium(4,"Al",23,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 672   AliMedium(5,"Kevlar",14,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 673   AliMedium(6,"Nomex",15,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 674   AliMedium(7,"Makrolon",16,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 675   AliMedium(8,"Mylar",18,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 676   AliMedium(9,"Tedlar",17,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 677   //
 678   AliMedium(10,"Prepreg1",19,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 679   AliMedium(11,"Prepreg2",20,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 680   AliMedium(12,"Prepreg3",21,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 681   AliMedium(13,"Epoxy",26,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 682
 683   AliMedium(14,"Cu",25,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 684   AliMedium(15,"Si",24,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 685   AliMedium(16,"G10",22,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 686   AliMedium(17,"Plexiglas",27,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 687   AliMedium(18,"Steel",29,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 688
 689 }
 690
 691 void AliTPC::GenerNoise(Int_t tablesize)
 692 {
 693   //
 694   //Generate table with noise
 695   //
 696   if (fTPCParam==0) {
 697     // error message
 698     fNoiseDepth=0;
 699     return;
 700   }
 701   if (fNoiseTable)  delete[] fNoiseTable;
 702   fNoiseTable = new Float_t[tablesize];
 703   fNoiseDepth = tablesize;
 704   fCurrentNoise =0; //!index of the noise in  the noise table
 705
 706   Float_t norm = fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac();
 707   for (Int_t i=0;i<tablesize;i++) fNoiseTable[i]= gRandom->Gaus(0,norm);
 708 }
 709
 710 Float_t AliTPC::GetNoise()
 711 {
 712   // get noise from table
 713   //  if ((fCurrentNoise%10)==0)
 714   //  fCurrentNoise= gRandom->Rndm()*fNoiseDepth;
 715   if (fCurrentNoise>=fNoiseDepth) fCurrentNoise=0;
 716   return fNoiseTable[fCurrentNoise++];
 717   //gRandom->Gaus(0, fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac());
 718 }
 719
 720
 721 Bool_t  AliTPC::IsSectorActive(Int_t sec) const
 722 {
 723   //
 724   // check if the sector is active
 725   if (!fActiveSectors) return kTRUE;
 726   else return fActiveSectors[sec];
 727 }
 728
 729 void    AliTPC::SetActiveSectors(Int_t * sectors, Int_t n)
 730 {
 731   // activate interesting sectors
 732   SetTreeAddress();//just for security
 733   if (!fActiveSectors) fActiveSectors = new Bool_t[fTPCParam->GetNSector()];
 734   for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kFALSE;
 735   for (Int_t i=0;i<n;i++)
 736     if ((sectors[i]>=0) && sectors[i]<fTPCParam->GetNSector())  fActiveSectors[sectors[i]]=kTRUE;
 737
 738 }
 739
 740 void    AliTPC::SetActiveSectors(Int_t flag)
 741 {
 742   //
 743   // activate sectors which were hitted by tracks
 744   //loop over tracks
 745   SetTreeAddress();//just for security
 746   if (fHitType==0) return;  // if Clones hit - not short volume ID information
 747   if (!fActiveSectors) fActiveSectors = new Bool_t[fTPCParam->GetNSector()];
 748   if (flag) {
 749     for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kTRUE;
 750     return;
 751   }
 752   for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kFALSE;
 753   TBranch * branch=0;
 754   if (TreeH() == 0x0)
 755    {
 756      AliFatal("Can not find TreeH in folder");
 757      return;
 758    }
 759   if (fHitType>1) branch = TreeH()->GetBranch("TPC2");
 760   else branch = TreeH()->GetBranch("TPC");
 761   Stat_t ntracks = TreeH()->GetEntries();
 762   // loop over all hits
 763   AliDebug(1,Form("Got %d tracks",ntracks));
 764
 765   for(Int_t track=0;track<ntracks;track++) {
 766     ResetHits();
 767     //
 768     if (fTrackHits && fHitType&4) {
 769       TBranch * br1 = TreeH()->GetBranch("fVolumes");
 770       TBranch * br2 = TreeH()->GetBranch("fNVolumes");
 771       br1->GetEvent(track);
 772       br2->GetEvent(track);
 773       Int_t *volumes = fTrackHits->GetVolumes();
 774       for (Int_t j=0;j<fTrackHits->GetNVolumes(); j++) {
 775         if (volumes[j]>-1 && volumes[j]<fTPCParam->GetNSector()) {
 776           fActiveSectors[volumes[j]]=kTRUE;
 777         }
 778         else {
 779             AliError(Form("Volume %d -> sector number %d is outside (0..%d)",
 780                           j,
 781                           volumes[j],
 782                           fTPCParam->GetNSector()));
 783         }
 784       }
 785     }
 786
 787     //
 788 //     if (fTrackHitsOld && fHitType&2) {
 789 //       TBranch * br = TreeH()->GetBranch("fTrackHitsInfo");
 790 //       br->GetEvent(track);
 791 //       AliObjectArray * ar = fTrackHitsOld->fTrackHitsInfo;
 792 //       for (UInt_t j=0;j<ar->GetSize();j++){
 793 //      fActiveSectors[((AliTrackHitsInfo*)ar->At(j))->fVolumeID] =kTRUE;
 794 //       }
 795 //     }
 796   }
 797 }
 798
 799
 800
 801
 802 //_____________________________________________________________________________
 803 void AliTPC::Digits2Raw()
 804 {
 805 // convert digits of the current event to raw data
 806
 807   static const Int_t kThreshold = 0;
 808
 809   fLoader->LoadDigits();
 810   TTree* digits = fLoader->TreeD();
 811   if (!digits) {
 812     AliError("No digits tree");
 813     return;
 814   }
 815
 816   AliSimDigits digarr;
 817   AliSimDigits* digrow = &digarr;
 818   digits->GetBranch("Segment")->SetAddress(&digrow);
 819
 820   const char* fileName = "AliTPCDDL.dat";
 821   AliTPCBuffer* buffer  = new AliTPCBuffer(fileName);
 822   //Verbose level
 823   // 0: Silent
 824   // 1: cout messages
 825   // 2: txt files with digits
 826   //BE CAREFUL, verbose level 2 MUST be used only for debugging and
 827   //it is highly suggested to use this mode only for debugging digits files
 828   //reasonably small, because otherwise the size of the txt files can reach
 829   //quickly several MB wasting time and disk space.
 830   buffer->SetVerbose(0);
 831
 832   Int_t nEntries = Int_t(digits->GetEntries());
 833   Int_t previousSector = -1;
 834   Int_t subSector = 0;
 835   for (Int_t i = 0; i < nEntries; i++) {
 836     digits->GetEntry(i);
 837     Int_t sector, row;
 838     fTPCParam->AdjustSectorRow(digarr.GetID(), sector, row);
 839     if(previousSector != sector) {
 840       subSector = 0;
 841       previousSector = sector;
 842     }
 843
 844     if (sector < 36) { //inner sector [0;35]
 845       if (row != 30) {
 846         //the whole row is written into the output file
 847         buffer->WriteRowBinary(kThreshold, digrow, 0, 0, 0,
 848                                sector, subSector, row);
 849       } else {
 850         //only the pads in the range [37;48] are written into the output file
 851         buffer->WriteRowBinary(kThreshold, digrow, 37, 48, 1,
 852                                sector, subSector, row);
 853         subSector = 1;
 854         //only the pads outside the range [37;48] are written into the output file
 855         buffer->WriteRowBinary(kThreshold, digrow, 37, 48, 2,
 856                                sector, subSector, row);
 857       }//end else
 858
 859     } else { //outer sector [36;71]
 860       if (row == 54) subSector = 2;
 861       if ((row != 27) && (row != 76)) {
 862         buffer->WriteRowBinary(kThreshold, digrow, 0, 0, 0,
 863                                sector, subSector, row);
 864       } else if (row == 27) {
 865         //only the pads outside the range [43;46] are written into the output file
 866         buffer->WriteRowBinary(kThreshold, digrow, 43, 46, 2,
 867                                  sector, subSector, row);
 868         subSector = 1;
 869         //only the pads in the range [43;46] are written into the output file
 870         buffer->WriteRowBinary(kThreshold, digrow, 43, 46, 1,
 871                                  sector, subSector, row);
 872       } else if (row == 76) {
 873         //only the pads outside the range [33;88] are written into the output file
 874         buffer->WriteRowBinary(kThreshold, digrow, 33, 88, 2,
 875                                sector, subSector, row);
 876         subSector = 3;
 877         //only the pads in the range [33;88] are written into the output file
 878         buffer->WriteRowBinary(kThreshold, digrow, 33, 88, 1,
 879                                  sector, subSector, row);
 880       }
 881     }//end else
 882   }//end for
 883
 884   delete buffer;
 885   fLoader->UnloadDigits();
 886
 887   AliTPCDDLRawData rawWriter;
 888   rawWriter.SetVerbose(0);
 889
 890   rawWriter.RawData(fileName);
 891   gSystem->Unlink(fileName);
 892
 893 }
 894
 895
 896 //_____________________________________________________________________________
 897 Bool_t AliTPC::Raw2SDigits(AliRawReader* rawReader){
 898   // Converts the TPC raw data into summable digits
 899   // The method is used for merging simulated and
 900   // real data events
 901   if (fLoader->TreeS() == 0x0 ) {
 902     fLoader->MakeTree("S");
 903   }
 904
 905   if(fDefaults == 0) SetDefaults();  // check if the parameters are set
 906
 907   //setup TPCDigitsArray
 908   if(GetDigitsArray()) delete GetDigitsArray();
 909
 910   AliTPCDigitsArray *arr = new AliTPCDigitsArray;
 911   arr->SetClass("AliSimDigits");
 912   arr->Setup(fTPCParam);
 913   arr->MakeTree(fLoader->TreeS());
 914
 915   SetDigitsArray(arr);
 916
 917   // set zero suppression to "0"
 918   fTPCParam->SetZeroSup(0);
 919
 920   // Loop over sectors
 921   const Int_t kmaxTime = fTPCParam->GetMaxTBin();
 922   const Int_t kNIS = fTPCParam->GetNInnerSector();
 923   const Int_t kNOS = fTPCParam->GetNOuterSector();
 924   const Int_t kNS = kNIS + kNOS;
 925
 926   Short_t** allBins = NULL; //array which contains the data for one sector
 927
 928   for(Int_t iSector = 0; iSector < kNS; iSector++) {
 929
 930     Int_t nRows = fTPCParam->GetNRow(iSector);
 931     Int_t nDDLs = 0, indexDDL = 0;
 932     if (iSector < kNIS) {
 933       nDDLs = 2;
 934       indexDDL = iSector * 2;
 935     }
 936     else {
 937       nDDLs = 4;
 938       indexDDL = (iSector-kNIS) * 4 + kNIS * 2;
 939     }
 940
 941     // Loas the raw data for corresponding DDLs
 942     rawReader->Reset();
 943     AliTPCRawStream input(rawReader);
 944     rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
 945
 946     // Alocate and init the array with the sector data
 947     allBins = new Short_t*[nRows];
 948     for (Int_t iRow = 0; iRow < nRows; iRow++) {
 949       Int_t maxPad = fTPCParam->GetNPads(iSector,iRow);
 950       Int_t maxBin = kmaxTime*maxPad;
 951       allBins[iRow] = new Short_t[maxBin];
 952       memset(allBins[iRow],0,sizeof(Short_t)*maxBin);
 953     }
 954
 955     // Begin loop over altro data
 956     while (input.Next()) {
 957
 958       if (input.GetSector() != iSector)
 959         AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",iSector,input.GetSector()));
 960
 961       Int_t iRow = input.GetRow();
 962       if (iRow < 0 || iRow >= nRows)
 963         AliFatal(Form("Pad-row index (%d) outside the range (%d -> %d) !",
 964                       iRow, 0, nRows -1));
 965       Int_t iPad = input.GetPad();
 966
 967       Int_t maxPad = fTPCParam->GetNPads(iSector,iRow);
 968
 969       if (iPad < 0 || iPad >= maxPad)
 970         AliFatal(Form("Pad index (%d) outside the range (%d -> %d) !",
 971                       iPad, 0, maxPad -1));
 972
 973       Int_t iTimeBin = input.GetTime();
 974       if ( iTimeBin < 0 || iTimeBin >= kmaxTime)
 975         AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
 976                       iTimeBin, 0, kmaxTime -1));
 977
 978       Int_t maxBin = kmaxTime*maxPad;
 979
 980       if (((iPad*kmaxTime+iTimeBin) >= maxBin) ||
 981           ((iPad*kmaxTime+iTimeBin) < 0))
 982         AliFatal(Form("Index outside the allowed range"
 983                       " Sector=%d Row=%d Pad=%d Timebin=%d"
 984                       " (Max.index=%d)",iSector,iRow,iPad,iTimeBin,maxBin));
 985
 986       allBins[iRow][iPad*kmaxTime+iTimeBin] = input.GetSignal();
 987
 988     } // End loop over altro data
 989
 990     // Now fill the digits array
 991     if (fDigitsArray->GetTree()==0) {
 992       AliFatal("Tree not set in fDigitsArray");
 993     }
 994
 995     for (Int_t iRow = 0; iRow < nRows; iRow++) {
 996       AliDigits * dig = fDigitsArray->CreateRow(iSector,iRow);
 997
 998       Int_t maxPad = fTPCParam->GetNPads(iSector,iRow);
 999       for(Int_t iPad = 0; iPad < maxPad; iPad++) {
1000         for(Int_t iTimeBin = 0; iTimeBin < kmaxTime; iTimeBin++) {
1001           Short_t q = allBins[iRow][iPad*kmaxTime + iTimeBin];
1002           if (q <= 0) continue;
1003           q *= 16;
1004           dig->SetDigitFast((Short_t)q,iTimeBin,iPad);
1005         }
1006       }
1007       fDigitsArray->StoreRow(iSector,iRow);
1008       Int_t ndig = dig->GetDigitSize();
1009
1010       AliDebug(10,
1011                Form("*** Sector, row, compressed digits %d %d %d ***\n",
1012                     iSector,iRow,ndig));
1013
1014       fDigitsArray->ClearRow(iSector,iRow);
1015
1016     } // end of the sector digitization
1017
1018     for (Int_t iRow = 0; iRow < nRows; iRow++)
1019       delete [] allBins[iRow];
1020
1021     delete [] allBins;
1022
1023   }
1024
1025   fLoader->WriteSDigits("OVERWRITE");
1026
1027   if(GetDigitsArray()) delete GetDigitsArray();
1028   SetDigitsArray(0x0);
1029
1030   return kTRUE;
1031 }
1032
1033 //______________________________________________________________________
1034 AliDigitizer* AliTPC::CreateDigitizer(AliRunDigitizer* manager) const
1035 {
1036   return new AliTPCDigitizer(manager);
1037 }
1038 //__
1039 void AliTPC::SDigits2Digits2(Int_t /*eventnumber*/)
1040 {
1041   //create digits from summable digits
1042   GenerNoise(500000); //create teble with noise
1043
1044   //conect tree with sSDigits
1045   TTree *t = fLoader->TreeS();
1046
1047   if (t == 0x0) {
1048     fLoader->LoadSDigits("READ");
1049     t = fLoader->TreeS();
1050     if (t == 0x0) {
1051       AliError("Can not get input TreeS");
1052       return;
1053     }
1054   }
1055
1056   if (fLoader->TreeD() == 0x0) fLoader->MakeTree("D");
1057
1058   AliSimDigits digarr, *dummy=&digarr;
1059   TBranch* sdb = t->GetBranch("Segment");
1060   if (sdb == 0x0) {
1061     AliError("Can not find branch with segments in TreeS.");
1062     return;
1063   }
1064
1065   sdb->SetAddress(&dummy);
1066
1067   Stat_t nentries = t->GetEntries();
1068
1069   // set zero suppression
1070
1071   fTPCParam->SetZeroSup(2);
1072
1073   // get zero suppression
1074
1075   Int_t zerosup = fTPCParam->GetZeroSup();
1076
1077   //make tree with digits
1078
1079   AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1080   arr->SetClass("AliSimDigits");
1081   arr->Setup(fTPCParam);
1082   arr->MakeTree(fLoader->TreeD());
1083
1084   AliTPCParam * par = fTPCParam;
1085
1086   //Loop over segments of the TPC
1087
1088   for (Int_t n=0; n<nentries; n++) {
1089     t->GetEvent(n);
1090     Int_t sec, row;
1091     if (!par->AdjustSectorRow(digarr.GetID(),sec,row)) {
1092       AliWarning(Form("Invalid segment ID ! %d",digarr.GetID()));
1093       continue;
1094     }
1095     if (!IsSectorActive(sec)) continue;
1096
1097     AliSimDigits * digrow =(AliSimDigits*) arr->CreateRow(sec,row);
1098     Int_t nrows = digrow->GetNRows();
1099     Int_t ncols = digrow->GetNCols();
1100
1101     digrow->ExpandBuffer();
1102     digarr.ExpandBuffer();
1103     digrow->ExpandTrackBuffer();
1104     digarr.ExpandTrackBuffer();
1105
1106
1107     Short_t * pamp0 = digarr.GetDigits();
1108     Int_t   * ptracks0 = digarr.GetTracks();
1109     Short_t * pamp1 = digrow->GetDigits();
1110     Int_t   * ptracks1 = digrow->GetTracks();
1111     Int_t  nelems =nrows*ncols;
1112     Int_t saturation = fTPCParam->GetADCSat() - 1;
1113     //use internal structure of the AliDigits - for speed reason
1114     //if you cahnge implementation
1115     //of the Alidigits - it must be rewriten -
1116     for (Int_t i= 0; i<nelems; i++){
1117       Float_t q = TMath::Nint(Float_t(*pamp0)/16.+GetNoise());
1118       if (q>zerosup){
1119         if (q>saturation) q=saturation;
1120         *pamp1=(Short_t)q;
1121
1122         ptracks1[0]=ptracks0[0];
1123         ptracks1[nelems]=ptracks0[nelems];
1124         ptracks1[2*nelems]=ptracks0[2*nelems];
1125       }
1126       pamp0++;
1127       pamp1++;
1128       ptracks0++;
1129       ptracks1++;
1130     }
1131
1132     arr->StoreRow(sec,row);
1133     arr->ClearRow(sec,row);
1134   }
1135
1136
1137   //write results
1138   fLoader->WriteDigits("OVERWRITE");
1139
1140   delete arr;
1141 }
1142 //__________________________________________________________________
1143 void AliTPC::SetDefaults(){
1144   //
1145   // setting the defaults
1146   //
1147
1148   // Set response functions
1149
1150   //
1151   AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::GetRunLoaderName());
1152   rl->CdGAFile();
1153   AliTPCParamSR *param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60");
1154   if(param){
1155     AliInfo("You are using 2 pad-length geom hits with 3 pad-lenght geom digits...");
1156     delete param;
1157     param = new AliTPCParamSR();
1158   }
1159   else {
1160     param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60_150x60");
1161   }
1162   if(!param){
1163     AliFatal("No TPC parameters found");
1164   }
1165
1166
1167   AliTPCPRF2D    * prfinner   = new AliTPCPRF2D;
1168   AliTPCPRF2D    * prfouter1   = new AliTPCPRF2D;
1169   AliTPCPRF2D    * prfouter2   = new AliTPCPRF2D;
1170   AliTPCRF1D     * rf    = new AliTPCRF1D(kTRUE);
1171   rf->SetGauss(param->GetZSigma(),param->GetZWidth(),1.);
1172   rf->SetOffset(3*param->GetZSigma());
1173   rf->Update();
1174
1175   TDirectory *savedir=gDirectory;
1176   TFile *f=TFile::Open("$ALICE_ROOT/TPC/AliTPCprf2d.root");
1177   if (!f->IsOpen())
1178     AliFatal("Can't open $ALICE_ROOT/TPC/AliTPCprf2d.root !");
1179
1180   TString s;
1181   prfinner->Read("prf_07504_Gati_056068_d02");
1182   //PH Set different names
1183   s=prfinner->GetGRF()->GetName();
1184   s+="in";
1185   prfinner->GetGRF()->SetName(s.Data());
1186
1187   prfouter1->Read("prf_10006_Gati_047051_d03");
1188   s=prfouter1->GetGRF()->GetName();
1189   s+="out1";
1190   prfouter1->GetGRF()->SetName(s.Data());
1191
1192   prfouter2->Read("prf_15006_Gati_047051_d03");
1193   s=prfouter2->GetGRF()->GetName();
1194   s+="out2";
1195   prfouter2->GetGRF()->SetName(s.Data());
1196
1197   f->Close();
1198   savedir->cd();
1199
1200   param->SetInnerPRF(prfinner);
1201   param->SetOuter1PRF(prfouter1);
1202   param->SetOuter2PRF(prfouter2);
1203   param->SetTimeRF(rf);
1204
1205   // set fTPCParam
1206
1207   SetParam(param);
1208
1209
1210   fDefaults = 1;
1211
1212 }
1213 //__________________________________________________________________
1214 void AliTPC::Hits2Digits()
1215 {
1216   //
1217   // creates digits from hits
1218   //
1219   if (!fTPCParam->IsGeoRead()){
1220     //
1221     // read transformation matrices for gGeoManager
1222     //
1223     fTPCParam->ReadGeoMatrices();
1224   }
1225
1226   fLoader->LoadHits("read");
1227   fLoader->LoadDigits("recreate");
1228   AliRunLoader* runLoader = fLoader->GetRunLoader();
1229
1230   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
1231     //PH    runLoader->GetEvent(iEvent);
1232     SetActiveSectors();
1233     Hits2Digits(iEvent);
1234   }
1235
1236   fLoader->UnloadHits();
1237   fLoader->UnloadDigits();
1238 }
1239 //__________________________________________________________________
1240 void AliTPC::Hits2Digits(Int_t eventnumber)
1241 {
1242  //----------------------------------------------------
1243  // Loop over all sectors for a single event
1244  //----------------------------------------------------
1245   AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::GetRunLoaderName());
1246   rl->GetEvent(eventnumber);
1247   if (fLoader->TreeH() == 0x0) {
1248     if(fLoader->LoadHits()) {
1249       AliError("Can not load hits.");
1250     }
1251   }
1252   SetTreeAddress();
1253
1254   if (fLoader->TreeD() == 0x0 ) {
1255     fLoader->MakeTree("D");
1256     if (fLoader->TreeD() == 0x0 ) {
1257       AliError("Can not get TreeD");
1258       return;
1259     }
1260   }
1261
1262   if(fDefaults == 0) SetDefaults();  // check if the parameters are set
1263   GenerNoise(500000); //create teble with noise
1264
1265   //setup TPCDigitsArray
1266
1267   if(GetDigitsArray()) delete GetDigitsArray();
1268
1269   AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1270   arr->SetClass("AliSimDigits");
1271   arr->Setup(fTPCParam);
1272
1273   arr->MakeTree(fLoader->TreeD());
1274   SetDigitsArray(arr);
1275
1276   fDigitsSwitch=0; // standard digits
1277
1278   for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++)
1279     if (IsSectorActive(isec)) {
1280       AliDebug(1,Form("Hits2Digits","Sector %d is active.",isec));
1281       Hits2DigitsSector(isec);
1282     }
1283     else {
1284       AliDebug(1,Form("Hits2Digits","Sector %d is NOT active.",isec));
1285     }
1286
1287   fLoader->WriteDigits("OVERWRITE");
1288
1289 //this line prevents the crash in the similar one
1290 //on the beginning of this method
1291 //destructor attempts to reset the tree, which is deleted by the loader
1292 //need to be redesign
1293   if(GetDigitsArray()) delete GetDigitsArray();
1294   SetDigitsArray(0x0);
1295
1296 }
1297
1298 //__________________________________________________________________
1299 void AliTPC::Hits2SDigits2(Int_t eventnumber)
1300 {
1301
1302   //-----------------------------------------------------------
1303   //   summable digits - 16 bit "ADC", no noise, no saturation
1304   //-----------------------------------------------------------
1305
1306   //----------------------------------------------------
1307   // Loop over all sectors for a single event
1308   //----------------------------------------------------
1309
1310   AliRunLoader* rl = fLoader->GetRunLoader();
1311
1312   rl->GetEvent(eventnumber);
1313   if (fLoader->TreeH() == 0x0) {
1314     if(fLoader->LoadHits()) {
1315       AliError("Can not load hits.");
1316       return;
1317     }
1318   }
1319   SetTreeAddress();
1320
1321
1322   if (fLoader->TreeS() == 0x0 ) {
1323     fLoader->MakeTree("S");
1324   }
1325
1326   if(fDefaults == 0) SetDefaults();
1327
1328   GenerNoise(500000); //create table with noise
1329   //setup TPCDigitsArray
1330
1331   if(GetDigitsArray()) delete GetDigitsArray();
1332
1333
1334   AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1335   arr->SetClass("AliSimDigits");
1336   arr->Setup(fTPCParam);
1337   arr->MakeTree(fLoader->TreeS());
1338
1339   SetDigitsArray(arr);
1340
1341   fDigitsSwitch=1; // summable digits
1342
1343     // set zero suppression to "0"
1344
1345   fTPCParam->SetZeroSup(0);
1346
1347   for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++)
1348     if (IsSectorActive(isec)) {
1349       Hits2DigitsSector(isec);
1350     }
1351
1352   fLoader->WriteSDigits("OVERWRITE");
1353
1354 //this line prevents the crash in the similar one
1355 //on the beginning of this method
1356 //destructor attempts to reset the tree, which is deleted by the loader
1357 //need to be redesign
1358   if(GetDigitsArray()) delete GetDigitsArray();
1359   SetDigitsArray(0x0);
1360 }
1361 //__________________________________________________________________
1362
1363 void AliTPC::Hits2SDigits()
1364 {
1365
1366   //-----------------------------------------------------------
1367   //   summable digits - 16 bit "ADC", no noise, no saturation
1368   //-----------------------------------------------------------
1369
1370   if (!fTPCParam->IsGeoRead()){
1371     //
1372     // read transformation matrices for gGeoManager
1373     //
1374     fTPCParam->ReadGeoMatrices();
1375   }
1376
1377   fLoader->LoadHits("read");
1378   fLoader->LoadSDigits("recreate");
1379   AliRunLoader* runLoader = fLoader->GetRunLoader();
1380
1381   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
1382     runLoader->GetEvent(iEvent);
1383     SetTreeAddress();
1384     SetActiveSectors();
1385     Hits2SDigits2(iEvent);
1386   }
1387
1388   fLoader->UnloadHits();
1389   fLoader->UnloadSDigits();
1390 }
1391 //_____________________________________________________________________________
1392
1393 void AliTPC::Hits2DigitsSector(Int_t isec)
1394 {
1395   //-------------------------------------------------------------------
1396   // TPC conversion from hits to digits.
1397   //-------------------------------------------------------------------
1398
1399   //-----------------------------------------------------------------
1400   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1401   //-----------------------------------------------------------------
1402
1403   //-------------------------------------------------------
1404   //  Get the access to the track hits
1405   //-------------------------------------------------------
1406
1407   // check if the parameters are set - important if one calls this method
1408   // directly, not from the Hits2Digits
1409
1410   if(fDefaults == 0) SetDefaults();
1411
1412   TTree *tH = TreeH(); // pointer to the hits tree
1413   if (tH == 0x0) {
1414     AliFatal("Can not find TreeH in folder");
1415     return;
1416   }
1417
1418   Stat_t ntracks = tH->GetEntries();
1419
1420   if( ntracks > 0){
1421
1422   //-------------------------------------------
1423   //  Only if there are any tracks...
1424   //-------------------------------------------
1425
1426     TObjArray **row;
1427
1428     Int_t nrows =fTPCParam->GetNRow(isec);
1429
1430     row= new TObjArray* [nrows+2]; // 2 extra rows for cross talk
1431
1432     MakeSector(isec,nrows,tH,ntracks,row);
1433
1434     //--------------------------------------------------------
1435     //   Digitize this sector, row by row
1436     //   row[i] is the pointer to the TObjArray of TVectors,
1437     //   each one containing electrons accepted on this
1438     //   row, assigned into tracks
1439     //--------------------------------------------------------
1440
1441     Int_t i;
1442
1443     if (fDigitsArray->GetTree()==0) {
1444       AliFatal("Tree not set in fDigitsArray");
1445     }
1446
1447     for (i=0;i<nrows;i++){
1448
1449       AliDigits * dig = fDigitsArray->CreateRow(isec,i);
1450
1451       DigitizeRow(i,isec,row);
1452
1453       fDigitsArray->StoreRow(isec,i);
1454
1455       Int_t ndig = dig->GetDigitSize();
1456
1457       AliDebug(10,
1458                Form("*** Sector, row, compressed digits %d %d %d ***\n",
1459                     isec,i,ndig));
1460
1461       fDigitsArray->ClearRow(isec,i);
1462
1463
1464     } // end of the sector digitization
1465
1466     for(i=0;i<nrows+2;i++){
1467       row[i]->Delete();
1468       delete row[i];
1469     }
1470
1471     delete [] row; // delete the array of pointers to TObjArray-s
1472
1473   } // ntracks >0
1474
1475 } // end of Hits2DigitsSector
1476
1477
1478 //_____________________________________________________________________________
1479 void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows)
1480 {
1481   //-----------------------------------------------------------
1482   // Single row digitization, coupling from the neighbouring
1483   // rows taken into account
1484   //-----------------------------------------------------------
1485
1486   //-----------------------------------------------------------------
1487   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1488   // Modified: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
1489   //-----------------------------------------------------------------
1490
1491   Float_t zerosup = fTPCParam->GetZeroSup();
1492
1493   fCurrentIndex[1]= isec;
1494
1495
1496   Int_t nofPads = fTPCParam->GetNPads(isec,irow);
1497   Int_t nofTbins = fTPCParam->GetMaxTBin();
1498   Int_t indexRange[4];
1499   //
1500   //  Integrated signal for this row
1501   //  and a single track signal
1502   //
1503
1504   TMatrixF *m1 = new TMatrixF(0,nofPads,0,nofTbins); // integrated
1505   TMatrixF *m2 = new TMatrixF(0,nofPads,0,nofTbins); // single
1506   //
1507   TMatrixF &total  = *m1;
1508
1509   //  Array of pointers to the label-signal list
1510
1511   Int_t nofDigits = nofPads*nofTbins; // number of digits for this row
1512   Float_t  **pList = new Float_t* [nofDigits];
1513
1514   Int_t lp;
1515   Int_t i1;
1516   for(lp=0;lp<nofDigits;lp++)pList[lp]=0; // set all pointers to NULL
1517   //
1518   //calculate signal
1519   //
1520   Int_t row1=irow;
1521   Int_t row2=irow+2;
1522   for (Int_t row= row1;row<=row2;row++){
1523     Int_t nTracks= rows[row]->GetEntries();
1524     for (i1=0;i1<nTracks;i1++){
1525       fCurrentIndex[2]= row;
1526       fCurrentIndex[3]=irow+1;
1527       if (row==irow+1){
1528         m2->Zero();  // clear single track signal matrix
1529         Float_t trackLabel = GetSignal(rows[row],i1,m2,m1,indexRange);
1530         GetList(trackLabel,nofPads,m2,indexRange,pList);
1531       }
1532       else   GetSignal(rows[row],i1,0,m1,indexRange);
1533     }
1534   }
1535
1536   Int_t tracks[3];
1537
1538   AliDigits *dig = fDigitsArray->GetRow(isec,irow);
1539   Int_t gi=-1;
1540   Float_t fzerosup = zerosup+0.5;
1541   for(Int_t it=0;it<nofTbins;it++){
1542     for(Int_t ip=0;ip<nofPads;ip++){
1543       gi++;
1544       Float_t q=total(ip,it);
1545       if(fDigitsSwitch == 0){
1546         q+=GetNoise();
1547         if(q <=fzerosup) continue; // do not fill zeros
1548         q = TMath::Nint(q);
1549         if(q >= fTPCParam->GetADCSat()) q = fTPCParam->GetADCSat() - 1;  // saturation
1550
1551       }
1552
1553       else {
1554         if(q <= 0.) continue; // do not fill zeros
1555         if(q>2000.) q=2000.;
1556         q *= 16.;
1557         q = TMath::Nint(q);
1558       }
1559
1560       //
1561       //  "real" signal or electronic noise (list = -1)?
1562       //
1563
1564       for(Int_t j1=0;j1<3;j1++){
1565         tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -2;
1566       }
1567
1568 //Begin_Html
1569 /*
1570   <A NAME="AliDigits"></A>
1571   using of AliDigits object
1572 */
1573 //End_Html
1574       dig->SetDigitFast((Short_t)q,it,ip);
1575       if (fDigitsArray->IsSimulated()) {
1576         ((AliSimDigits*)dig)->SetTrackIDFast(tracks[0],it,ip,0);
1577         ((AliSimDigits*)dig)->SetTrackIDFast(tracks[1],it,ip,1);
1578         ((AliSimDigits*)dig)->SetTrackIDFast(tracks[2],it,ip,2);
1579       }
1580
1581     } // end of loop over time buckets
1582   }  // end of lop over pads
1583
1584   //
1585   //  This row has been digitized, delete nonused stuff
1586   //
1587
1588   for(lp=0;lp<nofDigits;lp++){
1589     if(pList[lp]) delete [] pList[lp];
1590   }
1591
1592   delete [] pList;
1593
1594   delete m1;
1595   delete m2;
1596
1597 } // end of DigitizeRow
1598
1599 //_____________________________________________________________________________
1600
1601 Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr,
1602              TMatrixF *m1, TMatrixF *m2,Int_t *indexRange)
1603 {
1604
1605   //---------------------------------------------------------------
1606   //  Calculates 2-D signal (pad,time) for a single track,
1607   //  returns a pointer to the signal matrix and the track label
1608   //  No digitization is performed at this level!!!
1609   //---------------------------------------------------------------
1610
1611   //-----------------------------------------------------------------
1612   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1613   // Modified: Marian Ivanov
1614   //-----------------------------------------------------------------
1615
1616   TVector *tv;
1617
1618   tv = (TVector*)p1->At(ntr); // pointer to a track
1619   TVector &v = *tv;
1620
1621   Float_t label = v(0);
1622   Int_t centralPad = (fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3]-1))/2;
1623
1624   Int_t nElectrons = (tv->GetNrows()-1)/5;
1625   indexRange[0]=9999; // min pad
1626   indexRange[1]=-1; // max pad
1627   indexRange[2]=9999; //min time
1628   indexRange[3]=-1; // max time
1629
1630   TMatrixF &signal = *m1;
1631   TMatrixF &total = *m2;
1632   //
1633   //  Loop over all electrons
1634   //
1635   for(Int_t nel=0; nel<nElectrons; nel++){
1636     Int_t idx=nel*5;
1637     Float_t aval =  v(idx+4);
1638     Float_t eltoadcfac=aval*fTPCParam->GetTotalNormFac();
1639     Float_t xyz[4]={v(idx+1),v(idx+2),v(idx+3),v(idx+5)};
1640     Int_t n = ((AliTPCParamSR*)fTPCParam)->CalcResponseFast(xyz,fCurrentIndex,fCurrentIndex[3]);
1641
1642     Int_t *index = fTPCParam->GetResBin(0);
1643     Float_t *weight = & (fTPCParam->GetResWeight(0));
1644
1645     if (n>0) for (Int_t i =0; i<n; i++){
1646       Int_t pad=index[1]+centralPad;  //in digit coordinates central pad has coordinate 0
1647
1648       if (pad>=0){
1649         Int_t time=index[2];
1650         Float_t qweight = *(weight)*eltoadcfac;
1651
1652         if (m1!=0) signal(pad,time)+=qweight;
1653         total(pad,time)+=qweight;
1654         if (indexRange[0]>pad) indexRange[0]=pad;
1655         if (indexRange[1]<pad) indexRange[1]=pad;
1656         if (indexRange[2]>time) indexRange[2]=time;
1657         if (indexRange[3]<time) indexRange[3]=time;
1658
1659         index+=3;
1660         weight++;
1661
1662       }
1663     }
1664   } // end of loop over electrons
1665
1666   return label; // returns track label when finished
1667 }
1668
1669 //_____________________________________________________________________________
1670 void AliTPC::GetList(Float_t label,Int_t np,TMatrixF *m,
1671                      Int_t *indexRange, Float_t **pList)
1672 {
1673   //----------------------------------------------------------------------
1674   //  Updates the list of tracks contributing to digits for a given row
1675   //----------------------------------------------------------------------
1676
1677   //-----------------------------------------------------------------
1678   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1679   //-----------------------------------------------------------------
1680
1681   TMatrixF &signal = *m;
1682
1683   // lop over nonzero digits
1684
1685   for(Int_t it=indexRange[2];it<indexRange[3]+1;it++){
1686     for(Int_t ip=indexRange[0];ip<indexRange[1]+1;ip++){
1687
1688
1689       // accept only the contribution larger than 500 electrons (1/2 s_noise)
1690
1691       if(signal(ip,it)<0.5) continue;
1692
1693       Int_t globalIndex = it*np+ip; // globalIndex starts from 0!
1694
1695       if(!pList[globalIndex]){
1696
1697         //
1698         // Create new list (6 elements - 3 signals and 3 labels),
1699         //
1700
1701         pList[globalIndex] = new Float_t [6];
1702
1703         // set list to -1
1704
1705         *pList[globalIndex] = -1.;
1706         *(pList[globalIndex]+1) = -1.;
1707         *(pList[globalIndex]+2) = -1.;
1708         *(pList[globalIndex]+3) = -1.;
1709         *(pList[globalIndex]+4) = -1.;
1710         *(pList[globalIndex]+5) = -1.;
1711
1712         *pList[globalIndex] = label;
1713         *(pList[globalIndex]+3) = signal(ip,it);
1714       }
1715       else {
1716
1717         // check the signal magnitude
1718
1719         Float_t highest = *(pList[globalIndex]+3);
1720         Float_t middle = *(pList[globalIndex]+4);
1721         Float_t lowest = *(pList[globalIndex]+5);
1722
1723         //
1724         //  compare the new signal with already existing list
1725         //
1726
1727         if(signal(ip,it)<lowest) continue; // neglect this track
1728
1729         //
1730
1731         if (signal(ip,it)>highest){
1732           *(pList[globalIndex]+5) = middle;
1733           *(pList[globalIndex]+4) = highest;
1734           *(pList[globalIndex]+3) = signal(ip,it);
1735
1736           *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
1737           *(pList[globalIndex]+1) = *pList[globalIndex];
1738           *pList[globalIndex] = label;
1739         }
1740         else if (signal(ip,it)>middle){
1741           *(pList[globalIndex]+5) = middle;
1742           *(pList[globalIndex]+4) = signal(ip,it);
1743
1744           *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
1745           *(pList[globalIndex]+1) = label;
1746         }
1747         else{
1748           *(pList[globalIndex]+5) = signal(ip,it);
1749           *(pList[globalIndex]+2) = label;
1750         }
1751       }
1752
1753     } // end of loop over pads
1754   } // end of loop over time bins
1755
1756 }//end of GetList
1757 //___________________________________________________________________
1758 void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
1759                         Stat_t ntracks,TObjArray **row)
1760 {
1761
1762   //-----------------------------------------------------------------
1763   // Prepares the sector digitization, creates the vectors of
1764   // tracks for each row of this sector. The track vector
1765   // contains the track label and the position of electrons.
1766   //-----------------------------------------------------------------
1767
1768   //
1769   // The trasport of the electrons through TPC drift volume
1770   //    Drift (drift velocity + velocity map(not yet implemented)))
1771   //    Application of the random processes (diffusion, gas gain)
1772   //    Systematic effects (ExB effect in drift volume + ROCs)
1773   //
1774   // Algorithm:
1775   // Loop over primary electrons:
1776   //    Creation of the secondary electrons
1777   //    Loop over electrons (primary+ secondaries)
1778   //        Global coordinate frame:
1779   //          1. Skip electrons if attached
1780   //          2. ExB effect in drift volume
1781   //             a.) Simulation   calib->GetExB()->CorrectInverse(dxyz0,dxyz1);
1782   //             b.) Reconstruction -  calib->GetExB()->CorrectInverse(dxyz0,dxyz1);
1783   //          3. Generation of gas gain (Random - Exponential distribution)
1784   //          4. TransportElectron function (diffusion)
1785   //
1786   //        5. Conversion to the local coordinate frame  pad-row, pad, timebin
1787   //        6. Apply Time0 shift - AliTPCCalPad class
1788   //            a.) Plus sign in simulation
1789   //            b.) Minus sign in reconstruction
1790   // end of loop
1791   //
1792   //-----------------------------------------------------------------
1793   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1794   // Origin: Marian Ivanov,  marian.ivanov@cern.ch
1795   //-----------------------------------------------------------------
1796   AliTPCcalibDB* const calib=AliTPCcalibDB::Instance();
1797
1798   Float_t gasgain = fTPCParam->GetGasGain();
1799   Int_t i;
1800   Float_t xyz[5];
1801
1802   AliTPChit *tpcHit; // pointer to a sigle TPC hit
1803   //MI change
1804   TBranch * branch=0;
1805   if (fHitType>1) branch = TH->GetBranch("TPC2");
1806   else branch = TH->GetBranch("TPC");
1807
1808
1809   //----------------------------------------------
1810   // Create TObjArray-s, one for each row,
1811   // each TObjArray will store the TVectors
1812   // of electrons, one TVectors per each track.
1813   //----------------------------------------------
1814
1815   Int_t *nofElectrons = new Int_t [nrows+2]; // electron counter for each row
1816   TVector **tracks = new TVector* [nrows+2]; //pointers to the track vectors
1817
1818   for(i=0; i<nrows+2; i++){
1819     row[i] = new TObjArray;
1820     nofElectrons[i]=0;
1821     tracks[i]=0;
1822   }
1823
1824
1825
1826   //--------------------------------------------------------------------
1827   //  Loop over tracks, the "track" contains the full history
1828   //--------------------------------------------------------------------
1829
1830   Int_t previousTrack,currentTrack;
1831   previousTrack = -1; // nothing to store so far!
1832
1833   for(Int_t track=0;track<ntracks;track++){
1834     Bool_t isInSector=kTRUE;
1835     ResetHits();
1836     isInSector = TrackInVolume(isec,track);
1837     if (!isInSector) continue;
1838     //MI change
1839     branch->GetEntry(track); // get next track
1840
1841     //M.I. changes
1842
1843     tpcHit = (AliTPChit*)FirstHit(-1);
1844
1845     //--------------------------------------------------------------
1846     //  Loop over hits
1847     //--------------------------------------------------------------
1848
1849
1850     while(tpcHit){
1851
1852       Int_t sector=tpcHit->fSector; // sector number
1853       if(sector != isec){
1854         tpcHit = (AliTPChit*) NextHit();
1855         continue;
1856       }
1857
1858       // Remove hits which arrive before the TPC opening gate signal
1859       if(((fTPCParam->GetZLength(isec)-TMath::Abs(tpcHit->Z()))
1860           /fTPCParam->GetDriftV()+tpcHit->Time())<fTPCParam->GetGateDelay()) {
1861         tpcHit = (AliTPChit*) NextHit();
1862         continue;
1863       }
1864
1865       currentTrack = tpcHit->Track(); // track number
1866
1867       if(currentTrack != previousTrack){
1868
1869         // store already filled fTrack
1870
1871         for(i=0;i<nrows+2;i++){
1872           if(previousTrack != -1){
1873             if(nofElectrons[i]>0){
1874               TVector &v = *tracks[i];
1875               v(0) = previousTrack;
1876               tracks[i]->ResizeTo(5*nofElectrons[i]+1); // shrink if necessary
1877               row[i]->Add(tracks[i]);
1878             }
1879             else {
1880               delete tracks[i]; // delete empty TVector
1881               tracks[i]=0;
1882             }
1883           }
1884
1885           nofElectrons[i]=0;
1886           tracks[i] = new TVector(601); // TVectors for the next fTrack
1887
1888         } // end of loop over rows
1889
1890         previousTrack=currentTrack; // update track label
1891       }
1892
1893       Int_t qI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons)
1894
1895       //---------------------------------------------------
1896       //  Calculate the electron attachment probability
1897       //---------------------------------------------------
1898
1899
1900       Float_t time = 1.e6*(fTPCParam->GetZLength(isec)-TMath::Abs(tpcHit->Z()))
1901         /fTPCParam->GetDriftV();
1902       // in microseconds!
1903       Float_t attProb = fTPCParam->GetAttCoef()*
1904         fTPCParam->GetOxyCont()*time; //  fraction!
1905
1906       //-----------------------------------------------
1907       //  Loop over electrons
1908       //-----------------------------------------------
1909       Int_t index[3];
1910       index[1]=isec;
1911       for(Int_t nel=0;nel<qI;nel++){
1912         // skip if electron lost due to the attachment
1913         if((gRandom->Rndm(0)) < attProb) continue; // electron lost!
1914
1915         //
1916         // ExB effect
1917         //
1918         Double_t dxyz0[3],dxyz1[3];
1919         dxyz0[0]=tpcHit->X();
1920         dxyz0[1]=tpcHit->Y();
1921         dxyz0[2]=tpcHit->Z();
1922         if (calib->GetExB()){
1923           calib->GetExB()->CorrectInverse(dxyz0,dxyz1);
1924         }else{
1925           AliError("Not valid ExB calibration");
1926           dxyz1[0]=tpcHit->X();
1927           dxyz1[1]=tpcHit->Y();
1928           dxyz1[2]=tpcHit->Z();
1929         }
1930         xyz[0]=dxyz1[0];
1931         xyz[1]=dxyz1[1];
1932         xyz[2]=dxyz1[2];
1933         //
1934         //
1935         //
1936         // protection for the nonphysical avalanche size (10**6 maximum)
1937         //
1938         Double_t rn=TMath::Max(gRandom->Rndm(0),1.93e-22);
1939         xyz[3]= (Float_t) (-gasgain*TMath::Log(rn));
1940         index[0]=1;
1941
1942         TransportElectron(xyz,index);
1943         Int_t rowNumber;
1944         fTPCParam->GetPadRow(xyz,index);
1945         //
1946         // Add Time0 correction due unisochronity
1947         // xyz[0] - pad row coordinate
1948         // xyz[1] - pad coordinate
1949         // xyz[2] - is in now time bin coordinate system
1950         Float_t correction =0;
1951         if (calib->GetPadTime0()){
1952           if (!calib->GetPadTime0()->GetCalROC(isec)) continue;
1953           correction = calib->GetPadTime0()->GetCalROC(isec)->GetValue(TMath::Nint(xyz[0]),TMath::Nint(xyz[1]));
1954         }
1955         xyz[2]+=correction;
1956         //
1957         // Electron track time (for pileup simulation)
1958         xyz[4] = tpcHit->Time()/fTPCParam->GetTSample();
1959         // row 0 - cross talk from the innermost row
1960         // row fNRow+1 cross talk from the outermost row
1961         rowNumber = index[2]+1;
1962         //transform position to local digit coordinates
1963         //relative to nearest pad row
1964         if ((rowNumber<0)||rowNumber>fTPCParam->GetNRow(isec)+1) continue;
1965         Float_t x1,y1;
1966         if (isec <fTPCParam->GetNInnerSector()) {
1967           x1 = xyz[1]*fTPCParam->GetInnerPadPitchWidth();
1968           y1 = fTPCParam->GetYInner(rowNumber);
1969         }
1970         else{
1971           x1=xyz[1]*fTPCParam->GetOuterPadPitchWidth();
1972           y1 = fTPCParam->GetYOuter(rowNumber);
1973         }
1974         // gain inefficiency at the wires edges - linear
1975         x1=TMath::Abs(x1);
1976         y1-=1.;
1977         if(x1>y1) xyz[3]*=TMath::Max(1.e-6,(y1-x1+1.));
1978
1979         nofElectrons[rowNumber]++;
1980         //----------------------------------
1981         // Expand vector if necessary
1982         //----------------------------------
1983         if(nofElectrons[rowNumber]>120){
1984           Int_t range = tracks[rowNumber]->GetNrows();
1985           if((nofElectrons[rowNumber])>(range-1)/5){
1986
1987             tracks[rowNumber]->ResizeTo(range+500); // Add 100 electrons
1988           }
1989         }
1990
1991         TVector &v = *tracks[rowNumber];
1992         Int_t idx = 5*nofElectrons[rowNumber]-4;
1993         Real_t * position = &(((TVector&)v)(idx)); //make code faster
1994         memcpy(position,xyz,5*sizeof(Float_t));
1995
1996       } // end of loop over electrons
1997
1998       tpcHit = (AliTPChit*)NextHit();
1999
2000     } // end of loop over hits
2001   } // end of loop over tracks
2002
2003     //
2004     //   store remaining track (the last one) if not empty
2005     //
2006
2007   for(i=0;i<nrows+2;i++){
2008     if(nofElectrons[i]>0){
2009       TVector &v = *tracks[i];
2010       v(0) = previousTrack;
2011       tracks[i]->ResizeTo(5*nofElectrons[i]+1); // shrink if necessary
2012       row[i]->Add(tracks[i]);
2013     }
2014     else{
2015       delete tracks[i];
2016       tracks[i]=0;
2017     }
2018   }
2019
2020   delete [] tracks;
2021   delete [] nofElectrons;
2022
2023 } // end of MakeSector
2024
2025
2026 //_____________________________________________________________________________
2027 void AliTPC::Init()
2028 {
2029   //
2030   // Initialise TPC detector after definition of geometry
2031   //
2032   AliDebug(1,"*********************************************");
2033 }
2034
2035 //_____________________________________________________________________________
2036 void AliTPC::MakeBranch(Option_t* option)
2037 {
2038   //
2039   // Create Tree branches for the TPC.
2040   //
2041   AliDebug(1,"");
2042   Int_t buffersize = 4000;
2043   char branchname[10];
2044   sprintf(branchname,"%s",GetName());
2045
2046   const char *h = strstr(option,"H");
2047
2048   if ( h && (fHitType<=1) && (fHits == 0x0)) fHits = new TClonesArray("AliTPChit", 176);//skowron 20.06.03
2049
2050   AliDetector::MakeBranch(option);
2051
2052   const char *d = strstr(option,"D");
2053
2054   if (fDigits   && fLoader->TreeD() && d) {
2055     MakeBranchInTree(gAlice->TreeD(), branchname, &fDigits, buffersize, 0);
2056   }
2057
2058   if (fHitType>1) MakeBranch2(option,0); // MI change 14.09.2000
2059 }
2060
2061 //_____________________________________________________________________________
2062 void AliTPC::ResetDigits()
2063 {
2064   //
2065   // Reset number of digits and the digits array for this detector
2066   //
2067   fNdigits   = 0;
2068   if (fDigits)   fDigits->Clear();
2069 }
2070
2071
2072
2073 //_____________________________________________________________________________
2074 void AliTPC::SetSens(Int_t sens)
2075 {
2076
2077   //-------------------------------------------------------------
2078   // Activates/deactivates the sensitive strips at the center of
2079   // the pad row -- this is for the space-point resolution calculations
2080   //-------------------------------------------------------------
2081
2082   //-----------------------------------------------------------------
2083   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2084   //-----------------------------------------------------------------
2085
2086   fSens = sens;
2087 }
2088
2089
2090 void AliTPC::SetSide(Float_t side=0.)
2091 {
2092   // choice of the TPC side
2093
2094   fSide = side;
2095
2096 }
2097 //_____________________________________________________________________________
2098
2099 void AliTPC::TransportElectron(Float_t *xyz, Int_t *index)
2100 {
2101   //
2102   // electron transport taking into account:
2103   // 1. diffusion,
2104   // 2.ExB at the wires
2105   // 3. nonisochronity
2106   //
2107   // xyz and index must be already transformed to system 1
2108   //
2109
2110   fTPCParam->Transform1to2(xyz,index);
2111
2112   //add diffusion
2113   Float_t driftl=xyz[2];
2114   if(driftl<0.01) driftl=0.01;
2115   driftl=TMath::Sqrt(driftl);
2116   Float_t sigT = driftl*(fTPCParam->GetDiffT());
2117   Float_t sigL = driftl*(fTPCParam->GetDiffL());
2118   xyz[0]=gRandom->Gaus(xyz[0],sigT);
2119   xyz[1]=gRandom->Gaus(xyz[1],sigT);
2120   xyz[2]=gRandom->Gaus(xyz[2],sigL);
2121
2122   // ExB
2123
2124   if (fTPCParam->GetMWPCReadout()==kTRUE){
2125     Float_t dx = fTPCParam->Transform2to2NearestWire(xyz,index);
2126     xyz[1]+=dx*(fTPCParam->GetOmegaTau());
2127   }
2128   //add nonisochronity (not implemented yet)
2129
2130
2131 }
2132
2133 ClassImp(AliTPChit)
2134   //______________________________________________________________________
2135   AliTPChit::AliTPChit()
2136             :AliHit(),
2137              fSector(0),
2138              fPadRow(0),
2139              fQ(0),
2140              fTime(0)
2141 {
2142   //
2143   // default
2144   //
2145
2146 }
2147 //_____________________________________________________________________________
2148 AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits)
2149           :AliHit(shunt,track),
2150              fSector(0),
2151              fPadRow(0),
2152              fQ(0),
2153              fTime(0)
2154 {
2155   //
2156   // Creates a TPC hit object
2157   //
2158   fSector     = vol[0];
2159   fPadRow     = vol[1];
2160   fX          = hits[0];
2161   fY          = hits[1];
2162   fZ          = hits[2];
2163   fQ          = hits[3];
2164   fTime       = hits[4];
2165 }
2166
2167 //________________________________________________________________________
2168 // Additional code because of the AliTPCTrackHitsV2
2169
2170 void AliTPC::MakeBranch2(Option_t *option,const char */*file*/)
2171 {
2172   //
2173   // Create a new branch in the current Root Tree
2174   // The branch of fHits is automatically split
2175   // MI change 14.09.2000
2176   AliDebug(1,"");
2177   if (fHitType<2) return;
2178   char branchname[10];
2179   sprintf(branchname,"%s2",GetName());
2180   //
2181   // Get the pointer to the header
2182   const char *cH = strstr(option,"H");
2183   //
2184   if (fTrackHits   && TreeH() && cH && fHitType&4) {
2185     AliDebug(1,"Making branch for Type 4 Hits");
2186     TreeH()->Branch(branchname,"AliTPCTrackHitsV2",&fTrackHits,fBufferSize,99);
2187   }
2188
2189 //   if (fTrackHitsOld   && TreeH() && cH && fHitType&2) {
2190 //     AliDebug(1,"Making branch for Type 2 Hits");
2191 //     AliObjectBranch * branch = new AliObjectBranch(branchname,"AliTPCTrackHits",&fTrackHitsOld,
2192 //                                                    TreeH(),fBufferSize,99);
2193 //     TreeH()->GetListOfBranches()->Add(branch);
2194 //   }
2195 }
2196
2197 void AliTPC::SetTreeAddress()
2198 {
2199   //Sets tree address for hits
2200   if (fHitType<=1) {
2201     if (fHits == 0x0 ) fHits = new TClonesArray("AliTPChit", 176);//skowron 20.06.03
2202     AliDetector::SetTreeAddress();
2203   }
2204   if (fHitType>1) SetTreeAddress2();
2205 }
2206
2207 void AliTPC::SetTreeAddress2()
2208 {
2209   //
2210   // Set branch address for the TrackHits Tree
2211   //
2212   AliDebug(1,"");
2213
2214   TBranch *branch;
2215   char branchname[20];
2216   sprintf(branchname,"%s2",GetName());
2217   //
2218   // Branch address for hit tree
2219   TTree *treeH = TreeH();
2220   if ((treeH)&&(fHitType&4)) {
2221     branch = treeH->GetBranch(branchname);
2222     if (branch) {
2223       branch->SetAddress(&fTrackHits);
2224       AliDebug(1,"fHitType&4 Setting");
2225     }
2226     else
2227       AliDebug(1,"fHitType&4 Failed (can not find branch)");
2228
2229   }
2230  //  if ((treeH)&&(fHitType&2)) {
2231 //     branch = treeH->GetBranch(branchname);
2232 //     if (branch) {
2233 //       branch->SetAddress(&fTrackHitsOld);
2234 //       AliDebug(1,"fHitType&2 Setting");
2235 //     }
2236 //     else
2237 //       AliDebug(1,"fHitType&2 Failed (can not find branch)");
2238 //   }
2239   //set address to TREETR
2240
2241   TTree *treeTR = TreeTR();
2242   if (treeTR && fTrackReferences) {
2243     branch = treeTR->GetBranch(GetName());
2244     if (branch) branch->SetAddress(&fTrackReferences);
2245   }
2246
2247 }
2248
2249 void AliTPC::FinishPrimary()
2250 {
2251   if (fTrackHits &&fHitType&4)      fTrackHits->FlushHitStack();
2252   //  if (fTrackHitsOld && fHitType&2)  fTrackHitsOld->FlushHitStack();
2253 }
2254
2255
2256 void AliTPC::AddHit2(Int_t track, Int_t *vol, Float_t *hits)
2257 {
2258   //
2259   // add hit to the list
2260   Int_t rtrack;
2261   if (fIshunt) {
2262     int primary = gAlice->GetMCApp()->GetPrimary(track);
2263     gAlice->GetMCApp()->Particle(primary)->SetBit(kKeepBit);
2264     rtrack=primary;
2265   } else {
2266     rtrack=track;
2267     gAlice->GetMCApp()->FlagTrack(track);
2268   }
2269   if (fTrackHits && fHitType&4)
2270     fTrackHits->AddHitKartez(vol[0],rtrack, hits[0],
2271                              hits[1],hits[2],(Int_t)hits[3],hits[4]);
2272  //  if (fTrackHitsOld &&fHitType&2 )
2273 //     fTrackHitsOld->AddHitKartez(vol[0],rtrack, hits[0],
2274 //                                 hits[1],hits[2],(Int_t)hits[3]);
2275
2276 }
2277
2278 void AliTPC::ResetHits()
2279 {
2280   if (fHitType&1) AliDetector::ResetHits();
2281   if (fHitType>1) ResetHits2();
2282 }
2283
2284 void AliTPC::ResetHits2()
2285 {
2286   //
2287   //reset hits
2288   if (fTrackHits && fHitType&4) fTrackHits->Clear();
2289   // if (fTrackHitsOld && fHitType&2) fTrackHitsOld->Clear();
2290
2291 }
2292
2293 AliHit* AliTPC::FirstHit(Int_t track)
2294 {
2295   if (fHitType>1) return FirstHit2(track);
2296   return AliDetector::FirstHit(track);
2297 }
2298 AliHit* AliTPC::NextHit()
2299 {
2300   //
2301   // gets next hit
2302   //
2303   if (fHitType>1) return NextHit2();
2304
2305   return AliDetector::NextHit();
2306 }
2307
2308 AliHit* AliTPC::FirstHit2(Int_t track)
2309 {
2310   //
2311   // Initialise the hit iterator
2312   // Return the address of the first hit for track
2313   // If track>=0 the track is read from disk
2314   // while if track<0 the first hit of the current
2315   // track is returned
2316   //
2317   if(track>=0) {
2318     gAlice->ResetHits();
2319     TreeH()->GetEvent(track);
2320   }
2321   //
2322   if (fTrackHits && fHitType&4) {
2323     fTrackHits->First();
2324     return fTrackHits->GetHit();
2325   }
2326  //  if (fTrackHitsOld && fHitType&2) {
2327 //     fTrackHitsOld->First();
2328 //     return fTrackHitsOld->GetHit();
2329 //   }
2330
2331   else return 0;
2332 }
2333
2334 AliHit* AliTPC::NextHit2()
2335 {
2336   //
2337   //Return the next hit for the current track
2338
2339
2340 //   if (fTrackHitsOld && fHitType&2) {
2341 //     fTrackHitsOld->Next();
2342 //     return fTrackHitsOld->GetHit();
2343 //   }
2344   if (fTrackHits) {
2345     fTrackHits->Next();
2346     return fTrackHits->GetHit();
2347   }
2348   else
2349     return 0;
2350 }
2351
2352 void AliTPC::LoadPoints(Int_t)
2353 {
2354   //
2355   Int_t a = 0;
2356
2357   if(fHitType==1) AliDetector::LoadPoints(a);
2358   else LoadPoints2(a);
2359 }
2360
2361
2362 void AliTPC::RemapTrackHitIDs(Int_t *map)
2363 {
2364   //
2365   // remapping
2366   //
2367   if (!fTrackHits) return;
2368
2369 //   if (fTrackHitsOld && fHitType&2){
2370 //     AliObjectArray * arr = fTrackHitsOld->fTrackHitsInfo;
2371 //     for (UInt_t i=0;i<arr->GetSize();i++){
2372 //       AliTrackHitsInfo * info = (AliTrackHitsInfo *)(arr->At(i));
2373 //       info->fTrackID = map[info->fTrackID];
2374 //     }
2375 //   }
2376 //  if (fTrackHitsOld && fHitType&4){
2377   if (fTrackHits && fHitType&4){
2378     TClonesArray * arr = fTrackHits->GetArray();;
2379     for (Int_t i=0;i<arr->GetEntriesFast();i++){
2380       AliTrackHitsParamV2 * info = (AliTrackHitsParamV2 *)(arr->At(i));
2381       info->SetTrackID(map[info->GetTrackID()]);
2382     }
2383   }
2384 }
2385
2386 Bool_t   AliTPC::TrackInVolume(Int_t id,Int_t track)
2387 {
2388   //return bool information - is track in given volume
2389   //load only part of the track information
2390   //return true if current track is in volume
2391   //
2392   //  return kTRUE;
2393  //  if (fTrackHitsOld && fHitType&2) {
2394 //     TBranch * br = TreeH()->GetBranch("fTrackHitsInfo");
2395 //     br->GetEvent(track);
2396 //     AliObjectArray * ar = fTrackHitsOld->fTrackHitsInfo;
2397 //     for (UInt_t j=0;j<ar->GetSize();j++){
2398 //       if (  ((AliTrackHitsInfo*)ar->At(j))->fVolumeID==id) return kTRUE;
2399 //     }
2400 //   }
2401
2402   if (fTrackHits && fHitType&4) {
2403     TBranch * br1 = TreeH()->GetBranch("fVolumes");
2404     TBranch * br2 = TreeH()->GetBranch("fNVolumes");
2405     br2->GetEvent(track);
2406     br1->GetEvent(track);
2407     Int_t *volumes = fTrackHits->GetVolumes();
2408     Int_t nvolumes = fTrackHits->GetNVolumes();
2409     if (!volumes && nvolumes>0) {
2410       AliWarning(Form("Problematic track\t%d\t%d",track,nvolumes));
2411       return kFALSE;
2412     }
2413     for (Int_t j=0;j<nvolumes; j++)
2414       if (volumes[j]==id) return kTRUE;
2415   }
2416
2417   if (fHitType&1) {
2418     TBranch * br = TreeH()->GetBranch("fSector");
2419     br->GetEvent(track);
2420     for (Int_t j=0;j<fHits->GetEntriesFast();j++){
2421       if (  ((AliTPChit*)fHits->At(j))->fSector==id) return kTRUE;
2422     }
2423   }
2424   return kFALSE;
2425
2426 }
2427
2428 //_____________________________________________________________________________
2429 void AliTPC::LoadPoints2(Int_t)
2430 {
2431   //
2432   // Store x, y, z of all hits in memory
2433   //
2434   //  if (fTrackHits == 0 && fTrackHitsOld==0) return;
2435   if (fTrackHits == 0 ) return;
2436   //
2437   Int_t nhits =0;
2438   if (fHitType&4) nhits = fTrackHits->GetEntriesFast();
2439   //  if (fHitType&2) nhits = fTrackHitsOld->GetEntriesFast();
2440
2441   if (nhits == 0) return;
2442   Int_t tracks = gAlice->GetMCApp()->GetNtrack();
2443   if (fPoints == 0) fPoints = new TObjArray(tracks);
2444   AliHit *ahit;
2445   //
2446   Int_t *ntrk=new Int_t[tracks];
2447   Int_t *limi=new Int_t[tracks];
2448   Float_t **coor=new Float_t*[tracks];
2449   for(Int_t i=0;i<tracks;i++) {
2450     ntrk[i]=0;
2451     coor[i]=0;
2452     limi[i]=0;
2453   }
2454   //
2455   AliPoints *points = 0;
2456   Float_t *fp=0;
2457   Int_t trk;
2458   Int_t chunk=nhits/4+1;
2459   //
2460   // Loop over all the hits and store their position
2461   //
2462   ahit = FirstHit2(-1);
2463   while (ahit){
2464     trk=ahit->GetTrack();
2465     if(ntrk[trk]==limi[trk]) {
2466       //
2467       // Initialise a new track
2468       fp=new Float_t[3*(limi[trk]+chunk)];
2469       if(coor[trk]) {
2470         memcpy(fp,coor[trk],sizeof(Float_t)*3*limi[trk]);
2471         delete [] coor[trk];
2472       }
2473       limi[trk]+=chunk;
2474       coor[trk] = fp;
2475     } else {
2476       fp = coor[trk];
2477     }
2478     fp[3*ntrk[trk]  ] = ahit->X();
2479     fp[3*ntrk[trk]+1] = ahit->Y();
2480     fp[3*ntrk[trk]+2] = ahit->Z();
2481     ntrk[trk]++;
2482     ahit = NextHit2();
2483   }
2484
2485
2486
2487   //
2488   for(trk=0; trk<tracks; ++trk) {
2489     if(ntrk[trk]) {
2490       points = new AliPoints();
2491       points->SetMarkerColor(kYellow); //PH kYellow it the default in TPC
2492       points->SetMarkerSize(1);//PH Default size=1
2493       points->SetDetector(this);
2494       points->SetParticle(trk);
2495       points->SetPolyMarker(ntrk[trk],coor[trk],1); // Default style=1
2496       fPoints->AddAt(points,trk);
2497       delete [] coor[trk];
2498       coor[trk]=0;
2499     }
2500   }
2501   delete [] coor;
2502   delete [] ntrk;
2503   delete [] limi;
2504 }
2505
2506
2507 //_____________________________________________________________________________
2508 void AliTPC::LoadPoints3(Int_t)
2509 {
2510   //
2511   // Store x, y, z of all hits in memory
2512   // - only intersection point with pad row
2513   if (fTrackHits == 0) return;
2514   //
2515   Int_t nhits = fTrackHits->GetEntriesFast();
2516   if (nhits == 0) return;
2517   Int_t tracks = gAlice->GetMCApp()->GetNtrack();
2518   if (fPoints == 0) fPoints = new TObjArray(2*tracks);
2519   fPoints->Expand(2*tracks);
2520   AliHit *ahit;
2521   //
2522   Int_t *ntrk=new Int_t[tracks];
2523   Int_t *limi=new Int_t[tracks];
2524   Float_t **coor=new Float_t*[tracks];
2525   for(Int_t i=0;i<tracks;i++) {
2526     ntrk[i]=0;
2527     coor[i]=0;
2528     limi[i]=0;
2529   }
2530   //
2531   AliPoints *points = 0;
2532   Float_t *fp=0;
2533   Int_t trk;
2534   Int_t chunk=nhits/4+1;
2535   //
2536   // Loop over all the hits and store their position
2537   //
2538   ahit = FirstHit2(-1);
2539
2540   Int_t lastrow = -1;
2541   while (ahit){
2542     trk=ahit->GetTrack();
2543     Float_t  x[3]={ahit->X(),ahit->Y(),ahit->Z()};
2544     Int_t    index[3]={1,((AliTPChit*)ahit)->fSector,0};
2545     Int_t    currentrow = fTPCParam->GetPadRow(x,index) ;
2546     if (currentrow!=lastrow){
2547       lastrow = currentrow;
2548       //later calculate intersection point
2549       if(ntrk[trk]==limi[trk]) {
2550         //
2551         // Initialise a new track
2552         fp=new Float_t[3*(limi[trk]+chunk)];
2553         if(coor[trk]) {
2554           memcpy(fp,coor[trk],sizeof(Float_t)*3*limi[trk]);
2555           delete [] coor[trk];
2556         }
2557         limi[trk]+=chunk;
2558         coor[trk] = fp;
2559       } else {
2560         fp = coor[trk];
2561       }
2562       fp[3*ntrk[trk]  ] = ahit->X();
2563       fp[3*ntrk[trk]+1] = ahit->Y();
2564       fp[3*ntrk[trk]+2] = ahit->Z();
2565       ntrk[trk]++;
2566     }
2567     ahit = NextHit2();
2568   }
2569
2570   //
2571   for(trk=0; trk<tracks; ++trk) {
2572     if(ntrk[trk]) {
2573       points = new AliPoints();
2574       points->SetMarkerColor(kMagenta); //PH kYellow + 1 is kMagenta
2575       points->SetMarkerStyle(5);
2576       points->SetMarkerSize(0.2);
2577       points->SetDetector(this);
2578       points->SetParticle(trk);
2579       points->SetPolyMarker(ntrk[trk],coor[trk],30);
2580       fPoints->AddAt(points,tracks+trk);
2581       delete [] coor[trk];
2582       coor[trk]=0;
2583     }
2584   }
2585   delete [] coor;
2586   delete [] ntrk;
2587   delete [] limi;
2588 }
2589
2590
2591
2592 AliLoader* AliTPC::MakeLoader(const char* topfoldername)
2593 {
2594   //Makes TPC loader
2595   fLoader = new AliTPCLoader(GetName(),topfoldername);
2596   return fLoader;
2597 }
2598
2599 ////////////////////////////////////////////////////////////////////////
2600 AliTPCParam* AliTPC::LoadTPCParam(TFile *file) {
2601 //
2602 // load TPC paarmeters from a given file or create new if the object
2603 // is not found there
2604 // 12/05/2003 This method should be moved to the AliTPCLoader
2605 // and one has to decide where to store the TPC parameters
2606 // M.Kowalski
2607   char paramName[50];
2608   sprintf(paramName,"75x40_100x60_150x60");
2609   AliTPCParam *paramTPC=(AliTPCParam*)file->Get(paramName);
2610   if (paramTPC) {
2611     AliDebugClass(1,Form("TPC parameters %s found.",paramName));
2612   } else {
2613     AliWarningClass("TPC parameters not found. Create new (they may be incorrect)");
2614     paramTPC = new AliTPCParamSR;
2615   }
2616   return paramTPC;
2617
2618 // the older version of parameters can be accessed with this code.
2619 // In some cases, we have old parameters saved in the file but
2620 // digits were created with new parameters, it can be distinguish
2621 // by the name of TPC TreeD. The code here is just for the case
2622 // we would need to compare with old data, uncomment it if needed.
2623 //
2624 //  char paramName[50];
2625 //  sprintf(paramName,"75x40_100x60");
2626 //  AliTPCParam *paramTPC=(AliTPCParam*)in->Get(paramName);
2627 //  if (paramTPC) {
2628 //    cout<<"TPC parameters "<<paramName<<" found."<<endl;
2629 //  } else {
2630 //    sprintf(paramName,"75x40_100x60_150x60");
2631 //    paramTPC=(AliTPCParam*)in->Get(paramName);
2632 //    if (paramTPC) {
2633 //      cout<<"TPC parameters "<<paramName<<" found."<<endl;
2634 //    } else {
2635 //      cerr<<"TPC parameters not found. Create new (they may be incorrect)."
2636 //          <<endl;
2637 //      paramTPC = new AliTPCParamSR;
2638 //    }
2639 //  }
2640 //  return paramTPC;
2641
2642 }
2643
2644