43a4b3253b0a46b459c834aafab0980bf0accf87
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / photons / AliAnalysisTaskPMDSim.cxx
1 #include "TChain.h"
2 #include "TTree.h"
3 #include "TH1F.h"
4 #include "TH2F.h"
5 #include "TCanvas.h"
6 #include "AliAnalysisTask.h"
7 #include "AliAnalysisManager.h"
8 #include "AliESDEvent.h"
9 #include "AliESDInputHandler.h"
10 #include "AliESDPmdTrack.h"
11 #include "AliESDVertex.h"
12 #include "AliAnalysisTaskPMDSim.h"
13 #include <TParticle.h>
14 #include <AliLog.h>
15 #include <AliStack.h>
16 #include <AliHeader.h>
17 #include <AliGenEventHeader.h>
18 #include <AliHeader.h>
19 #include <AliGenPythiaEventHeader.h>
20 #include <AliMCEventHandler.h>
21 #include <AliMCEvent.h>
22 #include "AliTriggerAnalysis.h"
23 #include "AliPhysicsSelection.h"
24
25 // AnalysisTask For PMD
26 // Authors: Sudipan De, Subhash Singha
27
28 ClassImp(AliAnalysisTaskPMDSim)
29
30 //________________________________________________________________________
31 AliAnalysisTaskPMDSim::AliAnalysisTaskPMDSim(const char *name) 
32   : AliAnalysisTaskSE(name), 
33   fESD(0), 
34   fOutputList(0), 
35   fHistTotEvent(0),
36   fHistTotEventAfterPhySel(0),
37   fHistTotEventAfterVtx(0),
38   fVtxZ(0),
39   fHistXYPre(0),
40   fHistEtaPhM(0),
41   fHistEtaPhM1(0),
42   fHistEtaT(0),
43   fMultMeasured(0),
44   fMultMeasured1(0),
45   fMultTrue(0),
46   fMultCorr(0),
47   fMultCorr1(0)
48 {
49   for(Int_t i=0; i<10; i++){
50     fHistMultMeasEtaBinA[i] = 0;
51     fHistMultMeasEtaBinA1[i] = 0;
52     fHistMultTrueEtaBinA[i] = 0;
53     fHistMultCorrEtaBinA[i] = 0;
54     fHistMultCorrEtaBinA1[i] = 0;
55   }
56   
57   // Constructor
58   
59   // Define input and output slots here
60   // Input slot #0 works with a TChain
61   DefineInput(0, TChain::Class());
62   // Output slot #0 id reserved by the base class for AOD
63   // Output slot #1 writes into a TH1 container
64   DefineOutput(1, TList::Class());
65 }
66
67 //________________________________________________________________________
68 void AliAnalysisTaskPMDSim::UserCreateOutputObjects()
69 {
70   // Create histograms
71   // Called once
72   
73   fOutputList = new TList();
74   
75   Int_t kNbinsMult = 100;  Float_t XminMult = -0.5; 
76   Float_t XmaxMult = 99.5;//mult
77   Int_t kNbinsMultA = 50;  Float_t XminMultA = -0.5; 
78   Float_t XmaxMultA = 49.5;//mult
79   Int_t kNbinEtaPMDCov = 10; Float_t XminEtaPMDCov = 2.1; Float_t XmaxEtaPMDCov = 4.1;//etaPMD
80   Int_t kNBinsEvent = 10; Float_t XminEvent = 0; Float_t XmaxEvent = 10;
81   Int_t kNbinsXY = 200; Float_t XminXY = -100.0; Float_t XmaxXY  = 100.0;
82   fHistTotEvent = new TH1F("TotEvent","TotEvent",
83                            kNBinsEvent,XminEvent,XmaxEvent); 
84   fOutputList->Add(fHistTotEvent);
85   fHistTotEventAfterPhySel = new TH1F("TotEventAfterPhySel","TotEventAfterPhySel",
86                                       kNBinsEvent,XminEvent,XmaxEvent); 
87   fOutputList->Add(fHistTotEventAfterPhySel);
88   fHistTotEventAfterVtx = new TH1F("TotEventAfterVtx","TotEventAfterVtx",
89                                   kNBinsEvent,XminEvent,XmaxEvent); 
90   fOutputList->Add(fHistTotEventAfterVtx);
91   fVtxZ     = new TH1F("VtxZ","VtxZ",200,-50.0,50.0);
92   fOutputList->Add(fVtxZ);
93   fHistXYPre = new TH2F("XYPre","XYPre",kNbinsXY,XminXY,XmaxXY,
94                         kNbinsXY,XminXY,XmaxXY);
95   fOutputList->Add(fHistXYPre);
96   fHistEtaPhM = new TH1F("fHistEtaPhM","fHistEtaPhM",kNbinEtaPMDCov,XminEtaPMDCov,XmaxEtaPMDCov);
97   fOutputList->Add(fHistEtaPhM);
98   fHistEtaPhM1 = new TH1F("fHistEtaPhM1","fHistEtaPhM1",kNbinEtaPMDCov,XminEtaPMDCov,XmaxEtaPMDCov);
99   fOutputList->Add(fHistEtaPhM1);
100   fHistEtaT = new TH1F("fHistEtaT","fHistEtaT",kNbinEtaPMDCov,XminEtaPMDCov,XmaxEtaPMDCov);
101   fOutputList->Add(fHistEtaT);
102   fMultMeasured = new TH1F("MultM","MultM",kNbinsMult,XminMult,XmaxMult);
103   fOutputList->Add(fMultMeasured);
104   fMultMeasured1 = new TH1F("MultM1","MultM1",kNbinsMult,XminMult,XmaxMult);
105   fOutputList->Add(fMultMeasured1);
106   fMultTrue = new TH1F("MultT","MultT",kNbinsMult,XminMult,XmaxMult);
107   fOutputList->Add(fMultTrue);
108   fMultCorr = new TH2F("MultCorr","MultCorr",kNbinsMult,XminMult,XmaxMult,
109                        kNbinsMult,XminMult,XmaxMult);
110   fOutputList->Add(fMultCorr);
111   fMultCorr1 = new TH2F("MultCorr1","MultCorr1",kNbinsMult,XminMult,XmaxMult,
112                        kNbinsMult,XminMult,XmaxMult);
113   fOutputList->Add(fMultCorr1);
114    
115   
116   Char_t nameT[256], nameM[256], nameCorr[256],nameM1[256], nameCorr1[256];
117   for(Int_t i=0; i<10; i++){
118     sprintf(nameM,"MultM_EtaBin%d",i+1);
119     sprintf(nameM1,"MultM1_EtaBin%d",i+1);
120     sprintf(nameT,"MultT_EtaBin%d",i+1);
121     sprintf(nameCorr,"MultCorr_EtaBin%d",i+1);
122     sprintf(nameCorr1,"MultCorr1_EtaBin%d",i+1);
123     fHistMultMeasEtaBinA[i] = 
124       new TH1F(nameM,nameM,kNbinsMultA,XminMultA,XmaxMultA);
125     fHistMultMeasEtaBinA1[i] = 
126       new TH1F(nameM1,nameM1,kNbinsMultA,XminMultA,XmaxMultA);
127     fHistMultTrueEtaBinA[i] = 
128       new TH1F(nameT,nameT,kNbinsMultA,XminMultA,XmaxMultA);
129     fHistMultCorrEtaBinA[i] = 
130       new TH2F(nameCorr,nameCorr,kNbinsMultA,XminMultA,XmaxMultA,
131                kNbinsMultA,XminMultA,XmaxMultA);
132     fHistMultCorrEtaBinA1[i] = 
133       new TH2F(nameCorr1,nameCorr1,kNbinsMultA,XminMultA,XmaxMultA,
134                kNbinsMultA,XminMultA,XmaxMultA);
135     fOutputList->Add(fHistMultMeasEtaBinA[i]);
136     fOutputList->Add(fHistMultMeasEtaBinA1[i]);
137     fOutputList->Add(fHistMultTrueEtaBinA[i]);
138     fOutputList->Add(fHistMultCorrEtaBinA[i]);
139     fOutputList->Add(fHistMultCorrEtaBinA1[i]);
140   }//i loop
141   //fOutputList->Add(fHistPt);
142 }
143
144 //________________________________________________________________________
145 void AliAnalysisTaskPMDSim::UserExec(Option_t *) 
146 {
147   // Main loop
148   // Called for each event
149   Float_t MipCut1 = 432;//MPV=432=> 6*MPV=432
150   Float_t MipCut2 = 648;//MPV=432=> 6*MPV=648 
151   Float_t etacls=0., theta=0., rdist=0.;
152   Int_t PhotonCls = 0;
153   Int_t PhotonCls1 = 0;
154   Int_t PhotonClsAEtaBin[10] = {0};//#of photon measured with diff Eta bin
155   Int_t PhotonClsAEtaBin1[10] = {0};//#of photon measured with diff Eta bin
156   Int_t PhotonTrueAEtaBin[10] = {0};//#of photon Incident with diff Eta bin
157   // Post output data.
158   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
159   if (!fESD) {
160     printf("ERROR: fESD not available\n");
161     return;
162   }
163   fHistTotEvent->Fill(5);
164   //event selection
165   Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
166   if (! isSelected) return;
167   fHistTotEventAfterPhySel->Fill(5);
168   //Vertex selection
169   const AliESDVertex *vertex = fESD->GetPrimaryVertex();
170   Float_t Vz = vertex->GetZv();    
171   // Float_t Vx = vertex->GetXv();    
172   //Float_t Vy = vertex->GetYv();
173   Bool_t zVerStatus = vertex->GetStatus();
174   if(zVerStatus){
175     fVtxZ->Fill(Vz);
176     if(TMath::Abs(Vz)<10){   
177       fHistTotEventAfterVtx->Fill(5);   
178       Int_t ptracks = fESD->GetNumberOfPmdTracks();
179       for(Int_t kk=0;kk<ptracks;kk++){
180         AliESDPmdTrack *pmdtr = fESD->GetPmdTrack(kk);
181         Int_t   det   = pmdtr->GetDetector();
182         Float_t clsX  = pmdtr->GetClusterX();
183         Float_t clsY  = pmdtr->GetClusterY();
184         Float_t clsZ  = pmdtr->GetClusterZ();
185         clsZ -= Vz;
186         Float_t ncell = pmdtr->GetClusterCells();
187         Float_t adc   = pmdtr->GetClusterADC();
188         //Float_t pid   = pmdtr->GetClusterPID();
189         //Float_t isotag = pmdtr->GetClusterSigmaX();
190         //Int_t trno = pmdtr->GetClusterTrackNo();
191         //calculation of #eta
192         rdist = TMath::Sqrt(clsX*clsX + clsY*clsY);
193         if(clsZ!=0) theta = TMath::ATan2(rdist,clsZ);
194         etacls   = -TMath::Log(TMath::Tan(0.5*theta));
195         
196         if(det==0 && adc>0)fHistXYPre->Fill(clsX,clsY);
197         if(det==0 && adc>MipCut1 && ncell>2){
198           fHistEtaPhM->Fill(etacls);
199           if(etacls > 2.3 && etacls <= 3.9)  PhotonCls++;
200           if(etacls > 2.1 && etacls <= 2.3)  PhotonClsAEtaBin[0]++;
201           if(etacls > 2.3 && etacls <= 2.5)  PhotonClsAEtaBin[1]++;
202           if(etacls > 2.5 && etacls <= 2.7)  PhotonClsAEtaBin[2]++;
203           if(etacls > 2.7 && etacls <= 2.9)  PhotonClsAEtaBin[3]++;
204           if(etacls > 2.9 && etacls <= 3.1)  PhotonClsAEtaBin[4]++;
205           if(etacls > 3.1 && etacls <= 3.3)  PhotonClsAEtaBin[5]++;
206           if(etacls > 3.3 && etacls <= 3.5)  PhotonClsAEtaBin[6]++;
207           if(etacls > 3.5 && etacls <= 3.7)  PhotonClsAEtaBin[7]++;
208           if(etacls > 3.7 && etacls <= 3.9)  PhotonClsAEtaBin[8]++;
209           if(etacls > 3.9 && etacls <= 4.1)  PhotonClsAEtaBin[9]++;
210         }//if MipCut1
211         if( det==0 && adc>MipCut2 && ncell>2){
212           fHistEtaPhM1->Fill(etacls);
213           if(etacls > 2.3 && etacls <= 3.9)  PhotonCls1++;
214           if(etacls > 2.1 && etacls <= 2.3)  PhotonClsAEtaBin1[0]++;
215           if(etacls > 2.3 && etacls <= 2.5)  PhotonClsAEtaBin1[1]++;
216           if(etacls > 2.5 && etacls <= 2.7)  PhotonClsAEtaBin1[2]++;
217           if(etacls > 2.7 && etacls <= 2.9)  PhotonClsAEtaBin1[3]++;
218           if(etacls > 2.9 && etacls <= 3.1)  PhotonClsAEtaBin1[4]++;
219           if(etacls > 3.1 && etacls <= 3.3)  PhotonClsAEtaBin1[5]++;
220           if(etacls > 3.3 && etacls <= 3.5)  PhotonClsAEtaBin1[6]++;
221           if(etacls > 3.5 && etacls <= 3.7)  PhotonClsAEtaBin1[7]++;
222           if(etacls > 3.7 && etacls <= 3.9)  PhotonClsAEtaBin1[8]++;
223           if(etacls > 3.9 && etacls <= 4.1)  PhotonClsAEtaBin1[9]++;
224         }//if MipCut2
225       } //track loop 
226       //reading MC info.
227       AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*>
228         (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
229       if (!eventHandler) {
230         Printf("ERROR: Could not retrieve MC event handler");
231         return;
232       }
233       AliMCEvent* mcEvent = eventHandler->MCEvent();
234       if (!mcEvent) {
235         Printf("ERROR: Could not retrieve MC event");
236         return;
237       }
238       AliStack* stack = mcEvent->Stack();
239       if (!stack)
240         {
241           AliDebug(AliLog::kError, "Stack not available");
242           return;
243         }
244       if(stack){
245         Int_t nPrim  = stack->GetNprimary();
246         Int_t PhotonTrue = 0;
247         for (Int_t iMc = 0; iMc < nPrim; ++iMc)
248           {
249             TParticle *MPart = stack->Particle(iMc);
250             Int_t mpart  = MPart->GetPdgCode();
251             Float_t eta   = MPart->Eta();
252             if(mpart == 22){  
253               fHistEtaT->Fill(eta);
254               if(eta > 2.3 && eta <= 3.9)PhotonTrue++;
255               if(eta > 2.1 && eta <= 2.3)PhotonTrueAEtaBin[0]++;
256               if(eta > 2.3 && eta <= 2.5)PhotonTrueAEtaBin[1]++;
257               if(eta > 2.5 && eta <= 2.7)PhotonTrueAEtaBin[2]++;
258               if(eta > 2.7 && eta <= 2.9)PhotonTrueAEtaBin[3]++;
259               if(eta > 2.9 && eta <= 3.1)PhotonTrueAEtaBin[4]++;
260               if(eta > 3.1 && eta <= 3.3)PhotonTrueAEtaBin[5]++;
261               if(eta > 3.3 && eta <= 3.5)PhotonTrueAEtaBin[6]++;
262               if(eta > 3.5 && eta <= 3.7)PhotonTrueAEtaBin[7]++;
263               if(eta > 3.7 && eta <= 3.9)PhotonTrueAEtaBin[8]++;
264               if(eta > 3.9 && eta <= 4.1)PhotonTrueAEtaBin[9]++;
265             }//mpart
266           }//track loop
267         for(Int_t i=0; i<10; i++){
268           fHistMultMeasEtaBinA[i]->Fill(PhotonClsAEtaBin[i]);
269           fHistMultMeasEtaBinA1[i]->Fill(PhotonClsAEtaBin1[i]);
270           fHistMultTrueEtaBinA[i]->Fill(PhotonTrueAEtaBin[i]);
271           fHistMultCorrEtaBinA[i]->Fill(PhotonTrueAEtaBin[i],PhotonClsAEtaBin[i]);
272           fHistMultCorrEtaBinA1[i]->Fill(PhotonTrueAEtaBin[i],PhotonClsAEtaBin1[i]);
273         }//i loop
274         fMultMeasured->Fill(PhotonCls);
275         fMultMeasured1->Fill(PhotonCls1);
276         fMultTrue->Fill(PhotonTrue);
277         fMultCorr->Fill(PhotonTrue,PhotonCls);
278         fMultCorr1->Fill(PhotonTrue,PhotonCls1);
279       }//if stack
280     }//if Vz!= 0
281   }//Vz loop
282   PostData(1, fOutputList);
283 }//UserExec()     
284
285 //_______________________________________________________________________
286 void AliAnalysisTaskPMDSim::Terminate(Option_t *) 
287 {
288   // Draw result to the screen
289   // Called once at the end of the query
290   fOutputList = dynamic_cast<TList*> (GetOutputData(1));
291   if (!fOutputList) {
292     printf("ERROR: Output list not available\n");
293     return;
294   }
295 }//Terminate