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