TPC/AliTPCFast.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 //  fast TPC cluster simulation                                              //
  21 //                                                                          //
  22 ///////////////////////////////////////////////////////////////////////////////
  23
  24 #include <TParticle.h>
  25
  26 #include "AliRunLoader.h"
  27 #include "AliRun.h"
  28 #include "AliMC.h"
  29 #include "AliTPC.h"
  30 #include "AliTPCParam.h"
  31 #include "AliTPCClustersArray.h"
  32 #include "AliTPCClustersRow.h"
  33 #include "AliTPCcluster.h"
  34 #include "AliComplexCluster.h"
  35 #include "AliTPCFast.h"
  36 #include "AliLog.h"
  37
  38 ClassImp(AliTPCFast)
  39
  40
  41 //_____________________________________________________________________________
  42 void AliTPCFast::Hits2Clusters(AliRunLoader* runLoader) const
  43 {
  44   //--------------------------------------------------------
  45   // TPC simple cluster generator from hits
  46   // obtained from the TPC Fast Simulator
  47   // The point errors are taken from the parametrization
  48   //--------------------------------------------------------
  49
  50   //-----------------------------------------------------------------
  51   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
  52   //-----------------------------------------------------------------
  53   // Adopted to Marian's cluster data structure by I.Belikov, CERN,
  54   // Jouri.Belikov@cern.ch
  55   //----------------------------------------------------------------
  56
  57   /////////////////////////////////////////////////////////////////////////////
  58   //
  59   //---------------------------------------------------------------------
  60   //   ALICE TPC Cluster Parameters
  61   //--------------------------------------------------------------------
  62
  63
  64
  65   // Cluster width in rphi
  66   const Float_t kACrphi=0.18322;
  67   const Float_t kBCrphi=0.59551e-3;
  68   const Float_t kCCrphi=0.60952e-1;
  69   // Cluster width in z
  70   const Float_t kACz=0.19081;
  71   const Float_t kBCz=0.55938e-3;
  72   const Float_t kCCz=0.30428;
  73
  74
  75   AliLoader* loader = runLoader->GetLoader("TPCLoader");
  76   if (!loader) {
  77     AliError("No TPC loader found");
  78     return;
  79   }
  80   if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
  81   AliRun* aliRun = runLoader->GetAliRun();
  82   if (!aliRun) {
  83     AliError("Couldn't get AliRun object");
  84     return;
  85   }
  86   AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
  87   if (!tpc) {
  88     AliError("Couldn't get TPC detector");
  89     return;
  90   }
  91   AliTPCParam* param = tpc->GetParam();
  92   if (!param) {
  93     AliError("No TPC parameters available");
  94     return;
  95   }
  96
  97   //if(fDefaults == 0) SetDefaults();
  98
  99   Float_t sigmaRphi,sigmaZ,clRphi,clZ;
 100   //
 101   TParticle *particle; // pointer to a given particle
 102   AliTPChit *tpcHit; // pointer to a sigle TPC hit
 103   Int_t sector;
 104   Int_t ipart;
 105   Float_t xyz[5];
 106   Float_t pl,pt,tanth,rpad,ratio;
 107   Float_t cph,sph;
 108
 109   //---------------------------------------------------------------
 110   //  Get the access to the tracks
 111   //---------------------------------------------------------------
 112
 113   TTree *tH = loader->TreeH();
 114   if (tH == 0x0)
 115     AliFatal("Can not find TreeH in folder");
 116
 117   tpc->SetTreeAddress();
 118
 119   Stat_t ntracks = tH->GetEntries();
 120
 121   //Switch to the output file
 122
 123   if (loader->TreeR() == 0x0) loader->MakeTree("R");
 124
 125   AliDebug(1,Form("param->GetTitle() = %s",param->GetTitle()));
 126
 127   runLoader->CdGAFile();
 128   //param->Write(param->GetTitle());
 129
 130   AliTPCClustersArray carray;
 131   carray.Setup(param);
 132   carray.SetClusterType("AliTPCcluster");
 133   carray.MakeTree(loader->TreeR());
 134
 135   Int_t nclusters=0; //cluster counter
 136
 137   //------------------------------------------------------------
 138   // Loop over all sectors (72 sectors for 20 deg
 139   // segmentation for both lower and upper sectors)
 140   // Sectors 0-35 are lower sectors, 0-17 z>0, 17-35 z<0
 141   // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
 142   //
 143   // First cluster for sector 0 starts at "0"
 144   //------------------------------------------------------------
 145
 146   for(Int_t isec=0;isec<param->GetNSector();isec++){
 147     //MI change
 148     param->AdjustCosSin(isec,cph,sph);
 149
 150     //------------------------------------------------------------
 151     // Loop over tracks
 152     //------------------------------------------------------------
 153
 154     for(Int_t track=0;track<ntracks;track++){
 155       tpc->ResetHits();
 156       tH->GetEvent(track);
 157       //
 158       //  Get number of the TPC hits
 159       //
 160        tpcHit = (AliTPChit*)tpc->FirstHit(-1);
 161
 162       // Loop over hits
 163       //
 164        while(tpcHit){
 165
 166          if (tpcHit->fQ == 0.) {
 167            tpcHit = (AliTPChit*) tpc->NextHit();
 168            continue; //information about track (I.Belikov)
 169          }
 170         sector=tpcHit->fSector; // sector number
 171
 172        if(sector != isec){
 173          tpcHit = (AliTPChit*) tpc->NextHit();
 174          continue;
 175        }
 176         ipart=tpcHit->Track();
 177         particle=aliRun->GetMCApp()->Particle(ipart);
 178         pl=particle->Pz();
 179         pt=particle->Pt();
 180         if(pt < 1.e-9) pt=1.e-9;
 181         tanth=pl/pt;
 182         tanth = TMath::Abs(tanth);
 183         rpad=TMath::Sqrt(tpcHit->X()*tpcHit->X() + tpcHit->Y()*tpcHit->Y());
 184         ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
 185
 186         //   space-point resolutions
 187
 188         sigmaRphi=AliTPCcluster::SigmaY2(rpad,tanth,pt);
 189         sigmaZ   =AliTPCcluster::SigmaZ2(rpad,tanth   );
 190
 191         //   cluster widths
 192
 193         clRphi=kACrphi-kBCrphi*rpad*tanth+kCCrphi*ratio*ratio;
 194         clZ=kACz-kBCz*rpad*tanth+kCCz*tanth*tanth;
 195
 196         // temporary protection
 197
 198         if(sigmaRphi < 0.) sigmaRphi=0.4e-3;
 199         if(sigmaZ < 0.) sigmaZ=0.4e-3;
 200         if(clRphi < 0.) clRphi=2.5e-3;
 201         if(clZ < 0.) clZ=2.5e-5;
 202
 203         //
 204
 205         //
 206         // smearing --> rotate to the 1 (13) or to the 25 (49) sector,
 207         // then the inaccuracy in a X-Y plane is only along Y (pad row)!
 208         //
 209         Float_t xprim= tpcHit->X()*cph + tpcHit->Y()*sph;
 210         Float_t yprim=-tpcHit->X()*sph + tpcHit->Y()*cph;
 211         xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigmaRphi));   // y
 212           Float_t alpha=(isec < param->GetNInnerSector()) ?
 213           param->GetInnerAngle() : param->GetOuterAngle();
 214           Float_t ymax=xprim*TMath::Tan(0.5*alpha);
 215           if (TMath::Abs(xyz[0])>ymax) xyz[0]=yprim;
 216         xyz[1]=gRandom->Gaus(tpcHit->Z(),TMath::Sqrt(sigmaZ)); // z
 217           if (TMath::Abs(xyz[1])>param->GetZLength()) xyz[1]=tpcHit->Z();
 218         xyz[2]=sigmaRphi;                                     // fSigmaY2
 219         xyz[3]=sigmaZ;                                        // fSigmaZ2
 220         xyz[4]=tpcHit->fQ;                                    // q
 221
 222         AliTPCClustersRow *clrow=carray.GetRow(sector,tpcHit->fPadRow);
 223         if (!clrow) clrow=carray.CreateRow(sector,tpcHit->fPadRow);
 224
 225         Int_t tracks[3]={tpcHit->Track(), -1, -1};
 226         AliTPCcluster cluster(tracks,xyz);
 227
 228         clrow->InsertCluster(&cluster); nclusters++;
 229
 230         tpcHit = (AliTPChit*)tpc->NextHit();
 231
 232
 233       } // end of loop over hits
 234
 235     }   // end of loop over tracks
 236
 237     Int_t nrows=param->GetNRow(isec);
 238     for (Int_t irow=0; irow<nrows; irow++) {
 239         AliTPCClustersRow *clrow=carray.GetRow(isec,irow);
 240         if (!clrow) continue;
 241         carray.StoreRow(isec,irow);
 242         carray.ClearRow(isec,irow);
 243     }
 244
 245   } // end of loop over sectors
 246
 247   //  cerr<<"Number of made clusters : "<<nclusters<<"                        \n";
 248   loader->WriteRecPoints("OVERWRITE");
 249
 250
 251 } // end of function
 252
 253 //_________________________________________________________________
 254 void AliTPCFast::Hits2ExactClustersSector(AliRunLoader* runLoader,
 255                                           AliTPCClustersArray* clustersArray,
 256                                           Int_t isec) const
 257 {
 258   //--------------------------------------------------------
 259   //calculate exact cross point of track and given pad row
 260   //resulting values are expressed in "digit" coordinata
 261   //--------------------------------------------------------
 262
 263   //-----------------------------------------------------------------
 264   // Origin: Marian Ivanov  GSI Darmstadt, m.ivanov@gsi.de
 265   //-----------------------------------------------------------------
 266   //
 267   if (clustersArray==0){
 268     return;
 269   }
 270   AliLoader* loader = runLoader->GetLoader("TPCLoader");
 271   if (!loader) {
 272     AliError("No TPC loader found");
 273     return;
 274   }
 275   if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
 276   AliRun* aliRun = runLoader->GetAliRun();
 277   if (!aliRun) {
 278     AliError("Couldn't get AliRun object");
 279     return;
 280   }
 281   AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
 282   if (!tpc) {
 283     AliError("Couldn't get TPC detector");
 284     return;
 285   }
 286   AliTPCParam* param = tpc->GetParam();
 287   if (!param) {
 288     AliError("No TPC parameters available");
 289     return;
 290   }
 291   //
 292   TParticle *particle; // pointer to a given particle
 293   AliTPChit *tpcHit; // pointer to a sigle TPC hit
 294   //  Int_t sector,nhits;
 295   Int_t ipart;
 296   const Int_t kcmaxhits=30000;
 297   TVector * xxxx = new TVector(kcmaxhits*4);
 298   TVector & xxx = *xxxx;
 299   Int_t maxhits = kcmaxhits;
 300   //construct array for each padrow
 301   for (Int_t i=0; i<param->GetNRow(isec);i++)
 302     clustersArray->CreateRow(isec,i);
 303
 304   //---------------------------------------------------------------
 305   //  Get the access to the tracks
 306   //---------------------------------------------------------------
 307
 308   TTree *tH = loader->TreeH();
 309   if (tH == 0x0)
 310     AliFatal("Can not find TreeH in folder");
 311
 312   tpc->SetTreeAddress();
 313
 314   Stat_t ntracks = tH->GetEntries();
 315   Int_t npart = aliRun->GetMCApp()->GetNtrack();
 316   //MI change
 317   TBranch * branch=0;
 318   if (tpc->GetHitType()>1) branch = tH->GetBranch("TPC2");
 319   else branch = tH->GetBranch("TPC");
 320
 321   //------------------------------------------------------------
 322   // Loop over tracks
 323   //------------------------------------------------------------
 324
 325   for(Int_t track=0;track<ntracks;track++){
 326     Bool_t isInSector=kTRUE;
 327     tpc->ResetHits();
 328      isInSector = tpc->TrackInVolume(isec,track);
 329     if (!isInSector) continue;
 330     //MI change
 331     branch->GetEntry(track); // get next track
 332     //
 333     //  Get number of the TPC hits and a pointer
 334     //  to the particles
 335     // Loop over hits
 336     //
 337     Int_t currentIndex=0;
 338     Int_t lastrow=-1;  //last writen row
 339
 340     //M.I. changes
 341
 342     tpcHit = (AliTPChit*)tpc->FirstHit(-1);
 343     while(tpcHit){
 344
 345       Int_t sector=tpcHit->fSector; // sector number
 346       if(sector != isec){
 347         tpcHit = (AliTPChit*) tpc->NextHit();
 348         continue;
 349       }
 350
 351       ipart=tpcHit->Track();
 352       if (ipart<npart) particle=aliRun->GetMCApp()->Particle(ipart);
 353
 354       //find row number
 355
 356       Float_t  x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
 357       Int_t    index[3]={1,isec,0};
 358       Int_t    currentrow = param->GetPadRow(x,index) ;
 359       if (currentrow<0) {tpcHit = (AliTPChit*)tpc->NextHit(); continue;}
 360       if (lastrow<0) lastrow=currentrow;
 361       if (currentrow==lastrow){
 362         if ( currentIndex>=maxhits){
 363           maxhits+=kcmaxhits;
 364           xxx.ResizeTo(4*maxhits);
 365         }
 366         xxx(currentIndex*4)=x[0];
 367         xxx(currentIndex*4+1)=x[1];
 368         xxx(currentIndex*4+2)=x[2];
 369         xxx(currentIndex*4+3)=tpcHit->fQ;
 370         currentIndex++;
 371       }
 372       else
 373         if (currentIndex>2){
 374           Float_t sumx=0;
 375           Float_t sumx2=0;
 376           Float_t sumx3=0;
 377           Float_t sumx4=0;
 378           Float_t sumy=0;
 379           Float_t sumxy=0;
 380           Float_t sumx2y=0;
 381           Float_t sumz=0;
 382           Float_t sumxz=0;
 383           Float_t sumx2z=0;
 384           Float_t sumq=0;
 385           for (Int_t index=0;index<currentIndex;index++){
 386             Float_t x,x2,x3,x4;
 387             x=x2=x3=x4=xxx(index*4);
 388             x2*=x;
 389             x3*=x2;
 390             x4*=x3;
 391             sumx+=x;
 392             sumx2+=x2;
 393             sumx3+=x3;
 394             sumx4+=x4;
 395             sumy+=xxx(index*4+1);
 396             sumxy+=xxx(index*4+1)*x;
 397             sumx2y+=xxx(index*4+1)*x2;
 398             sumz+=xxx(index*4+2);
 399             sumxz+=xxx(index*4+2)*x;
 400             sumx2z+=xxx(index*4+2)*x2;
 401             sumq+=xxx(index*4+3);
 402           }
 403           Float_t centralPad = (param->GetNPads(isec,lastrow)-1)/2;
 404           Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
 405             sumx2*(sumx*sumx3-sumx2*sumx2);
 406
 407           Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
 408             sumx2*(sumxy*sumx3-sumx2y*sumx2);
 409           Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
 410             sumx2*(sumxz*sumx3-sumx2z*sumx2);
 411
 412           Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
 413             sumx2*(sumx*sumx2y-sumx2*sumxy);
 414           Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
 415             sumx2*(sumx*sumx2z-sumx2*sumxz);
 416
 417           if (TMath::Abs(det)<0.00001){
 418              tpcHit = (AliTPChit*)tpc->NextHit();
 419             continue;
 420           }
 421
 422           Float_t y=detay/det+centralPad;
 423           Float_t z=detaz/det;
 424           Float_t by=detby/det; //y angle
 425           Float_t bz=detbz/det; //z angle
 426           sumy/=Float_t(currentIndex);
 427           sumz/=Float_t(currentIndex);
 428
 429           AliTPCClustersRow * row = (clustersArray->GetRow(isec,lastrow));
 430           if (row!=0) {
 431             AliComplexCluster* cl = new((AliComplexCluster*)row->Append()) AliComplexCluster ;
 432             //    AliComplexCluster cl;
 433             cl->fX=z;
 434             cl->fY=y;
 435             cl->fQ=sumq;
 436             cl->fSigmaX2=bz;
 437             cl->fSigmaY2=by;
 438             cl->fTracks[0]=ipart;
 439           }
 440           currentIndex=0;
 441           lastrow=currentrow;
 442         } //end of calculating cluster for given row
 443
 444
 445       tpcHit = (AliTPChit*)tpc->NextHit();
 446     } // end of loop over hits
 447   }   // end of loop over tracks
 448   //write padrows to tree
 449   for (Int_t ii=0; ii<param->GetNRow(isec);ii++) {
 450     clustersArray->StoreRow(isec,ii);
 451     clustersArray->ClearRow(isec,ii);
 452   }
 453   xxxx->Delete();
 454
 455 }
 456
 457
 458 void AliTPCFast::FindTrackHitsIntersection(AliRunLoader* runLoader,
 459                                            AliTPCClustersArray* clustersArray) const
 460 {
 461
 462   //
 463   //fill clones array with intersection of current point with the
 464   //middle of the row
 465   AliLoader* loader = runLoader->GetLoader("TPCLoader");
 466   if (!loader) {
 467     AliError("No TPC loader found");
 468     return;
 469   }
 470   if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
 471   AliRun* aliRun = runLoader->GetAliRun();
 472   if (!aliRun) {
 473     AliError("Couldn't get AliRun object");
 474     return;
 475   }
 476   AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
 477   if (!tpc) {
 478     AliError("Couldn't get TPC detector");
 479     return;
 480   }
 481   AliTPCParam* param = tpc->GetParam();
 482   if (!param) {
 483     AliError("No TPC parameters available");
 484     return;
 485   }
 486
 487   Int_t sector;
 488   Int_t ipart;
 489
 490   const Int_t kcmaxhits=30000;
 491   TVector * xxxx = new TVector(kcmaxhits*4);
 492   TVector & xxx = *xxxx;
 493   Int_t maxhits = kcmaxhits;
 494
 495   //
 496   AliTPChit * tpcHit=0;
 497   tpcHit = (AliTPChit*)tpc->FirstHit2(-1);
 498   Int_t currentIndex=0;
 499   Int_t lastrow=-1;  //last writen row
 500
 501   while (tpcHit){
 502     if (tpcHit==0) continue;
 503     sector=tpcHit->fSector; // sector number
 504     ipart=tpcHit->Track();
 505
 506     //find row number
 507
 508     Float_t  x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
 509     Int_t    index[3]={1,sector,0};
 510     Int_t    currentrow = param->GetPadRow(x,index) ;
 511     if (currentrow<0) continue;
 512     if (lastrow<0) lastrow=currentrow;
 513     if (currentrow==lastrow){
 514       if ( currentIndex>=maxhits){
 515         maxhits+=kcmaxhits;
 516         xxx.ResizeTo(4*maxhits);
 517       }
 518       xxx(currentIndex*4)=x[0];
 519       xxx(currentIndex*4+1)=x[1];
 520       xxx(currentIndex*4+2)=x[2];
 521       xxx(currentIndex*4+3)=tpcHit->fQ;
 522       currentIndex++;
 523     }
 524     else
 525       if (currentIndex>2){
 526         Float_t sumx=0;
 527         Float_t sumx2=0;
 528         Float_t sumx3=0;
 529         Float_t sumx4=0;
 530         Float_t sumy=0;
 531         Float_t sumxy=0;
 532         Float_t sumx2y=0;
 533         Float_t sumz=0;
 534         Float_t sumxz=0;
 535         Float_t sumx2z=0;
 536         Float_t sumq=0;
 537         for (Int_t index=0;index<currentIndex;index++){
 538           Float_t x,x2,x3,x4;
 539           x=x2=x3=x4=xxx(index*4);
 540           x2*=x;
 541           x3*=x2;
 542           x4*=x3;
 543           sumx+=x;
 544           sumx2+=x2;
 545           sumx3+=x3;
 546           sumx4+=x4;
 547           sumy+=xxx(index*4+1);
 548           sumxy+=xxx(index*4+1)*x;
 549           sumx2y+=xxx(index*4+1)*x2;
 550           sumz+=xxx(index*4+2);
 551           sumxz+=xxx(index*4+2)*x;
 552           sumx2z+=xxx(index*4+2)*x2;
 553           sumq+=xxx(index*4+3);
 554         }
 555         Float_t centralPad = (param->GetNPads(sector,lastrow)-1)/2;
 556         Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
 557           sumx2*(sumx*sumx3-sumx2*sumx2);
 558
 559         Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
 560           sumx2*(sumxy*sumx3-sumx2y*sumx2);
 561         Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
 562           sumx2*(sumxz*sumx3-sumx2z*sumx2);
 563
 564         Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
 565           sumx2*(sumx*sumx2y-sumx2*sumxy);
 566         Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
 567           sumx2*(sumx*sumx2z-sumx2*sumxz);
 568
 569         Float_t y=detay/det+centralPad;
 570         Float_t z=detaz/det;
 571         Float_t by=detby/det; //y angle
 572         Float_t bz=detbz/det; //z angle
 573         sumy/=Float_t(currentIndex);
 574         sumz/=Float_t(currentIndex);
 575
 576         AliComplexCluster cl;
 577         cl.fX=z;
 578         cl.fY=y;
 579         cl.fQ=sumq;
 580         cl.fSigmaX2=bz;
 581         cl.fSigmaY2=by;
 582         cl.fTracks[0]=ipart;
 583
 584         AliTPCClustersRow * row = (clustersArray->GetRow(sector,lastrow));
 585         if (row!=0) row->InsertCluster(&cl);
 586         currentIndex=0;
 587         lastrow=currentrow;
 588       } //end of calculating cluster for given row
 589
 590   } // end of loop over hits
 591   xxxx->Delete();
 592
 593 }
 594