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