d915e30d956b937c75fcd82611308f44d2b97c93
[u/mrichter/AliRoot.git] / ITS / AliITSFindClustersV2.C
1 #ifndef __CINT__
2   #include <iostream.h>
3
4   #include "AliRun.h"
5   #include "AliITS.h"
6   #include "AliITSgeom.h"
7   #include "AliITSRecPoint.h"
8   #include "AliITSclusterV2.h"
9   #include "AliITSsimulationFastPoints.h"
10
11   #include "TFile.h"
12   #include "TTree.h"
13   #include "TParticle.h"
14 #endif
15
16 //Int_t AliITSFindClustersV2() {
17 Int_t AliITSFindClustersV2(char SlowOrFast='s') {
18
19 /****************************************************************
20  *  Just something to start with                                *
21  ****************************************************************/
22    cerr<<"Looking for clusters...\n";
23    cout<<"!!!! SlowOrFast =  "<<SlowOrFast<<endl;
24 ///////////////// Conversion AliITSRecPoint -> AliITSclusterV2 //////////////
25
26    //    for the fast recpoints
27   if(SlowOrFast=='f') {
28    cout<<"22 !!!! SlowOrFast =  "<<SlowOrFast<<endl; 
29    /* 
30    if (gAlice) {delete gAlice; gAlice=0;}
31
32    TFile *in=TFile::Open("galice.root","update");
33    if (!in->IsOpen()) {
34       cerr<<"Can't open galice.root !\n"; 
35       return 1;
36    }
37    
38    if (!(gAlice=(AliRun*)in->Get("gAlice"))) {
39       cerr<<"Can't find gAlice !\n";
40       return 2;
41    }
42
43    Int_t ev=0;
44    gAlice->GetEvent(ev);
45
46    AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
47    if (!ITS) {
48       cerr<<"Can't find the ITS !\n";
49       return 3;
50    }
51
52    gAlice->MakeTree("R"); ITS->MakeBranch("R",0);
53 //////////////// Taken from ITSHitsToFastPoints.C ///////////////////////
54    ITS->SetSimulationModel(0,new AliITSsimulationFastPoints());
55    ITS->SetSimulationModel(1,new AliITSsimulationFastPoints());
56    ITS->SetSimulationModel(2,new AliITSsimulationFastPoints());
57    Int_t nsignal=25;
58    Int_t size=-1;
59    Int_t bgr_ev=Int_t(ev/nsignal);
60    ITS->HitsToFastRecPoints(ev,bgr_ev,size," ","All"," ");
61 //////////////////////////////////////////////////////////////////////////
62
63    delete gAlice; gAlice=0;
64    in->Close();
65    */
66   }    // end of fast recpoints
67
68
69    // for the slow points
70
71    /*TFile */in=TFile::Open("galice.root");
72
73    if (gAlice) {delete gAlice; gAlice=0;}
74
75    if (!(gAlice=(AliRun*)in->Get("gAlice"))) {
76       cerr<<"Can't find gAlice !\n";
77       return 4;
78    }
79
80    gAlice->GetEvent(0);
81
82    /*AliITS */ITS  = (AliITS*)gAlice->GetModule("ITS");
83    if (!ITS) {
84       cerr<<"Can't find the ITS !\n";
85       return 5;
86    }
87    AliITSgeom *geom=ITS->GetITSgeom();
88  
89    TFile *out=TFile::Open("AliITSclustersV2.root","new");
90    if (!out->IsOpen()) {
91       cerr<<"Delete old AliITSclustersV2.root !\n"; 
92       return 6;
93    }
94    geom->Write();
95
96    TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
97    //TTree *cTree=new TTree("cTree","ITS clusters");
98    TTree *cTree=new TTree("TreeC_ITS_0","ITS clusters");
99    cTree->Branch("Clusters",&clusters);
100
101    TTree *pTree=gAlice->TreeR();
102    if (!pTree) {
103       cerr<<"Can't get TreeR !\n";
104       return 7;
105    }
106    TBranch *branch=pTree->GetBranch("ITSRecPoints");
107    if (!branch) { 
108       cerr<<"Can't get ITSRecPoints branch !\n";
109       return 8;
110    }
111    TClonesArray *points=new TClonesArray("AliITSRecPoint",10000);
112    branch->SetAddress(&points);
113
114    TClonesArray &cl=*clusters;
115    Int_t nclusters=0;
116    Int_t nentr=(Int_t)pTree->GetEntries();
117
118    cerr<<"Number of entries: "<<nentr<<endl;
119
120    Float_t kmip; // ADC->mip normalization factor for the SDD and SSD 
121
122    for (Int_t i=0; i<nentr; i++) {
123        points->Clear();
124        pTree->GetEvent(i);
125        Int_t ncl=points->GetEntriesFast(); if (ncl==0){cTree->Fill();continue;}
126        Int_t lay,lad,det; geom->GetModuleId(i,lay,lad,det);
127        Float_t x,y,zshift; geom->GetTrans(lay,lad,det,x,y,zshift); 
128        Double_t rot[9];    geom->GetRotMatrix(lay,lad,det,rot);
129        Double_t yshift = x*rot[0] + y*rot[1];
130        Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1);
131        nclusters+=ncl;
132
133        kmip=1.;
134        if(lay<5&&lay>2){kmip=280.;}; // b.b.
135        if(lay<7&&lay>4){kmip=38.;};
136        //cout<<"i,ncl ="<<i<<","<<ncl<<endl;
137
138        for (Int_t j=0; j<ncl; j++) {
139           AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(j);
140           Float_t lp[5];
141           lp[0]=-p->GetX()-yshift; if (lay==1) lp[0]=-lp[0];
142           lp[1]=p->GetZ()+zshift;
143           lp[2]=p->GetSigmaX2();
144           lp[3]=p->GetSigmaZ2();
145           lp[4]=p->GetQ(); lp[4]/=kmip;
146           Int_t lab[6]; 
147           lab[0]=p->GetLabel(0);lab[1]=p->GetLabel(1);lab[2]=p->GetLabel(2);
148           lab[3]=ndet;
149           Int_t label=lab[0];
150
151           //if(label<=0) cout<<" !!!!!!!lab="<<lab[0]<<" "<<endl;   
152
153   if(label>=0) {  // b.b.
154           TParticle *part=(TParticle*)gAlice->Particle(label);
155           label=-3;
156           while (part->P() < 0.005) {
157              Int_t m=part->GetFirstMother();
158              if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;}
159              label=m;
160              part=(TParticle*)gAlice->Particle(label);
161           }
162   }
163           if      (lab[1]<0) lab[1]=label;
164           else if (lab[2]<0) lab[2]=label;
165           else cerr<<"No empty labels !\n";
166
167           new(cl[j]) AliITSclusterV2(lab,lp);
168        }
169        cTree->Fill(); clusters->Delete();
170        points->Delete();
171    }
172    cTree->Write();
173
174    cerr<<"Number of clusters: "<<nclusters<<endl;
175
176    delete cTree; delete clusters; delete points;
177
178    delete gAlice; gAlice=0;
179
180    in->Close();
181    out->Close();
182
183    return 0;
184
185 }
186
187
188
189
190
191
192
193
194
195
196
197
198
199