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