TPC/AliTPC.cxx

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