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