8f5f13df0cf45dc767303f137eb88c3992fc239a
[u/mrichter/AliRoot.git] / FMD / analysis / AliFMDAnalysisTaskSharing.cxx
1  
2 #include <TROOT.h>
3 #include <TSystem.h>
4 #include <TInterpreter.h>
5 #include <TChain.h>
6 #include <TFile.h>
7 #include <TList.h>
8 #include <iostream>
9 #include <TMath.h>
10 //#include "AliFMDDebug.h"
11 #include "AliFMDAnalysisTaskSharing.h"
12 #include "AliAnalysisManager.h"
13 #include "AliESDFMD.h"
14 //#include "AliFMDGeometry.h"
15 #include "AliMCEventHandler.h"
16 #include "AliMCEvent.h"
17 #include "AliStack.h"
18 #include "AliESDVertex.h"
19 #include "AliMultiplicity.h"
20 #include "AliFMDAnaParameters.h"
21 //#include "/home/canute/ALICE/AliRoot/PWG0/AliPWG0Helper.h"
22 //#include "AliFMDParameters.h"
23 #include "AliGenEventHeader.h"
24 #include "AliHeader.h"
25 #include "AliStack.h"
26 #include "AliMCParticle.h"
27
28 ClassImp(AliFMDAnalysisTaskSharing)
29
30 //_____________________________________________________________________
31 AliFMDAnalysisTaskSharing::AliFMDAnalysisTaskSharing()
32 : fDebug(0),
33   fESD(0x0),
34   foutputESDFMD(),
35   fEnergy(0),
36   fNstrips(0),
37   fSharedThis(kFALSE),
38   fSharedPrev(kFALSE),
39   fDiagList(0),
40   fStandalone(kTRUE),
41   fEsdVertex(0),
42   fStatus(kTRUE)
43 {
44   // Default constructor
45   DefineInput (0, AliESDEvent::Class());
46   DefineOutput(0, AliESDFMD::Class());
47   DefineOutput(1, AliESDVertex::Class());
48   DefineOutput(2, AliESDEvent::Class());
49   DefineOutput(3, TList::Class());
50 }
51 //_____________________________________________________________________
52 AliFMDAnalysisTaskSharing::AliFMDAnalysisTaskSharing(const char* name, Bool_t SE):
53     AliAnalysisTask(name, "AnalysisTaskFMD"),
54     fDebug(0),
55     fESD(0x0),
56     foutputESDFMD(),
57     fEnergy(0),
58     fNstrips(0),
59     fSharedThis(kFALSE),
60     fSharedPrev(kFALSE),
61     fDiagList(0),
62     fStandalone(kTRUE),
63     fEsdVertex(0),
64     fStatus(kTRUE)
65 {
66   fStandalone = SE;
67   if(fStandalone) {
68     DefineInput (0, AliESDEvent::Class());
69     DefineOutput(0, AliESDFMD::Class());
70     DefineOutput(1, AliESDVertex::Class());
71     DefineOutput(2, AliESDEvent::Class());
72     DefineOutput(3, TList::Class());
73   }
74 }
75 //_____________________________________________________________________
76 void AliFMDAnalysisTaskSharing::CreateOutputObjects()
77 {
78   if(!foutputESDFMD)
79     foutputESDFMD = new AliESDFMD();
80   
81   if(!fEsdVertex)
82     fEsdVertex    = new AliESDVertex();
83   //Diagnostics
84   if(!fDiagList)
85     fDiagList = new TList();
86   
87   fDiagList->SetName("Sharing diagnostics");
88   
89   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
90   TH2F* hBg   = pars->GetBackgroundCorrection(1, 'I', 0);
91   TH1F* hPrimary = new TH1F("hMultvsEtaNoCuts","hMultvsEtaNoCuts",
92                             hBg->GetNbinsX(),
93                             hBg->GetXaxis()->GetXmin(),
94                             hBg->GetXaxis()->GetXmax());
95   hPrimary->Sumw2();
96   fDiagList->Add(hPrimary);
97   TH1F* hPrimVertexBin = 0;
98   for(Int_t i = 0; i< pars->GetNvtxBins(); i++) {
99     
100     hPrimVertexBin = new TH1F(Form("primmult_NoCuts_vtxbin%d",i),
101                               Form("primmult_NoCuts_vtxbin%d",i),
102                               hBg->GetNbinsX(),
103                               hBg->GetXaxis()->GetXmin(),
104                               hBg->GetXaxis()->GetXmax());
105     hPrimVertexBin->Sumw2();
106     fDiagList->Add(hPrimVertexBin);
107     
108   }
109   
110   for(Int_t det = 1; det<=3; det++) {
111     Int_t nRings = (det==1 ? 1 : 2);
112     
113     for(Int_t iring = 0;iring<nRings; iring++) {
114       Char_t ringChar = (iring == 0 ? 'I' : 'O');
115       TH1F* hEdist        = new TH1F(Form("Edist_before_sharing_FMD%d%c", det, ringChar),
116                                      Form("Edist_before_sharing_FMD%d%c", det, ringChar),
117                                      1000,0,25);
118       TH1F* hEdist_after  = new TH1F(Form("Edist_after_sharing_FMD%d%c", det, ringChar),
119                                      Form("Edist_after_sharing_FMD%d%c", det, ringChar),
120                                      1000,0,25);
121       
122       
123       //TH1F* hNstripsHit    = new TH1F(Form("N_strips_hit_FMD%d%c",det,ringChar),
124       //                                     Form("N_strips_hit_FMD%d%c",det,ringChar),
125       //                                     25,0,25);
126       fDiagList->Add(hEdist);
127       fDiagList->Add(hEdist_after);
128       //fDiagList->Add(hNstripsHit);
129
130     }
131   }
132   TH1F*  nMCevents = new TH1F("nMCEventsNoCuts","nMCEventsNoCuts",pars->GetNvtxBins(),0,pars->GetNvtxBins());
133   
134   fDiagList->Add(nMCevents);
135   
136 }
137 //_____________________________________________________________________
138 void AliFMDAnalysisTaskSharing::ConnectInputData(Option_t */*option*/)
139 {
140   if(fStandalone)
141     fESD = (AliESDEvent*)GetInputData(0);
142 }
143 //_____________________________________________________________________
144 void AliFMDAnalysisTaskSharing::Exec(Option_t */*option*/)
145 {
146
147   AliESD* old = fESD->GetAliESDOld();
148   if (old) {
149     fESD->CopyFromOldESD();
150   }
151   
152   foutputESDFMD->Clear();
153   
154   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
155   Double_t vertex[3];
156   pars->GetVertex(fESD,vertex);
157   fEsdVertex->SetXYZ(vertex);
158   
159   // Process primaries here to get true MC distribution
160   if(pars->GetProcessPrimary())
161     ProcessPrimary();
162   
163   Bool_t isTriggered = pars->IsEventTriggered(fESD);
164   
165   if(!isTriggered) {
166     fStatus = kFALSE;
167     return;
168   }
169   else
170     fStatus = kTRUE;
171   
172   if(vertex[0] == 0 && vertex[1] == 0 && vertex[2] == 0) {
173     
174     fStatus = kFALSE;
175     return;
176   }
177   else
178     fStatus = kTRUE;
179   
180   
181   const AliMultiplicity* testmult = fESD->GetMultiplicity();
182   
183   Int_t nTrackLets = testmult->GetNumberOfTracklets();
184   if(nTrackLets < 1000) foutputESDFMD->SetUniqueID(kTRUE);
185   else foutputESDFMD->SetUniqueID(kFALSE);
186   
187   AliESDFMD* fmd = fESD->GetFMDData();
188   
189   if (!fmd) return;
190   Int_t nHits = 0;
191   for(UShort_t det=1;det<=3;det++) {
192     Int_t nRings = (det==1 ? 1 : 2);
193     for (UShort_t ir = 0; ir < nRings; ir++) {
194       Char_t   ring = (ir == 0 ? 'I' : 'O');
195       UShort_t nsec = (ir == 0 ? 20  : 40);
196       UShort_t nstr = (ir == 0 ? 512 : 256);
197       
198       TH1F* hEdist = (TH1F*)fDiagList->FindObject(Form("Edist_before_sharing_FMD%d%c",det,ring));
199       
200       for(UShort_t sec =0; sec < nsec;  sec++) {
201         fSharedThis      = kFALSE;
202         fSharedPrev      = kFALSE;
203         fEnergy = 0;
204         fNstrips = 0;
205         for(UShort_t strip = 0; strip < nstr; strip++) {
206           foutputESDFMD->SetMultiplicity(det,ring,sec,strip,0.);
207           Float_t mult = fmd->Multiplicity(det,ring,sec,strip);
208           
209           if(mult == AliESDFMD::kInvalidMult || mult == 0) continue;
210           
211           //Double_t eta  = EtaFromStrip(det,ring,sec,strip,vertex[2]);//fmd->Eta(det,ring,sec,strip);
212           //Double_t eta = fmd->Eta(det,ring,sec,strip);
213           Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex[2]);
214           //std::cout<<EtaFromStrip(det,ring,sec,strip,vertex[2]) <<"    "<<fmd->Eta(det,ring,sec,strip)<<std::endl;
215           
216           hEdist->Fill(mult);
217           if(fmd->IsAngleCorrected())
218             mult = mult/TMath::Cos(Eta2Theta(eta));
219           Float_t Eprev = 0;
220           Float_t Enext = 0;
221           if(strip != 0)
222             if(fmd->Multiplicity(det,ring,sec,strip-1) != AliESDFMD::kInvalidMult) {
223               Eprev = fmd->Multiplicity(det,ring,sec,strip-1);
224               if(fmd->IsAngleCorrected())
225                 Eprev = Eprev/TMath::Cos(Eta2Theta(fmd->Eta(det,ring,sec,strip-1)));
226             }
227           if(strip != nstr - 1)
228             if(fmd->Multiplicity(det,ring,sec,strip+1) != AliESDFMD::kInvalidMult) {
229               Enext = fmd->Multiplicity(det,ring,sec,strip+1);
230               if(fmd->IsAngleCorrected())
231                 Enext = Enext/TMath::Cos(Eta2Theta(fmd->Eta(det,ring,sec,strip+1)));
232             }
233           
234           Float_t merged_energy = GetMultiplicityOfStrip(mult,eta,Eprev,Enext,det,ring,sec,strip);
235
236           if(merged_energy > 0 )
237             nHits++;
238           foutputESDFMD->SetMultiplicity(det,ring,sec,strip,merged_energy);
239           foutputESDFMD->SetEta(det,ring,sec,strip,eta);
240           
241         }
242       }
243     }
244   }
245   
246    
247   if(fStandalone) {
248     PostData(0, foutputESDFMD); 
249     PostData(1, fEsdVertex); 
250     PostData(2, fESD); 
251     PostData(3, fDiagList); 
252   }
253 }
254 //_____________________________________________________________________
255 Float_t AliFMDAnalysisTaskSharing::GetMultiplicityOfStrip(Float_t mult,
256                                                           Float_t eta,
257                                                           Float_t Eprev,
258                                                           Float_t Enext,
259                                                           UShort_t   det,
260                                                           Char_t  ring,
261                                                           UShort_t /*sec*/,
262                                                           UShort_t /*strip*/) {
263   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
264  
265   Float_t merged_energy = 0;
266   //Float_t nParticles = 0;
267   Float_t cutLow  = 0.15;
268   if(ring == 'I')
269     cutLow = 0.1;
270   
271   //cutLow = 0;
272   //AliFMDParameters* recopars = AliFMDParameters::Instance();
273   //cutLow = (5*recopars->GetPedestalWidth(det,ring,sec,strip))/(recopars->GetPulseGain(det,ring,sec,strip)*recopars->GetDACPerMIP());
274   
275   
276   
277   Float_t cutHigh = pars->GetMPV(det,ring,eta) - 3*pars->GetSigma(det,ring,eta);
278   // Float_t cutPart = pars->GetMPV(det,ring,eta) - 5*pars->GetSigma(det,ring,eta);
279   Float_t Etotal  = mult;
280   
281
282     //std::cout<<det<<ring<<"   "<<sec<<"    "<<strip<<"   "<<cutLow<<std::endl;
283   if(fSharedThis) {
284     fSharedThis      = kFALSE;
285     fSharedPrev      = kTRUE;
286     return 0.;
287   }
288   
289   /*  if(mult < 0.33*pars->GetMPV(det,ring,eta)) {
290     fSharedThis      = kFALSE;
291     fSharedPrev      = kFALSE;
292     return 0;
293     }*/
294   if(mult<Enext && Enext>cutHigh && foutputESDFMD->GetUniqueID() == kTRUE)
295     {
296       fSharedThis      = kFALSE;
297       fSharedPrev      = kFALSE;
298       return 0;
299     }
300   if(mult > 15)
301     {
302       //   std::cout<<"rejecting hit in FMD "<<det<<" "<<ring<<std::endl;
303       fSharedThis      = kFALSE;
304       fSharedPrev      = kFALSE;
305       return 0;
306     }
307   
308   if(Eprev > cutLow && Eprev < cutHigh && !fSharedPrev ) {
309     Etotal += Eprev;
310   }
311   
312   if(Enext > cutLow && Enext < cutHigh ) {
313     Etotal += Enext;
314     fSharedThis      = kTRUE;
315   }
316   TH1F* hEdist = (TH1F*)fDiagList->FindObject(Form("Edist_after_sharing_FMD%d%c",det,ring));
317   hEdist->Fill(Etotal);
318   
319   Etotal = Etotal*TMath::Cos(Eta2Theta(eta));
320   if(Etotal > 0) {
321     
322     merged_energy = Etotal;
323     fSharedPrev      = kTRUE;
324     // if(det == 1 && ring =='I')
325     // std::cout<<Form("Merged signals %f %f %f into %f , %f in strip %d, sec %d, ring %c, det %d",Eprev, mult, Enext, Etotal/TMath::Cos(Eta2Theta(eta)),Etotal,strip,sec,ring,det )<<std::endl;
326   }
327     else{// if(Etotal > 0) {
328       //if(det == 3 && ring =='I')
329       //        std::cout<<Form("NO HIT  for  %f %f %f into %f , %f in strip %d, sec %d, ring %c, det %d, cuts %f , %f",Eprev, mult, Enext, Etotal/TMath::Cos(Eta2Theta(eta)),Etotal,strip,sec,ring,det,cutPart,cutHigh )<<std::endl;
330     fSharedThis      = kFALSE;
331     fSharedPrev      = kFALSE;
332   }
333   // merged_energy = mult;
334   
335   return merged_energy; 
336   //}  
337 }
338
339 //_____________________________________________________________________
340 Float_t AliFMDAnalysisTaskSharing::Eta2Theta(Float_t eta) {
341
342   Float_t theta = 2*TMath::ATan(TMath::Exp(-1*eta));
343   
344   if(eta < 0)
345     theta = theta-TMath::Pi();
346   
347   //  std::cout<<"From eta2Theta: "<<theta<<"   "<<eta<<std::endl;
348   return theta;
349   
350
351
352 }
353
354
355
356 //_____________________________________________________________________
357 void AliFMDAnalysisTaskSharing::ProcessPrimary() {
358   
359   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
360   AliMCEvent* mcEvent = eventHandler->MCEvent();
361   if(!mcEvent)
362     return;
363   
364   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
365   
366   AliMCParticle* particle = 0;
367   
368   AliStack* stack = mcEvent->Stack();
369   
370   TH1F* hPrimary = (TH1F*)fDiagList->FindObject("hMultvsEtaNoCuts");
371   AliHeader* header            = mcEvent->Header();
372   AliGenEventHeader* genHeader = header->GenEventHeader();
373   
374   
375   
376   TArrayF vertex;
377   genHeader->PrimaryVertex(vertex);
378   if(TMath::Abs(vertex.At(2)) > pars->GetVtxCutZ())
379     return;
380   
381   Double_t delta           = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
382   Double_t vertexBinDouble = (vertex.At(2) + pars->GetVtxCutZ()) / delta;
383   Int_t    vertexBin       = (Int_t)vertexBinDouble;
384   
385   Bool_t firstTrack = kTRUE;
386   
387   Int_t nTracks = stack->GetNprimary();
388   TH1F* nMCevents = (TH1F*)fDiagList->FindObject("nMCEventsNoCuts");
389   for(Int_t i = 0 ;i<nTracks;i++) {
390     particle = mcEvent->GetTrack(i);
391     if(!particle)
392       continue;
393     
394     if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
395       hPrimary->Fill(particle->Eta());
396       
397
398       TH1F* hPrimVtxBin = (TH1F*)fDiagList->FindObject(Form("primmult_NoCuts_vtxbin%d",vertexBin));
399       hPrimVtxBin->Fill(particle->Eta());
400       
401       if(firstTrack) {
402         nMCevents->Fill(vertexBin);
403         firstTrack = kFALSE;
404       }
405     
406     }
407   }
408
409 }
410
411 //_____________________________________________________________________
412 //
413 // EOF
414 //