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