Typo corrected (HP)
[u/mrichter/AliRoot.git] / ITS / AliITSFindClustersV2.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 /*
17 $Log$
18 Revision 1.1  2002/09/09 17:36:05  nilsen
19 new TTask to replace non-working AliITSFindClusterV2.C macro.
20
21 */
22 #include <TROOT.h>
23 #include <TFile.h>
24 #include <TTree.h>
25 #include <TBranch.h>
26 #include <TClonesArray.h>
27 #include <TParticle.h>
28
29 #include "AliRun.h"
30 #include "AliHeader.h"
31
32 #include "AliITS.h"
33 #include "AliITSRecPoint.h"
34 #include "AliITSFindClustersV2.h"
35 #include "AliITSclusterV2.h"
36 #include "AliITSgeom.h"
37
38 ClassImp(AliITSFindClustersV2)
39
40 //______________________________________________________________________
41 AliITSFindClustersV2::AliITSFindClustersV2(){
42     // Default constructor.
43     // Inputs:
44     //  none.
45     // Outputs:
46     //   none.
47     // Return:
48     //    A zero-ed constructed AliITSFindClustersV2 class.
49
50     fAr          = 0;
51     fDeletfAr     = kFALSE; // fAr=0 dont't delete it.
52     fGeom        = 0;
53     fInFileName  = 0;
54     fOutFileName = 0;
55     fIn          = 0;
56     fOut         = 0;
57     fSlowFast    = kFALSE;  // slow simulation
58     fInit        = kFALSE;  // Init failed
59 }
60 //______________________________________________________________________
61 AliITSFindClustersV2::AliITSFindClustersV2(const TString infile,
62                                            const TString outfile){
63     // Standard constructor.
64     // Inputs:
65     //  const TString infile   Input file name where the RecPoints are
66     //                         to be read from.
67     //  const TString outfile  Output file where V2 tracking clulsters
68     //                         are to be written. if =="" writen to the
69     //                         same file as input.
70     // Outputs:
71     //   none.
72     // Return:
73     //    A standardly constructed AliITSFindClustersV2 class.
74
75     fAr          = 0;
76     fDeletfAr    = kFALSE; // fAr=0 dont't delete it.
77     fGeom        = 0;
78     fInFileName  = 0;
79     fOutFileName = 0;
80     fIn          = 0;
81     fOut         = 0;
82     fSlowFast    = kFALSE;  // slow simulation
83     fInit        = kFALSE;  // Init failed
84
85     fInFileName = new TString(infile);
86     if(outfile.CompareTo("")!=0){
87         fOutFileName = new TString(outfile);
88     } // end if outfile.CompareeTo("")!=0
89     
90     if(fOutFileName!=0){
91         fIn = new TFile(fInFileName->Data(),"READ");
92         fOut = new TFile(fOutFileName->Data(),"UPDATE");
93     }else{ // open fIn file for update
94         fIn = new TFile(fInFileName->Data(),"UPDATE");
95     } // end if
96
97     fAr  = (AliRun*) fIn->Get("gAlice");
98     if(!fAr){
99         Warning("AliITSFindClusterV2",
100                 "Can't fine gAlice in file %s. Aborting.",fIn->GetName());
101         return;
102     } // end if !fAr
103     fDeletfAr = kTRUE; // Since gAlice was read in, delete it.
104
105     AliITS *its = (AliITS*) fAr->GetModule("ITS");
106     if(!its){
107         Warning("AliITSFindClusterV2",
108                 "Can't fine the ITS in gAlice. Aborting.");
109         return;
110     } // end if !its
111     fGeom = its->GetITSgeom();
112     if(!fGeom){
113         Warning("AliITSFindClusterV2",
114                 "Can't fine the ITS geometry in gAlice. Aborting.");
115         return;
116     } // end if !fGeom
117
118     if(fOut) fOut->cd();
119     fInit = kTRUE;
120 }
121 //______________________________________________________________________
122 AliITSFindClustersV2::AliITSFindClustersV2(TFile *in,
123                                            TFile *out){
124     // Standard constructor.
125     // Inputs:
126     //  const TFile *in   Input file where the RecPoints are
127     //                         to be read from.
128     //  const Tfile *out  Output file where V2 tracking clulsters
129     //                         are to be written. if ==0 writen to the
130     //                         same file as input.
131     // Outputs:
132     //   none.
133     // Return:
134     //    A standardly constructed AliITSFindClustersV2 class.
135
136     fAr          = 0;
137     fDeletfAr    = kFALSE; // fAr=0 dont't delete it.
138     fGeom        = 0;
139     fInFileName  = 0;
140     fOutFileName = 0;
141     fIn          = 0;
142     fOut         = 0;
143     fSlowFast    = kFALSE;  // slow simulation
144     fInit        = kFALSE;  // Init failed
145
146     fIn  = in;
147     fOut = out;
148     fAr  = (AliRun*) fIn->Get("gAlice");
149     if(!fAr){
150         Warning("AliITSFindClusterV2",
151                 "Can't fine gAlice in file %s. Aborting.",fIn->GetName());
152         return;
153     } // end if !fAr
154     fDeletfAr = kTRUE; // Since gAlice was read in, delete it.
155     AliITS *its = (AliITS*) fAr->GetModule("ITS");
156     if(!its){
157         Warning("AliITSFindClusterV2",
158                 "Can't fine the ITS in gAlice. Aborting.");
159         return;
160     } // end if !its
161     fGeom = its->GetITSgeom();
162     if(!fGeom){
163         Warning("AliITSFindClusterV2",
164                 "Can't fine the ITS geometry in gAlice. Aborting.");
165         return;
166     } // end if !fGeom
167
168     if(fOut) fOut->cd();
169     fInit = kTRUE;
170 }
171 //______________________________________________________________________
172 AliITSFindClustersV2::AliITSFindClustersV2(AliRun *ar,
173                                            const TString outfile){
174     // Standard constructor.
175     // Inputs:
176     //  const TString infile   Input file name where the RecPoints are
177     //                         to be read from.
178     //  const TString outfile  Output file where V2 tracking clulsters
179     //                         are to be written. if =="" writen to the
180     //                         same file as input.
181     // Outputs:
182     //   none.
183     // Return:
184     //    A standardly constructed AliITSFindClustersV2 class.
185
186     fAr          = 0;
187     fDeletfAr    = kFALSE; // fAr=0 dont't delete it.
188     fGeom        = 0;
189     fInFileName  = 0;
190     fOutFileName = 0;
191     fIn          = 0;
192     fOut         = 0;
193     fSlowFast    = kFALSE;  // slow simulation
194     fInit        = kFALSE;  // Init failed
195
196     if(outfile.CompareTo("")!=0){
197         fOutFileName = new TString(outfile);
198     } // end if outfile.CompareeTo("")!=0
199     
200     if(fOutFileName!=0){
201         fOut = new TFile(fOutFileName->Data(),"UPDATE");
202     } // end if
203
204     fAr  = ar;
205     if(!fAr){
206         Warning("AliITSFindClusterV2",
207                 "ar==0. Aborting.");
208         return;
209     } // end if !fAr
210     AliITS *its = (AliITS*) fAr->GetModule("ITS");
211     if(!its){
212         Warning("AliITSFindClusterV2",
213                 "Can't fine the ITS in gAlice. Aborting.");
214         return;
215     } // end if !its
216     fGeom = its->GetITSgeom();
217     if(!fGeom){
218         Warning("AliITSFindClusterV2",
219                 "Can't fine the ITS geometry in gAlice. Aborting.");
220         return;
221     } // end if !fGeom
222
223     if(fOut) fOut->cd();
224     fInit = kTRUE;
225 }
226 //______________________________________________________________________
227 AliITSFindClustersV2::~AliITSFindClustersV2(){
228     // Default constructor.
229     // Inputs:
230     //  none.
231     // Outputs:
232     //   none.
233     // Return:
234     //    A destroyed AliITSFindclustersV2 class.
235
236     fGeom = 0; // Deleted when AliRun/ITS is deleted.
237     if(fInFileName!=0){ // input file opened by AliITSFindClustersV2
238         if(fIn!=0){
239             if(fIn->IsOpen()) fIn->Close();
240             delete fIn;
241             fIn = 0;
242         } // end if
243         delete fInFileName;
244         fInFileName = 0;
245     } // end if
246
247     if(fOutFileName!=0){ // input file opened by AliITSFindClustersV2
248         if(fOut!=0){
249             if(fOut->IsOpen()) fOut->Close();
250             delete fOut;
251             fOut = 0;
252         } // end if
253         delete fOutFileName;
254         fOutFileName = 0;
255     } // end if
256     if(fDeletfAr && !fAr){
257         cout << "deleting fAr."<< endl;
258         delete fAr;
259         fAr = 0;
260         cout << "fAr deleted OK."<< endl;
261     } // end if fDeletfAr
262 }
263 //______________________________________________________________________ 
264 void AliITSFindClustersV2::Exec(const Option_t *opt){
265     // Main FindclustersV2 function.
266     // Inputs:
267     //      Option_t * opt   list of subdetector to digitize. =0 all.
268     // Outputs:
269     //      none.
270     // Return:
271     //      none.
272     Char_t name[50];
273
274     if(!fInit){
275         Warning("Exec","Initilization not succesfull. Aborting.");
276         return;
277     } // end if !fInit
278
279     fGeom->Write();
280
281     fAr->GetEvent(0);
282     TTree *pTree = fAr->TreeR();
283     if(!pTree){
284         Warning("Exec","Error getting TreeR. TreeR=%p",pTree);
285         return;
286     } // end if !pTree
287     TBranch *branch = 0;
288     if(fSlowFast) sprintf(name,"ITSRecPointsF");
289     else sprintf(name,"ITSRecPoints");
290     branch = pTree->GetBranch(name);
291     if(!branch){
292         Warning("Exec","Can't find branch %s in TreeR fSlowFast=%d",
293                 name,fSlowFast);
294         return;
295     } // end if !branch
296     TClonesArray *points = new TClonesArray("AliITSRecPoint",10000);
297     branch->SetAddress(&points);
298     Int_t nentr = (Int_t) branch->GetEntries();
299
300     if(fOut) fOut->cd();
301     TClonesArray *cluster = new TClonesArray("AliITSclusterV2",10000);
302     sprintf(name,"TreeC_ITS_%d",fAr->GetHeader()->GetEvent());
303     TTree *cTree = new TTree(name,"ITS clusters");
304     cTree->Branch("Clusters",&cluster);
305     TClonesArray &cl = *cluster;
306
307     Float_t lp[5];
308     Int_t lab[6];
309     Int_t i,j,lay,lad,det,nclusters=0,ncl;
310     Float_t kmip,x,y,zshift,yshift;
311     Double_t rot[9];
312     AliITSRecPoint *p;
313     TParticle *part;
314
315     for(i=0;i<nentr;i++){
316         points->Clear();
317         branch->GetEvent(i);
318         ncl = points->GetEntriesFast();
319         if(ncl<=0) {
320             cTree->Fill();
321             continue;
322         } // end if ncl<=0
323         nclusters += ncl;
324         fGeom->GetModuleId(i,lay,lad,det);
325         fGeom->GetTrans(i,x,y,zshift);
326         fGeom->GetRotMatrix(i,rot);
327         yshift = x*rot[0] + y*rot[1];
328         kmip = 1.0; // fGeom->GetModuletype(i)==kSPD
329         if(fGeom->GetModuleType(i)==kSDD) kmip = 280.0;
330         if(fGeom->GetModuleType(i)==kSSD) kmip = 38.0;
331         for(j=0;j<ncl;j++){
332             p = (AliITSRecPoint*)(points->UncheckedAt(j));
333             lp[0] = - (p->GetX() + yshift);
334             if(lay==1) lp[0] = -lp[0];
335             lp[1] = p->GetZ() + zshift;
336             lp[2] = p->GetSigmaX2();
337             lp[3] = p->GetSigmaZ2();
338             lp[4] = p->GetQ()/kmip;
339             lab[0] = p->GetLabel(0);
340             lab[1] = p->GetLabel(1);
341             lab[2] = p->GetLabel(2);
342             lab[3] = i;
343             lad = lab[0];
344             if(lad>=0){
345                 part = (TParticle*) fAr->Particle(lad);
346                 lad = -3;
347                 while(part->P() < 0.005){
348                     if(part->GetFirstMother()<0){
349                         Warning("Exec","Primary momentum: %e",part->P());
350                         break;
351                     } // end if part->GetFirstMother()<0
352                     lad = part->GetFirstMother();
353                     part = (TParticle*) fAr->Particle(lad);
354                 } // end while part->P() < 0.005
355                 if(lab[1]<0) lab[1] = lad;
356                 else if(lab[2]<0) lab[2] = lad;
357                 else Warning("Exec","No empty lables!");
358             } // end if lab>=0
359             new(cl[j]) AliITSclusterV2(lab,lp);
360         } // end for j
361         cTree->Fill();
362         cluster->Delete();
363         points->Delete();
364     } // end for i
365     cTree->Write();
366     
367     delete cTree;
368     delete cluster;
369     delete points;
370 }