Fixed grep syntax on OS X in the Analysis Plugin
[u/mrichter/AliRoot.git] / ANALYSIS / AliPPVsMultUtils.cxx
1 /*****************************************************************
2  
3  General utility class for doing pp vs multiplicity studies
4  
5  --- Functionality being added on demand.
6  
7  --- Please report any bugs, comments, suggestions to:
8  david.dobrigkeit.chinellato@cern.ch
9  
10  *****************************************************************/
11
12 #include "AliVEvent.h"
13 #include "AliESDEvent.h"
14 #include "AliAODEvent.h"
15 #include "AliVVertex.h"
16 #include "AliLog.h"
17 #include "AliAODVertex.h"
18 #include "AliVTrack.h"
19 #include "AliVEvent.h"
20 #include <TMatrixDSym.h>
21 #include <TMath.h>
22 #include "AliESDUtils.h"
23 #include "AliPPVsMultUtils.h"
24 #include <TFile.h>
25
26
27 ClassImp(AliPPVsMultUtils)
28
29 //______________________________________________________________________
30 AliPPVsMultUtils::AliPPVsMultUtils():TObject(),
31 fRunNumber(0),
32 fCalibrationLoaded(0),
33 fBoundaryHisto_V0M(0),
34 fBoundaryHisto_V0A(0),
35 fBoundaryHisto_V0C(0),
36 fBoundaryHisto_V0MEq(0),
37 fBoundaryHisto_V0AEq(0),
38 fBoundaryHisto_V0CEq(0)
39 {
40     // Default contructor
41 }
42
43
44 //_____________________________________________________________________________
45 AliPPVsMultUtils::AliPPVsMultUtils(const AliPPVsMultUtils &c) : TObject(c),
46 fRunNumber(0),
47 fCalibrationLoaded(0),
48 fBoundaryHisto_V0M(0),
49 fBoundaryHisto_V0A(0),
50 fBoundaryHisto_V0C(0),
51 fBoundaryHisto_V0MEq(0),
52 fBoundaryHisto_V0AEq(0),
53 fBoundaryHisto_V0CEq(0)
54 {
55     //
56     // copy constructor - untested
57     //
58     ((AliPPVsMultUtils &) c).Copy(*this);
59 }
60
61 //_____________________________________________________________________________
62 AliPPVsMultUtils &AliPPVsMultUtils::operator=(const AliPPVsMultUtils &c)
63 {
64     //
65     // Assignment operator - untested
66     //
67     
68     if (this != &c) ((AliPPVsMultUtils &) c).Copy(*this);
69     return *this;
70 }
71
72 //______________________________________________________________________
73 Float_t AliPPVsMultUtils::GetMultiplicityPercentile(AliVEvent *event, TString lMethod)
74 // Function to get multiplicity quantiles 
75 {
76     Int_t lRequestedRunNumber = event->GetRunNumber();
77     if( lRequestedRunNumber != fRunNumber ) fCalibrationLoaded = LoadCalibration( lRequestedRunNumber );
78     
79     if ( !fCalibrationLoaded ){ 
80       return -1000; //Return absurd value (hopefully noone will use this...) 
81     }
82     
83     Float_t lreturnval = -1;
84     
85     //Get VZERO Information for multiplicity later
86     AliVVZERO* esdV0 = event->GetVZEROData();
87     if (!esdV0) {
88         AliError("AliVVZERO not available");
89         return -1;
90     }
91     
92     // VZERO PART
93     Float_t  multV0A  = 0;            //  multiplicity from V0 reco side A
94     Float_t  multV0C  = 0;            //  multiplicity from V0 reco side C
95     Float_t  multV0AEq  = 0;          //  multiplicity from V0 reco side A
96     Float_t  multV0CEq  = 0;          //  multiplicity from V0 reco side C
97     Float_t  multV0ACorr  = 0;            //  multiplicity from V0 reco side A
98     Float_t  multV0CCorr  = 0;            //  multiplicity from V0 reco side C
99     
100     //Non-Equalized Signal: copy of multV0ACorr and multV0CCorr from AliCentralitySelectionTask
101     //Getters for uncorrected multiplicity
102     multV0A=esdV0->GetMTotV0A();
103     multV0C=esdV0->GetMTotV0C();
104    
105     //Getting around to the SPD vertex -> typecast to ESD/AOD
106     const AliVVertex *lPrimarySPDVtx = NULL;
107     /* get ESD vertex SPD */
108     if (event->InheritsFrom("AliESDEvent")) {
109       AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(event);
110       if (!esdevent) return kFALSE;
111       lPrimarySPDVtx = esdevent->GetPrimaryVertexSPD();
112     }
113     /* get AOD vertex SPD */
114     else if (event->InheritsFrom("AliAODEvent")) {
115       AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(event);
116       if (!aodevent) return kFALSE;
117       lPrimarySPDVtx = aodevent->GetPrimaryVertexSPD();
118     }
119    
120     //Get Z vertex position of SPD vertex (why not Tracking if available?)
121     Float_t zvtx = lPrimarySPDVtx->GetZ();
122     
123     //Acquire Corrected multV0A
124     multV0ACorr = AliESDUtils::GetCorrV0A(multV0A,zvtx);
125     multV0CCorr = AliESDUtils::GetCorrV0C(multV0C,zvtx);
126     
127     // Equalized signals // From AliCentralitySelectionTask // Updated
128     for(Int_t iCh = 32; iCh < 64; ++iCh) {
129         Double_t mult = event->GetVZEROEqMultiplicity(iCh);
130         multV0AEq += mult;
131     }
132     for(Int_t iCh = 0; iCh < 32; ++iCh) {
133         Double_t mult = event->GetVZEROEqMultiplicity(iCh);
134         multV0CEq += mult;
135     }
136     
137     if ( lMethod == "V0M" ) lreturnval =
138         fBoundaryHisto_V0M -> GetBinContent( fBoundaryHisto_V0M->FindBin(multV0ACorr+multV0CCorr) );
139     if ( lMethod == "V0A" ) lreturnval =
140         fBoundaryHisto_V0A -> GetBinContent( fBoundaryHisto_V0A->FindBin(multV0ACorr) );
141     if ( lMethod == "V0C" ) lreturnval =
142         fBoundaryHisto_V0C -> GetBinContent( fBoundaryHisto_V0C->FindBin(multV0CCorr) );
143     //equalized
144     if ( lMethod == "V0MEq" ) lreturnval =
145         fBoundaryHisto_V0MEq -> GetBinContent( fBoundaryHisto_V0MEq->FindBin(multV0AEq+multV0CEq) );
146     if ( lMethod == "V0AEq" ) lreturnval =
147         fBoundaryHisto_V0AEq -> GetBinContent( fBoundaryHisto_V0AEq->FindBin(multV0AEq) );
148     if ( lMethod == "V0CEq" ) lreturnval =
149         fBoundaryHisto_V0CEq -> GetBinContent( fBoundaryHisto_V0CEq->FindBin(multV0CEq) );
150     
151     return lreturnval;
152 }
153
154 //______________________________________________________________________
155 Bool_t AliPPVsMultUtils::LoadCalibration(Int_t lLoadThisCalibration)
156 //To be called if starting analysis on a new run
157 {
158     //If Histograms exist, de-allocate to prevent memory leakage
159     if( fBoundaryHisto_V0M ) {fBoundaryHisto_V0M->Delete(); fBoundaryHisto_V0M = 0x0; }
160     if( fBoundaryHisto_V0A ) {fBoundaryHisto_V0A->Delete(); fBoundaryHisto_V0A = 0x0; }
161     if( fBoundaryHisto_V0C ) {fBoundaryHisto_V0C->Delete(); fBoundaryHisto_V0C = 0x0; }
162     if( fBoundaryHisto_V0MEq ) {fBoundaryHisto_V0MEq->Delete(); fBoundaryHisto_V0MEq = 0x0; }
163     if( fBoundaryHisto_V0AEq ) {fBoundaryHisto_V0AEq->Delete(); fBoundaryHisto_V0AEq = 0x0; }
164     if( fBoundaryHisto_V0CEq ) {fBoundaryHisto_V0CEq->Delete(); fBoundaryHisto_V0CEq = 0x0; }
165     
166     AliInfo(Form( "Loading calibration file for run %i",lLoadThisCalibration) );
167     TFile *lCalibFile_V0M = 0x0;
168     TFile *lCalibFile_V0A = 0x0;
169     TFile *lCalibFile_V0C = 0x0;
170     TFile *lCalibFile_V0MEq = 0x0;
171     TFile *lCalibFile_V0AEq = 0x0;
172     TFile *lCalibFile_V0CEq = 0x0;
173     
174     //AliInfo("Calling TFile::Open");
175     lCalibFile_V0M = TFile::Open("$ALICE_ROOT/PWGLF/STRANGENESS/Cascades/corrections/calibration_V0M.root");
176     lCalibFile_V0A = TFile::Open("$ALICE_ROOT/PWGLF/STRANGENESS/Cascades/corrections/calibration_V0A.root");
177     lCalibFile_V0C = TFile::Open("$ALICE_ROOT/PWGLF/STRANGENESS/Cascades/corrections/calibration_V0C.root");
178     lCalibFile_V0MEq = TFile::Open("$ALICE_ROOT/PWGLF/STRANGENESS/Cascades/corrections/calibration_V0MEq.root");
179     lCalibFile_V0AEq = TFile::Open("$ALICE_ROOT/PWGLF/STRANGENESS/Cascades/corrections/calibration_V0AEq.root");
180     lCalibFile_V0CEq = TFile::Open("$ALICE_ROOT/PWGLF/STRANGENESS/Cascades/corrections/calibration_V0CEq.root");
181
182     //AliInfo("Casting");
183     //check memory consumption later... 
184     fBoundaryHisto_V0M   = dynamic_cast<TH1F *>(lCalibFile_V0M  -> Get(Form("histocalib%i",lLoadThisCalibration)) );
185     fBoundaryHisto_V0A   = dynamic_cast<TH1F *>(lCalibFile_V0A  -> Get(Form("histocalib%i",lLoadThisCalibration)) );
186     fBoundaryHisto_V0C   = dynamic_cast<TH1F *>(lCalibFile_V0C  -> Get(Form("histocalib%i",lLoadThisCalibration)) );
187     fBoundaryHisto_V0MEq   = dynamic_cast<TH1F *>(lCalibFile_V0MEq  -> Get(Form("histocalib%i",lLoadThisCalibration)) );
188     fBoundaryHisto_V0AEq   = dynamic_cast<TH1F *>(lCalibFile_V0AEq  -> Get(Form("histocalib%i",lLoadThisCalibration)) );
189     fBoundaryHisto_V0CEq   = dynamic_cast<TH1F *>(lCalibFile_V0CEq  -> Get(Form("histocalib%i",lLoadThisCalibration)) );
190     
191     if ( !fBoundaryHisto_V0M   || !fBoundaryHisto_V0A   || !fBoundaryHisto_V0C ||
192         !fBoundaryHisto_V0MEq || !fBoundaryHisto_V0AEq || !fBoundaryHisto_V0CEq ){
193         AliInfo(Form("No calibration for run %i exists at the moment!",lLoadThisCalibration));
194         return kFALSE; //return denial
195     }
196     
197     //Careful with manual cleanup if needed: to be implemented
198     fBoundaryHisto_V0M->SetDirectory(0);
199     fBoundaryHisto_V0A->SetDirectory(0);
200     fBoundaryHisto_V0C->SetDirectory(0);
201     fBoundaryHisto_V0MEq->SetDirectory(0);
202     fBoundaryHisto_V0AEq->SetDirectory(0);
203     fBoundaryHisto_V0CEq->SetDirectory(0);
204     
205     //AliInfo("Closing");
206     
207     if( lCalibFile_V0M ) lCalibFile_V0M->Close();
208     if( lCalibFile_V0A ) lCalibFile_V0A->Close();
209     if( lCalibFile_V0C ) lCalibFile_V0C->Close();
210     if( lCalibFile_V0MEq ) lCalibFile_V0MEq->Close();
211     if( lCalibFile_V0AEq ) lCalibFile_V0AEq->Close();
212     if( lCalibFile_V0CEq ) lCalibFile_V0CEq->Close();
213     
214     fRunNumber = lLoadThisCalibration; //Loaded!
215     AliInfo(Form("Finished loading calibration for run %i",lLoadThisCalibration));
216     return kTRUE;
217 }