2151dbf1fd26d6a23fe5bbbae0833c4a6c072c39
[u/mrichter/AliRoot.git] / PWG4 / totEt / AliAnalysisHadEtReconstructed.cxx
1 //_________________________________________________________________________
2 //  Utility Class for transverse energy studies, charged hadrons
3 //  Base class for ESD analysis
4 //  - reconstruction output
5 // implementation file
6 //
7 //Created by Christine Nattrass, Rebecca Scott, Irakli Martashvili
8 //University of Tennessee at Knoxville
9 //_________________________________________________________________________
10
11 #include <TROOT.h>
12 #include <TSystem.h>
13 #include <TInterpreter.h>
14 #include "AliAnalysisHadEtReconstructed.h"
15 #include "AliAnalysisEtCuts.h"
16 #include "AliESDtrack.h"
17 #include "AliESDCaloCluster.h"
18 #include "AliVEvent.h"
19 #include "AliESDEvent.h"
20 #include "AliESDtrackCuts.h"
21 #include "AliESDpid.h"
22 #include "AliVParticle.h"
23 #include <iostream>
24 #include "AliAnalysisHadEtCorrections.h"
25 #include "TFile.h"
26 #include "TString.h"
27
28 using namespace std;
29
30 ClassImp(AliAnalysisHadEtReconstructed);
31
32
33 AliAnalysisHadEtReconstructed::AliAnalysisHadEtReconstructed() :
34         AliAnalysisHadEt()
35         ,corrections(0)
36         ,fConfigFile("ConfigHadEtAnalysis.C")
37         ,fCorrTotEtFullAcceptanceTPC(0)
38         ,fCorrTotEtFullAcceptanceITS(0)
39         ,fCorrHadEtFullAcceptanceTPC(0)
40         ,fCorrHadEtFullAcceptanceITS(0)
41         ,fCorrTotEtEMCALAcceptanceTPC(0)
42         ,fCorrTotEtEMCALAcceptanceITS(0)
43         ,fCorrHadEtEMCALAcceptanceTPC(0)
44         ,fCorrHadEtEMCALAcceptanceITS(0)
45         ,fCorrTotEtPHOSAcceptanceTPC(0)
46         ,fCorrTotEtPHOSAcceptanceITS(0)
47         ,fCorrHadEtPHOSAcceptanceTPC(0)
48         ,fCorrHadEtPHOSAcceptanceITS(0)
49         ,fRawEtFullAcceptanceTPC(0)
50         ,fRawEtFullAcceptanceITS(0)
51         ,fRawEtEMCALAcceptanceTPC(0)
52         ,fRawEtEMCALAcceptanceITS(0)
53         ,fRawEtPHOSAcceptanceTPC(0)
54         ,fRawEtPHOSAcceptanceITS(0)
55 {
56
57 }
58
59 AliAnalysisHadEtReconstructed::~AliAnalysisHadEtReconstructed() 
60 {
61 }
62
63 Int_t AliAnalysisHadEtReconstructed::AnalyseEvent(AliVEvent* ev)
64 { // analyse ESD event
65     ResetEventValues();
66     fRawEtFullAcceptanceTPC=0.0;
67     fRawEtFullAcceptanceITS=0.0;
68     fRawEtEMCALAcceptanceTPC=0.0;
69     fRawEtEMCALAcceptanceITS=0.0;
70     fRawEtPHOSAcceptanceTPC=0.0;
71     fRawEtPHOSAcceptanceITS=0.0;
72
73     AliESDEvent *realEvent = dynamic_cast<AliESDEvent*>(ev);
74     //for PID
75     AliESDpid *pID = new AliESDpid();
76     pID->MakePID(realEvent);
77
78     TString *strTPC = new TString("TPC");
79     TString *strITS = new TString("ITS");
80     TString *strTPCITS = new TString("TPCITS");
81     for(Int_t cutset=0;cutset<2;cutset++){
82       TString *cutName;
83       TObjArray* list;
84       switch(cutset){
85       case 0:
86         cutName = strTPC;
87         list = fEsdtrackCutsTPC->GetAcceptedTracks(realEvent);
88         break;
89       case 1:
90         cutName = strITS;
91         list = fEsdtrackCutsITS->GetAcceptedTracks(realEvent);
92         break;
93       case 2:
94         cutName = strTPCITS;
95         list = fEsdtrackCutsITSTPC->GetAcceptedTracks(realEvent);
96         break;
97       default:
98         cerr<<"Error:  cannot fill histograms!"<<endl;
99         return -1;
100       }
101       Int_t nGoodTracks = list->GetEntries();
102       for (Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++)
103         {
104
105
106           AliESDtrack *track = dynamic_cast<AliESDtrack*> (list->At(iTrack));
107           if (!track)
108             {
109               Printf("ERROR: Could not get track %d", iTrack);
110               continue;
111             }
112           else{
113             Float_t nSigmaPion,nSigmaProton,nSigmaKaon,nSigmaElectron;
114             
115             if(cutset!=1){
116               nSigmaPion = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kPion));
117               nSigmaProton = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kProton));
118               nSigmaKaon = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kKaon));
119               nSigmaElectron = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kElectron));
120             }
121             else{
122               nSigmaPion = TMath::Abs(pID->NumberOfSigmasITS(track,AliPID::kPion));
123               nSigmaProton = TMath::Abs(pID->NumberOfSigmasITS(track,AliPID::kProton));
124               nSigmaKaon = TMath::Abs(pID->NumberOfSigmasITS(track,AliPID::kKaon));
125               nSigmaElectron = TMath::Abs(pID->NumberOfSigmasITS(track,AliPID::kElectron));
126             }
127             bool isPion = (nSigmaPion<3.0 && nSigmaProton>2.0 && nSigmaKaon>2.0);
128             bool isElectron = (nSigmaElectron<2.0 && nSigmaPion>4.0 && nSigmaProton>3.0 && nSigmaKaon>3.0);
129             bool isKaon = (nSigmaPion>3.0 && nSigmaProton>2.0 && nSigmaKaon<2.0);
130             bool isProton = (nSigmaPion>3.0 && nSigmaProton<2.0 && nSigmaKaon>2.0);
131
132             //bool IsElectron = false;
133             bool unidentified = (!isProton && !isKaon && !isElectron);
134             Float_t dEdx = track->GetTPCsignal();
135             if(cutset==1) dEdx = track->GetITSsignal();
136             FillHisto2D(Form("dEdxDataAll%s",cutName->Data()),track->P(),dEdx,1.0);
137
138             //if(!(corrections->GetEfficiencyPionTPC())) cerr<<"Uh-oh!  No histogram!"<<endl;
139
140             Float_t corrBkgd=0.0;
141             Float_t corrNotID=0.0;
142             Float_t corrNoID = corrections->GetNotIDCorrectionNoPID(track->Pt());
143             Float_t corrEff = 0.0;
144             Float_t corrEffNoID = 0.0;
145             if(cutset==0){//TPC
146               corrBkgd = corrections->GetBackgroundCorrectionTPC(track->Pt());
147               corrEffNoID = corrections->GetTPCEfficiencyCorrectionHadron(track->Pt());
148               corrNotID = corrections->GetNotIDCorrectionTPC(track->Pt());
149             }
150             if(cutset==1){//ITS
151               corrBkgd = corrections->GetBackgroundCorrectionITS(track->Pt());
152               //corrEffNoID = corrections->GetITSEfficiencyCorrectionHadron(track->Pt());
153               corrNotID = corrections->GetNotIDCorrectionITS(track->Pt());
154             }
155             Float_t et = 0.0;
156             Float_t etpartialcorrected = 0.0;
157             Float_t etpartialcorrectedNoID = corrNoID*corrBkgd*corrEffNoID*Et(track->P(),track->Theta(),fPiPlusCode,track->Charge());
158             FillHisto2D(Form("EtDataRaw%sNoID",cutName->Data()),track->Pt(),track->Eta(),etpartialcorrectedNoID);
159
160             if(isPion){
161               FillHisto2D(Form("dEdxDataPion%s",cutName->Data()),track->P(),dEdx,1.0);
162               et = Et(track->P(),track->Theta(),fPiPlusCode,track->Charge());
163               if(cutset==0){corrEff = corrections->GetTPCEfficiencyCorrectionPion(track->Pt());}
164               //else{corrEff = corrections->GetITSEfficiencyCorrectionPion(track->Pt());}
165               etpartialcorrected = et*corrBkgd*corrEff;
166               
167               if(track->Charge()>0.0){
168                 FillHisto2D(Form("EtDataRaw%sPiPlus",cutName->Data()),track->Pt(),track->Eta(),et);
169                 FillHisto2D(Form("EtDataCorrected%sPiPlus",cutName->Data()),track->Pt(),track->Eta(),etpartialcorrected);
170               }
171               else{
172                 FillHisto2D(Form("EtDataRaw%sPiMinus",cutName->Data()),track->Pt(),track->Eta(),et);
173                 FillHisto2D(Form("EtDataCorrected%sPiMinus",cutName->Data()),track->Pt(),track->Eta(),etpartialcorrected);
174               }
175             }
176             if(isKaon){
177               FillHisto2D(Form("dEdxDataKaon%s",cutName->Data()),track->P(),dEdx,1.0);
178               et = Et(track->P(),track->Theta(),fKPlusCode,track->Charge());
179               if(cutset==0){corrEff = corrections->GetTPCEfficiencyCorrectionKaon(track->Pt());}
180               //else{corrEff = corrections->GetITSEfficiencyCorrectionKaon(track->Pt());}
181               etpartialcorrected = et*corrBkgd*corrEff;
182               
183               if(track->Charge()>0.0){
184                 FillHisto2D(Form("EtDataRaw%sKPlus",cutName->Data()),track->Pt(),track->Eta(),et);
185                 FillHisto2D(Form("EtDataCorrected%sKPlus",cutName->Data()),track->Pt(),track->Eta(),etpartialcorrected);
186               }
187               else{
188                 FillHisto2D(Form("EtDataRaw%sKMinus",cutName->Data()),track->Pt(),track->Eta(),et);
189                 FillHisto2D(Form("EtDataCorrected%sKMinus",cutName->Data()),track->Pt(),track->Eta(),etpartialcorrected);
190               }
191             }
192             if(isProton){
193               FillHisto2D(Form("dEdxDataProton%s",cutName->Data()),track->P(),dEdx,1.0);
194               et = Et(track->P(),track->Theta(),fProtonCode,track->Charge());
195               if(cutset==0){corrEff = corrections->GetTPCEfficiencyCorrectionProton(track->Pt());}
196               //else{corrEff = corrections->GetITSEfficiencyCorrectionProton(track->Pt());}
197               etpartialcorrected = et*corrBkgd*corrEff;
198               
199               if(track->Charge()>0.0){
200                 FillHisto2D(Form("EtDataRaw%sProton",cutName->Data()),track->Pt(),track->Eta(),et);
201                 FillHisto2D(Form("EtDataCorrected%sProton",cutName->Data()),track->Pt(),track->Eta(),etpartialcorrected);
202               }
203               else{
204                 FillHisto2D(Form("EtDataRaw%sAntiProton",cutName->Data()),track->Pt(),track->Eta(),et);
205                 FillHisto2D(Form("EtDataCorrected%sAntiProton",cutName->Data()),track->Pt(),track->Eta(),etpartialcorrected);
206               }
207             }
208             if(isElectron){
209               FillHisto2D(Form("dEdxDataProton%s",cutName->Data()),track->P(),dEdx,1.0);
210               //et = Et(track->P(),track->Theta(),fPiPlusCode,track->Charge());
211             }
212             if(unidentified){
213               FillHisto2D(Form("dEdxDataUnidentified%s",cutName->Data()),track->P(),dEdx,1.0);
214               et = Et(track->P(),track->Theta(),fPiPlusCode,track->Charge());
215               etpartialcorrected = et*corrBkgd*corrEffNoID*corrNotID;
216               FillHisto2D(Form("EtDataCorrected%sUnidentified",cutName->Data()),track->Pt(),track->Eta(),etpartialcorrected);
217             }
218           }
219         }
220     }
221     return 1;
222 }
223
224 bool AliAnalysisHadEtReconstructed::CheckGoodVertex(AliVParticle* track)
225 { // check vertex
226
227     Float_t bxy = 999.;
228     Float_t bz = 999.;
229     dynamic_cast<AliESDtrack*>(track)->GetImpactParametersTPC(bxy,bz);
230
231     bool status = (TMath::Abs(track->Xv()) < fCuts->GetReconstructedVertexXCut()) && 
232       (TMath::Abs(track->Yv()) < fCuts->GetReconstructedVertexYCut()) && 
233       (TMath::Abs(track->Zv()) < fCuts->GetReconstructedVertexZCut()) && 
234       (TMath::Abs(bxy) < fCuts->GetReconstructedIPxyCut()) && 
235       (TMath::Abs(bz) < fCuts->GetReconstructedIPzCut()); 
236
237     return status;
238 }
239
240 void AliAnalysisHadEtReconstructed::Init()
241 { // Init
242     AliAnalysisHadEt::Init();
243
244   if (fConfigFile.Length()) {
245     gROOT->LoadMacro(fConfigFile);
246     corrections = (AliAnalysisHadEtCorrections *) gInterpreter->ProcessLine("ConfigHadEtAnalysis()");
247     fCorrTotEtFullAcceptanceTPC = corrections->GetConstantCorrections(kTRUE,0.15,"Full");
248     fCorrTotEtFullAcceptanceITS = corrections->GetConstantCorrections(kTRUE,0.1,"Full");
249     fCorrHadEtFullAcceptanceTPC = corrections->GetConstantCorrections(kFALSE,0.15,"Full");
250     fCorrHadEtFullAcceptanceITS = corrections->GetConstantCorrections(kFALSE,0.1,"Full");
251     fCorrTotEtEMCALAcceptanceTPC = corrections->GetConstantCorrections(kTRUE,0.15,"EMCAL");
252     fCorrTotEtEMCALAcceptanceITS = corrections->GetConstantCorrections(kTRUE,0.1,"EMCAL");
253     fCorrHadEtEMCALAcceptanceTPC = corrections->GetConstantCorrections(kFALSE,0.15,"EMCAL");
254     fCorrHadEtEMCALAcceptanceITS = corrections->GetConstantCorrections(kFALSE,0.1,"EMCAL");
255     fCorrTotEtPHOSAcceptanceTPC = corrections->GetConstantCorrections(kTRUE,0.15,"PHOS");
256     fCorrTotEtPHOSAcceptanceITS = corrections->GetConstantCorrections(kTRUE,0.1,"PHOS");
257     fCorrHadEtPHOSAcceptanceTPC = corrections->GetConstantCorrections(kFALSE,0.15,"PHOS");
258     fCorrHadEtPHOSAcceptanceITS = corrections->GetConstantCorrections(kFALSE,0.1,"PHOS");
259
260   }
261 }
262
263
264 void AliAnalysisHadEtReconstructed::CreateHistograms(){
265
266   TString *strTPC = new TString("TPC");
267   TString *strITS = new TString("ITS");
268   TString *strTPCITS = new TString("TPCITS");
269   for(Int_t i=0;i<2;i++){
270     TString *cutName;
271     Float_t maxPtdEdx = 10;
272     Float_t mindEdx = 35;
273     Float_t maxdEdx = 150.0;
274     switch(i){
275     case 0:
276       cutName = strTPC;
277       break;
278     case 1:
279       cutName = strITS;
280       maxPtdEdx = 5;
281       maxdEdx = 500.0;
282       break;
283     case 2:
284       cutName = strTPCITS;
285       break;
286     default:
287       cerr<<"Error:  cannot make histograms!"<<endl;
288       return;
289     }
290
291     CreateEtaPtHisto2D(Form("EtDataRaw%sPiPlus",cutName->Data()),"Raw reconstructed E_{T} from identified #pi^{+}");
292     CreateEtaPtHisto2D(Form("EtDataRaw%sPiMinus",cutName->Data()),"Raw reconstructed E_{T} from identified #pi^{-}");
293     CreateEtaPtHisto2D(Form("EtDataRaw%sKPlus",cutName->Data()),"Raw reconstructed E_{T} from identified K^{+}");
294 //     CreateEtaPtHisto2D(Form("EtDataRaw%sEMinus",cutName->Data()),"Raw reconstructed E_{T} from identified e^{-}");
295 //     CreateEtaPtHisto2D(Form("EtDataRaw%sEPlus",cutName->Data()),"Raw reconstructed E_{T} from identified e^{+}");
296     CreateEtaPtHisto2D(Form("EtDataRaw%sKMinus",cutName->Data()),"Raw reconstructed E_{T} from identified K^{-}");
297     CreateEtaPtHisto2D(Form("EtDataRaw%sProton",cutName->Data()),"Raw reconstructed E_{T} from identified p");
298     CreateEtaPtHisto2D(Form("EtDataRaw%sAntiProton",cutName->Data()),"Raw reconstructed E_{T} from identified #bar{p}");
299     CreateEtaPtHisto2D(Form("EtDataRaw%sUnidentified",cutName->Data()),"Raw reconstructed E_{T} from unidentified particles using real mass");
300     CreateEtaPtHisto2D(Form("EtDataRaw%sNoID",cutName->Data()),"Raw reconstructed E_{T} from unidentified particles using real mass");
301
302     CreateEtaPtHisto2D(Form("EtDataCorrected%sPiPlus",cutName->Data()),"Corrected reconstructed E_{T} from identified #pi^{+}");
303     CreateEtaPtHisto2D(Form("EtDataCorrected%sPiMinus",cutName->Data()),"Corrected reconstructed E_{T} from identified #pi^{-}");
304     CreateEtaPtHisto2D(Form("EtDataCorrected%sKPlus",cutName->Data()),"Corrected reconstructed E_{T} from identified K^{+}");
305 //     CreateEtaPtHisto2D(Form("EtDataCorrected%sEMinus",cutName->Data()),"Corrected reconstructed E_{T} from identified e^{-}");
306 //     CreateEtaPtHisto2D(Form("EtDataCorrected%sEPlus",cutName->Data()),"Corrected reconstructed E_{T} from identified e^{+}");
307     CreateEtaPtHisto2D(Form("EtDataCorrected%sKMinus",cutName->Data()),"Corrected reconstructed E_{T} from identified K^{-}");
308     CreateEtaPtHisto2D(Form("EtDataCorrected%sProton",cutName->Data()),"Corrected reconstructed E_{T} from identified p");
309     CreateEtaPtHisto2D(Form("EtDataCorrected%sAntiProton",cutName->Data()),"Corrected reconstructed E_{T} from identified #bar{p}");
310     CreateEtaPtHisto2D(Form("EtDataCorrected%sUnidentified",cutName->Data()),"Corrected reconstructed E_{T} from unidentified particles using real mass");
311     CreateEtaPtHisto2D(Form("EtDataCorrected%sNoID",cutName->Data()),"Corrected reconstructed E_{T} from unidentified particles using real mass");
312
313
314     CreateEtaPtHisto2D(Form("EtNData%sPiPlus",cutName->Data()),"Number of reconstructed #pi^{+}");
315     CreateEtaPtHisto2D(Form("EtNData%sPiMinus",cutName->Data()),"Number of reconstructed #pi^{-}");
316     CreateEtaPtHisto2D(Form("EtNData%sKPlus",cutName->Data()),"Number of reconstructed K^{+}");
317     CreateEtaPtHisto2D(Form("EtNData%sKMinus",cutName->Data()),"Number of reconstructed K^{-}");
318     CreateEtaPtHisto2D(Form("EtNData%sProton",cutName->Data()),"Number of reconstructed p");
319     CreateEtaPtHisto2D(Form("EtNData%sAntiProton",cutName->Data()),"Number of reconstructed #bar{p}");
320     CreateEtaPtHisto2D(Form("EtNData%sUnidentified",cutName->Data()),"Number of Reconstructed unidentified particles");
321
322     CreateHisto2D(Form("dEdxDataAll%s",cutName->Data()),"dE/dx for all particles","momentum (GeV/c)","dE/dx",400,0.0,maxPtdEdx,200,mindEdx,maxdEdx);
323     CreateHisto2D(Form("dEdxDataPion%s",cutName->Data()),"dE/dx for #pi^{#pm}","momentum (GeV/c)","dE/dx",400,0.0,maxPtdEdx,200,mindEdx,maxdEdx);
324     CreateHisto2D(Form("dEdxDataKaon%s",cutName->Data()),"dE/dx for K^{#pm}","momentum (GeV/c)","dE/dx",400,0.0,maxPtdEdx,200,mindEdx,maxdEdx);
325     CreateHisto2D(Form("dEdxDataProton%s",cutName->Data()),"dE/dx for p(#bar{p})","momentum (GeV/c)","dE/dx",400,0.0,maxPtdEdx,200,mindEdx,maxdEdx);
326     CreateHisto2D(Form("dEdxDataElectron%s",cutName->Data()),"dE/dx for e^{#pm}","momentum (GeV/c)","dE/dx",400,0.0,maxPtdEdx,200,mindEdx,maxdEdx);
327     CreateHisto2D(Form("dEdxDataUnidentified%s",cutName->Data()),"dE/dx for unidentified particles","momentum (GeV/c)","dE/dx",400,0.0,maxPtdEdx,200,mindEdx,maxdEdx);
328   }
329 }