Transition to NewIO
[u/mrichter/AliRoot.git] / FASTSIM / AliFastMuonTriggerEff.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.2.2.1  2003/04/15 15:57:31  hristov
19 Merging with v3-09-08
20
21 Revision 1.2  2003/01/21 10:49:17  morsch
22 Corrected path for vettotpara.root
23
24 Revision 1.1  2003/01/06 10:13:09  morsch
25 First commit.
26
27 */
28
29
30 #include "AliFastMuonTriggerEff.h"
31 #include <TFile.h>
32 #include <TTree.h>
33 #include <TROOT.h>
34 #include <stdlib.h>
35
36 // Debugging flag
37 //#define MYTRIGDEBUG
38
39 ClassImp(AliFastMuonTriggerEff);
40
41 AliFastMuonTriggerEff::AliFastMuonTriggerEff():
42     AliFastResponse("Efficiency", "Muon Trigger Efficiency")
43 {
44 //
45 // Default constructor
46 //
47     SetCut(kLow);
48     fZones=0;
49 }
50
51
52 void AliFastMuonTriggerEff::Init()
53 {
54 //
55 //  Initialization
56 //
57     printf("Initializing %s / %s", fName.Data(), fTitle.Data());
58 //
59
60     fPhiMin   = -90. ;
61     fPhiMax   =  90. ;
62     fDphi     =   9. ;
63     fThetaMin =   2. ;
64     fThetaMax =   9. ;
65     fDtheta   =   0.7;
66     fDpt      =   0.1;
67
68     InitTree();
69 }
70
71 void AliFastMuonTriggerEff::InitTree()
72 {
73 //
74 //  Initialize tables from tree
75 //
76     TTree          *chain;   // pointer to the analyzed TTree or TChain
77     Int_t           nSim, nZona;
78     Double_t        pmin, pmax, tmin, tmax, phimin, phimax, bkg;
79     Double_t        len50, hen50, leff, heff;
80     Double_t        vLPt[50];
81     Double_t        vHPt[50];
82     Char_t file[100]="$(ALICE_ROOT)/FASTSIM/data/vettorpara.root";
83 //
84 //  Avoid memory leak in case of reinitialization
85     if(fZones!=0) {
86         printf("\nWarning: reinitialization of an object of class: AliFastMuonTriggerEff\n");
87         for (Int_t i=0; i<fZones; i++) {
88             if(fEffLow [i])delete[] fEffLow[i];
89             if(fEffHigh[i])delete[] fEffHigh[i];
90         }
91         if(fEffLow) {
92             delete[] fEffLow;
93             fEffLow=0;
94         }
95         if(fEffHigh) {
96             delete[] fEffHigh;
97             fEffHigh=0;
98         }
99     }
100     printf("AliFastMuonTriggerEff: Initialization\n");
101     TFile *f = new TFile(file);
102     if(f->IsZombie()) {
103         printf("Cannot open file: %s\n",file);
104         return;
105     }
106     f->ls();
107     
108     TTree* tree = (TTree*)gDirectory->Get("fitPar");
109
110 //
111 //
112     chain = tree;
113     chain->SetMakeClass(1);
114     
115     chain->SetBranchAddress("nSim",&nSim);
116     chain->SetBranchAddress("nZona",&nZona);
117     chain->SetBranchAddress("ptmin",&pmin);
118     chain->SetBranchAddress("ptmax",&pmax);
119     chain->SetBranchAddress("Thetamin",&tmin);
120     chain->SetBranchAddress("Thetamax",&tmax);
121     chain->SetBranchAddress("Phimin",&phimin);
122     chain->SetBranchAddress("Phimax",&phimax);
123     chain->SetBranchAddress("Bkg",&bkg);
124     chain->SetBranchAddress("EffLPt",vLPt);
125     chain->SetBranchAddress("EffHPt",vHPt);
126     chain->SetBranchAddress("Pt0.5Low",&len50);
127     chain->SetBranchAddress("Pt0.5High",&hen50);
128     chain->SetBranchAddress("EffLowMax",&leff);
129     chain->SetBranchAddress("EffHighMax",&heff);
130 // 
131 //
132     Int_t nentries = Int_t(chain->GetEntries());
133     Int_t nbytes = 0, nb = 0;
134
135 // Count the number of zones of the parametrization
136     Int_t nzone0=0, nzone1=0;
137     for (Int_t jentry=0; jentry<nentries; jentry++) {
138         nb = chain->GetEntry(jentry);
139         if(bkg==0.) {
140              if(nSim==0)nzone0++;
141              if(nSim==1)nzone1++;
142         }
143     }
144     
145     printf("Trigger parametrisation for %d zones for Pt: 0. 5. GeV/c\n",nzone0);
146     printf("and %d zones extended to 10. GeV/c\n",nzone1);
147     fZones=nzone0+nzone1;
148 //    printf("Ciao\n");
149     if(fZones<=0){
150         printf("Number of zones must be greater than 0\n");
151         exit(6);
152     }
153     
154     fEffLow =new Float_t*[fZones];
155     fEffHigh=new Float_t*[fZones];
156     for (Int_t i=0; i<fZones; i++) {
157         fEffLow [i]=new Float_t[50];
158         fEffHigh[i]=new Float_t[50];
159     }
160     
161 //  Initialize look-up table to standard values
162     Int_t isim, itheta, iphi;
163     for (isim=0; isim<2; isim++) {
164         for (itheta=0; itheta<10; itheta++) {
165             for (iphi=0; iphi<20; iphi++) {
166                 fLook[isim][itheta][iphi]=0;
167             }
168         }
169     }
170
171 //  Loading Trigger efficiencies
172     Int_t myzone=0;
173 #ifdef MYTRIGDEBUG
174             printf("Filling nSim nZona pmin pmax tmin tmax phimin phimax: ....\n");
175 #endif
176     for (Int_t jentry=0; jentry<nentries; jentry++) {
177 //      Int_t ientry = LoadTree(jentry); 
178         nb = chain->GetEntry(jentry);   
179         nbytes += nb;
180 #ifdef MYTRIGDEBUG
181        printf("Getting entry %d... ",jentry);
182 #endif
183 // For the time being it works with background 0
184         if ((nSim == 0 || nSim == 1)&&bkg==0.) {
185 #ifdef MYTRIGDEBUG
186             printf("Filling %d %d %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f: ",
187             nSim,nZona,pmin,pmax,tmin,tmax,phimin,phimax);
188 #endif
189             for (Int_t k = 0; k < 50; k++) {
190                 fEffLow [myzone][k] = vLPt[k];
191                 fEffHigh[myzone][k] = vHPt[k];
192 #ifdef MYTRIGDEBUG
193                 if(k<15)printf(" %5.3f",vLPt[k]);
194 #endif
195             }
196 #ifdef MYTRIGDEBUG
197             printf("\n");
198 #endif
199             myzone++;
200             iphi=Int_t(nZona/10)-10;
201             itheta=nZona%10;
202             if(iphi<0||iphi>19||itheta<0||itheta>9) {
203                 printf("The format of data file is not consistent\nwith this version of the code\n");
204                 printf("This should never happen: iphi %d, itheta: %d\n",iphi,itheta);
205                 exit(7);
206             }
207             fLook[nSim][itheta][iphi]=myzone;
208         } else {
209             printf("Skipping entry with nSim=%d, bkg=%f\n",nSim,bkg);
210         }
211     }
212 #ifdef MYTRIGDEBUG
213     printf("This is the content of the LUT after first step of initialization\n");
214     for(isim=0; isim<2; isim++) {
215         printf("isim=%d\n",isim);
216         for(iphi=0; iphi<20; iphi++) {
217             printf("iphi: %2d:",iphi);
218             for(itheta=0; itheta<10; itheta++) {
219                 printf(" %4d",fLook[isim][itheta][iphi]);
220             }
221             printf("\n");
222         }
223     }
224 #endif
225 // Filling look=up table for the zones where the extended simulation does
226 // not exists
227     for(iphi=0; iphi<20; iphi++) {
228         for(itheta=0; itheta<10; itheta++) {
229             if(fLook[0][itheta][iphi]==0) {
230                 printf("Missing entry isim=%d itheta=%d iphi=%d\n",isim,itheta,iphi);
231                 exit(8);
232             }
233             if(fLook[0][itheta][iphi]<0||fLook[0][itheta][iphi]>fZones) {
234                 printf("Problem with entry isim=%d itheta=%d iphi=%d\n",isim,itheta,iphi);
235                 exit(9);
236             }
237             if(fLook[1][itheta][iphi]==0) {
238                  fLook[1][itheta][iphi]=-fLook[0][itheta][iphi];
239             }
240         }
241     }
242 #ifdef MYTRIGDEBUG
243     for(isim=0; isim<2; isim++) {
244         printf("isim=%d\n",isim);
245         for(iphi=0; iphi<20; iphi++) {
246             printf("iphi: %2d:",iphi);
247             for(itheta=0; itheta<10; itheta++) {
248                 printf(" %4d",fLook[isim][itheta][iphi]);
249             }
250             printf("\n");
251         }
252     }
253         for(iphi=0; iphi<20; iphi++) {
254             for(itheta=0; itheta<10; itheta++) {
255     for(isim=0; isim<2; isim++) {
256                 printf("%d %2d %2d %4d:",isim,iphi,itheta,fLook[isim][itheta][iphi]);
257                 if(fLook[isim][itheta][iphi]>0) {
258                     myzone=fLook[isim][itheta][iphi]-1;
259                     for(Int_t ipt=0; ipt<20; ipt++) {
260                         //printf(" %5.3f",fEffLow[myzone][ipt]);
261                         printf(" %5.3f",fEffHigh[myzone][ipt]);
262                     }
263                     printf(" ...");
264                     for(Int_t ipt=40; ipt<50; ipt++) {
265                         //printf(" %5.3f",fEffLow[myzone][ipt]);
266                         printf(" %5.3f",fEffHigh[myzone][ipt]);
267                     }               
268                 }
269                 printf("\n");
270             }
271         }
272     }    
273 #endif
274     f->Close();
275 }
276
277 void AliFastMuonTriggerEff::Evaluate(Float_t charge, Float_t pt,
278                 Float_t theta, Float_t phi,Float_t& effLow, Float_t& effHigh)
279 {
280 //
281 //  Trigger efficiency for pt, theta, phi (low and high cut)
282 //
283     effLow=0.;
284     effHigh=0.;
285     if(fZones==0) {
286         printf("Call to uninitialized object of class: AliFastMuonTriggerEff\n");
287         return;
288     }
289     if(pt<0) {
290         printf("Warning: pt: %f < 0. GeV/c\n",pt);
291         return; 
292     }
293     Int_t iPt   = Int_t(pt/fDpt)%50;
294     Int_t iSim  = Int_t(pt/fDpt)/50;
295     Int_t iPhi  = Int_t((phi-fPhiMin)/fDphi);
296     if(phi<fPhiMin)iPhi=iPhi-1;
297     Int_t iTheta = Int_t((theta-fThetaMin)/fDtheta);
298 #ifdef MYTRIGDEBUG
299     printf("iSim iPt iTheta iPhi: %d %d %d %d\n",iSim,iPt,iTheta,iPhi);
300 #endif
301     iPhi=iPhi-40*(iPhi/40);
302 #ifdef MYTRIGDEBUG
303     printf("1: iPhi converted to: %d for angle equivalence\n",iPhi);
304 #endif
305     if(iPhi<0)iPhi=-iPhi-1;
306     if(iPhi>19)iPhi=39-iPhi;
307 #ifdef MYTRIGDEBUG
308     printf("2: iPhi converted to: %d for the symmetry of the spectrometer\n",iPhi);
309 #endif
310     if(charge==1.){
311     } else if(charge==-1.) {
312     iPhi=19-iPhi;
313 #ifdef MYTRIGDEBUG
314     printf("3: iPhi converted to: %d for charge symmetry\n",iPhi);
315 #endif
316     } else {
317         printf("Warning: not understand charge: %f\n",charge);
318         return;
319     }
320     if(iTheta<0||iTheta>9) {
321         printf("Warning: theta: %f outside acceptance\n",theta);
322         return;
323     }
324     if(iPt<0) {
325         printf("Warning: what do you mean with pt: %f <0?\n",pt);
326         return;
327     }
328     if(iSim>=fSim) {
329         iSim=fSim-1;
330         iPt=49;
331 #ifdef MYTRIGDEBUG
332     printf("4: iSim iPt converted to: %d %d (last zone)\n",iSim,iPt);
333 #endif
334     }
335     Int_t iLook=fLook[iSim][iTheta][iPhi];
336     if(iLook<0) {
337 #ifdef MYTRIGDEBUG
338     printf("5: iLook iPt: %d %d converted to: ",iLook,iPt);
339 #endif
340         iLook=-iLook-1;
341         iPt=49;
342 #ifdef MYTRIGDEBUG
343     printf("%d %d from look up table contents\n",iLook,iPt);
344 #endif
345     } else {
346         iLook=iLook-1;
347     }
348     effLow=fEffLow[iLook][iPt];
349     effHigh=fEffHigh[iLook][iPt];
350 #ifdef MYTRIGDEBUG
351     printf("6: charge, iSim, iTheta, iPhi, iPt: %f %d %d %d %d effLow: %f, effHigh: %f\n",
352                charge,iSim,iTheta,iPhi,iPt,effLow,effHigh);
353 #endif
354     
355     //fEffLow [iPhi][iTheta][iPt];
356     //fEffHigh[iPhi][iTheta][iPt];    
357 }
358
359 Float_t AliFastMuonTriggerEff::Evaluate(Float_t charge, Float_t pt,
360                    Float_t theta, Float_t phi)
361 {
362     //
363     // Trigger efficiency for pt, theta, phi depending of fCut
364     // 
365     if(fZones==0) {
366         printf("Call to uninitialized object of class: AliFastMuonTriggerEff\n");
367         return 0.;
368     }
369     Float_t eff;
370     Float_t effLow, effHigh;
371     
372     Evaluate(charge,pt,theta,phi,effLow,effHigh);
373     if (fCut == kLow) 
374         eff  = effLow;
375     else
376         eff  = effHigh;
377
378     return eff;
379 }