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