Transition to NewIO
[u/mrichter/AliRoot.git] / ITS / AliITSFindClustersV2.C
1 /****************************************************************
2 *  This macro converts AliITSRecPoint(s) to AliITSclusterV2(s)  *
3 *           Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch      *
4 *****************************************************************/
5
6 #ifndef __CINT__
7   #include <Riostream.h>
8
9   #include "AliRun.h"
10   #include "AliRunLoader.h"
11   #include "AliLoader.h"
12   #include "AliITS.h"
13   #include "AliITSgeom.h"
14   #include "AliITSclustererV2.h"
15
16   #include "TStopwatch.h"
17 #endif
18
19 Int_t AliITSFindClustersV2(Char_t SlowOrFast='f')
20 {
21
22     cerr<<"AliITSRecPoint(s) -> AliITSclusterV2(s)...\n";
23     
24     if (gAlice) 
25      {
26       delete gAlice->GetRunLoader();
27       delete gAlice; 
28       gAlice=0;
29      }
30  
31     AliRunLoader* rl = AliRunLoader::Open("galice.root");
32     if (rl == 0x0)
33      {
34       cerr<<"AliITSHits2DigitsDefault.C : Can not open session RL=NULL"
35            << endl;
36        return 3;
37      }
38      
39     Int_t retval = rl->LoadgAlice();
40     if (retval)
41      {
42       cerr<<"AliITSHits2DigitsDefault.C : LoadgAlice returned error"
43            << endl;
44        delete rl;
45        return 3;
46      }
47     gAlice=rl->GetAliRun();
48     rl->LoadHeader();
49     retval = rl->LoadKinematics();
50     if (retval)
51      {
52       cerr<<"AliITSHits2DigitsDefault.C : LoadKinematics returned error"
53            << endl;
54        delete rl;
55        return 3;
56      }
57     
58     AliITSLoader* gime = (AliITSLoader*)rl->GetLoader("ITSLoader");
59     if (gime == 0x0)
60      {
61       cerr<<"AliITSHits2DigitsDefault.C : can not get ITS loader"
62            << endl;
63      }
64
65    rl->GetEvent(0);
66
67    AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
68    if (!ITS) { cerr<<"Can't find the ITS !\n"; delete rl; return 3; }
69    AliITSgeom *geom=ITS->GetITSgeom();
70  
71    TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
72   
73    gime->LoadRawClusters("recreate");
74
75    if (SlowOrFast=='f') 
76     {
77        gime->SetRecPointsFileName("ITS.FastRecPoints.root");
78     }
79    if (gime->LoadRecPoints())
80     {
81       cerr<<"Load Rec Pints returned error !\n"; 
82       delete rl;
83       return 4;
84     }
85
86    TClonesArray *points = new TClonesArray("AliITSRecPoint",10000);
87
88    Float_t lp[5]; 
89    Int_t lab[6];
90    
91    Int_t iEvent;
92    for (iEvent = 0; iEvent< rl->GetNumberOfEvents() ; iEvent++) 
93     {
94
95      rl->GetEvent(iEvent);
96
97      TTree *cTree = gime->TreeC();
98      if (cTree == 0x0)  
99       {
100        gime->MakeTree("C");
101        cTree = gime->TreeC();
102       }
103    
104      cTree->Branch("Clusters",&clusters);
105      
106      TTree *pTree=gime->TreeR();
107      if (pTree == 0x0) 
108       {
109         cerr<<"Can not get TreeR !\n"; delete rl;return 5;
110       }
111      TBranch *branch = 0;
112      if (SlowOrFast=='f') {
113        branch = pTree->GetBranch("ITSRecPointsF");
114      }
115      else {
116        branch = pTree->GetBranch("ITSRecPoints");
117      }
118      if (!branch) 
119       { 
120        cerr<<"Can't get ITSRecPoints branch !\n"; 
121        delete rl;
122        return 6; 
123       }
124     
125      branch->SetAddress(&points);
126   
127      AliStack* stack = rl->Stack();
128      if (stack == 0x0)
129       {
130        cerr<<"AliITSFindClustersV2.C : Can not get stack"
131            << endl;
132        delete rl;
133        return 3;
134       }
135
136      TClonesArray &cl=*clusters;
137      Int_t nclusters=0;
138      Int_t nentr=(Int_t)branch->GetEntries();
139
140      cerr<<"Number of entries: "<<nentr<<endl;
141
142      for (Int_t i=0; i<nentr; i++) 
143       {
144        points->Clear();
145        clusters->Clear();
146        branch->GetEvent(i);
147        Int_t ncl=points->GetEntriesFast(); if (ncl==0){cTree->Fill();continue;}
148        Int_t lay,lad,det; geom->GetModuleId(i,lay,lad,det);
149        if ( (lay<0) || (lad<0) || (det<0))
150         {
151           ::Error("AliITSFindClustersV2.C","No such a module %d",i);
152           continue;
153         }
154        Float_t x,y,zshift; geom->GetTrans(lay,lad,det,x,y,zshift); 
155        Double_t rot[9];    geom->GetRotMatrix(lay,lad,det,rot);
156        Double_t yshift = x*rot[0] + y*rot[1];
157        Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1);
158        nclusters+=ncl;
159
160        Float_t kmip=1; // ADC->mip normalization factor for the SDD and SSD 
161        if(lay==4 || lay==3){kmip=280.;};
162        if(lay==6 || lay==5){kmip=38.;};
163
164        for (Int_t j=0; j<ncl; j++) 
165         {
166           AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(j);
167           //Float_t lp[5];
168           lp[0]=-p->GetX()-yshift; if (lay==1) lp[0]=-lp[0];
169           lp[1]=p->GetZ()+zshift;
170           lp[2]=p->GetSigmaX2();
171           lp[3]=p->GetSigmaZ2();
172           lp[4]=p->GetQ(); lp[4]/=kmip;
173           //Int_t lab[6]; 
174           lab[0]=p->GetLabel(0);lab[1]=p->GetLabel(1);lab[2]=p->GetLabel(2);
175           lab[3]=ndet;
176
177           Int_t label=lab[0];
178           if (label>=0) {
179              TParticle *part=(TParticle*)stack->Particle(label);
180              if (part == 0x0)
181               cerr<<"Can not get particle with label "<<label<<endl;
182              label=-3;
183              while (part->P() < 0.005) {
184                 Int_t m=part->GetFirstMother();
185                 if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;}
186                 label=m;
187                 part=(TParticle*)gAlice->Particle(label);
188              }
189              if      (lab[1]<0) lab[1]=label;
190              else if (lab[2]<0) lab[2]=label;
191              else cerr<<"No empty labels !\n";
192           }
193
194           new(cl[j]) AliITSclusterV2(lab,lp);
195         }
196        cTree->Fill(); 
197 //       clusters->Delete(); points->Delete();
198      }
199     gime->WriteRawClusters("OVERWRITE");
200     cerr<<"Number of clusters: "<<nclusters<<endl;
201     
202    }
203    
204    delete clusters;
205    delete points;
206    delete rl;
207    return 0;
208
209 }
210
211
212
213
214
215
216
217
218
219
220
221
222
223