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