Next round of upgrades and cleanups of code. The code is now more streamlined and...
[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    
279   // Float_t cutPart = pars->GetMPV(det,ring,eta) - 5*pars->GetSigma(det,ring,eta);
280   Float_t Etotal  = mult;
281   
282
283     //std::cout<<det<<ring<<"   "<<sec<<"    "<<strip<<"   "<<cutLow<<std::endl;
284   if(fSharedThis) {
285     fSharedThis      = kFALSE;
286     fSharedPrev      = kTRUE;
287     return 0.;
288   }
289   
290   /*  if(mult < 0.33*pars->GetMPV(det,ring,eta)) {
291     fSharedThis      = kFALSE;
292     fSharedPrev      = kFALSE;
293     return 0;
294     }*/
295   if(mult<Enext && Enext>cutHigh && foutputESDFMD->GetUniqueID() == kTRUE)
296     {
297       fSharedThis      = kFALSE;
298       fSharedPrev      = kFALSE;
299       return 0;
300     }
301   if(mult > 15)
302     {
303       //   std::cout<<"rejecting hit in FMD "<<det<<" "<<ring<<std::endl;
304       fSharedThis      = kFALSE;
305       fSharedPrev      = kFALSE;
306       return 0;
307     }
308   
309   if(Eprev > cutLow && Eprev < cutHigh && !fSharedPrev ) {
310     Etotal += Eprev;
311   }
312   
313   if(Enext > cutLow && Enext < cutHigh ) {
314     Etotal += Enext;
315     fSharedThis      = kTRUE;
316   }
317   TH1F* hEdist = (TH1F*)fDiagList->FindObject(Form("Edist_after_sharing_FMD%d%c",det,ring));
318   hEdist->Fill(Etotal);
319   
320   Etotal = Etotal*TMath::Cos(Eta2Theta(eta));
321   if(Etotal > 0) {
322     
323     merged_energy = Etotal;
324     fSharedPrev      = kTRUE;
325     // if(det == 1 && ring =='I')
326     // 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;
327   }
328     else{// if(Etotal > 0) {
329       //if(det == 3 && ring =='I')
330       //        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;
331     fSharedThis      = kFALSE;
332     fSharedPrev      = kFALSE;
333   }
334   // merged_energy = mult;
335   
336   return merged_energy; 
337   //}  
338 }
339
340 //_____________________________________________________________________
341 Float_t AliFMDAnalysisTaskSharing::Eta2Theta(Float_t eta) {
342
343   Float_t theta = 2*TMath::ATan(TMath::Exp(-1*eta));
344   
345   if(eta < 0)
346     theta = theta-TMath::Pi();
347   
348   //  std::cout<<"From eta2Theta: "<<theta<<"   "<<eta<<std::endl;
349   return theta;
350   
351
352
353 }
354
355
356
357 //_____________________________________________________________________
358 void AliFMDAnalysisTaskSharing::ProcessPrimary() {
359   
360   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
361   AliMCEvent* mcEvent = eventHandler->MCEvent();
362   if(!mcEvent)
363     return;
364   
365   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
366   
367   AliMCParticle* particle = 0;
368   
369   AliStack* stack = mcEvent->Stack();
370   
371   TH1F* hPrimary = (TH1F*)fDiagList->FindObject("hMultvsEtaNoCuts");
372   AliHeader* header            = mcEvent->Header();
373   AliGenEventHeader* genHeader = header->GenEventHeader();
374   
375   
376   
377   TArrayF vertex;
378   genHeader->PrimaryVertex(vertex);
379   if(TMath::Abs(vertex.At(2)) > pars->GetVtxCutZ())
380     return;
381   
382   Double_t delta           = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
383   Double_t vertexBinDouble = (vertex.At(2) + pars->GetVtxCutZ()) / delta;
384   Int_t    vertexBin       = (Int_t)vertexBinDouble;
385   
386   Bool_t firstTrack = kTRUE;
387   
388   Int_t nTracks = stack->GetNprimary();
389   TH1F* nMCevents = (TH1F*)fDiagList->FindObject("nMCEventsNoCuts");
390   for(Int_t i = 0 ;i<nTracks;i++) {
391     particle = mcEvent->GetTrack(i);
392     if(!particle)
393       continue;
394     
395     if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
396       hPrimary->Fill(particle->Eta());
397       
398
399       TH1F* hPrimVtxBin = (TH1F*)fDiagList->FindObject(Form("primmult_NoCuts_vtxbin%d",vertexBin));
400       hPrimVtxBin->Fill(particle->Eta());
401       
402       if(firstTrack) {
403         nMCevents->Fill(vertexBin);
404         firstTrack = kFALSE;
405       }
406     
407     }
408   }
409
410 }
411
412 //_____________________________________________________________________
413 //
414 // EOF
415 //