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