wagon to the PWGPP train
[u/mrichter/AliRoot.git] / T0 / AliT0CalibAnalysisTask.cxx
1 #include "TChain.h"
2 #include "TTree.h"
3 #include "TList.h"
4 #include "TH1F.h"
5 #include "TH2F.h"
6 #include "TCanvas.h"
7 #include "TMath.h"
8 #include "TObjString.h"
9
10 #include "AliAnalysisTask.h"
11 #include "AliPhysicsSelection.h"
12 #include "AliAnalysisManager.h"
13
14 #include "AliESDEvent.h"
15 #include "AliESDfriend.h"
16 #include "AliESDInputHandler.h"
17 #include "AliESDTZEROfriend.h"
18 #include "AliT0CalibAnalysisTask.h"
19 #include "AliESDVZERO.h"
20 #include "AliMultiplicity.h"
21 #include "AliTriggerAnalysis.h"
22 #include "AliESDpid.h"
23 #include "AliESDtrackCuts.h"
24 #include "AliCentrality.h"
25 #include "AliVEvent.h"
26
27 ClassImp(AliT0CalibAnalysisTask)
28
29 //________________________________________________________________________
30 AliT0CalibAnalysisTask::AliT0CalibAnalysisTask(const char *name) 
31   : AliAnalysisTaskSE(name), fESD(0), fOutputList(0), 
32   fT0OutTree(0x0), fEvent(-99999),   fOrbit(-99999), fBC(-99999),
33   fTrackletSPD(-99999), fClustersSPD(-99999), fNcont(-99999), fNcontTPC(-99999),
34   fVertex(-99999), fVertexcalc(-99999), fVertexSPD(-99999), fVertexTPC(-99999),
35   fMeanAC(-99999), fMeanA(-99999),fMeanC(-99999), fMeanACcalc(-99999), fMeanAcalc(-99999), fMeanCcalc(-99999),
36   fMultV0A(-99999), fMultV0C(-99999),fTimeV0A(-99999),fTimeV0C(-99999), fSumampA(-99999), fSumampC(-99999),
37   ftimestamp(0), fSep2(0),
38   fZDCcut(kFALSE), fT0Trigger(0), fpileup(kFALSE), fTrigger(0x0),
39   fT0_amplitude(0x0), fT0_time(0x0),
40   fcentralityV0M(0), fcentralityZDC(0), fcentralityTRK(0), fcentralityCLA(0),
41   fESDtracks(0),
42   fMultiplicity(-99999),
43   fTriggerinput(0x0), fZDCbg(kFALSE),
44   fTOFtracks(0), fT0tofTrack(0),
45   fESDpid(new AliESDpid())
46 {
47   // Constructor
48
49   // Define input and output slots here
50   // Input slot #0 works with a TChain
51   DefineInput(0, TChain::Class());
52   // Output slot #0 id reserved by the base class for AOD
53   // Output slot #1 writes into a TH1 container
54   DefineOutput(1, TList::Class());
55 }
56
57 //________________________________________________________________________
58 AliT0CalibAnalysisTask::~AliT0CalibAnalysisTask() 
59 {
60 // Destructor
61    if (fT0OutTree) delete fT0OutTree ;
62 }
63
64 //_____________________________________________________________________________
65 Bool_t AliT0CalibAnalysisTask::UserNotify()
66 {
67 //
68 // Calls the mother class Notify()
69 //
70
71 //  AliDebug(AliLog::kDebug,"<-");
72   // AliDebug(AliLog::kDebug,"->");
73
74   return AliAnalysisTaskSE::UserNotify();
75 }
76
77 //________________________________________________________________________
78 void AliT0CalibAnalysisTask::UserCreateOutputObjects()
79 {
80   // Create histograms
81   // Called once
82    fOutputList = new TList();
83   fOutputList->SetOwner(); // Will delete the histos on cleanup
84
85   fT0OutTree=new TTree("t0tree","None here");
86   fT0OutTree->Branch("fNevent", &fEvent);
87   fT0OutTree->Branch("fOrbit", &fOrbit);
88   fT0OutTree->Branch("fBC", &fBC);
89   fT0OutTree->Branch("fNcont", &fNcont, "fNcont/I");
90   fT0OutTree->Branch("fNcontTPC", &fNcontTPC);
91   fT0OutTree->Branch("fVertexTPC", &fVertexTPC);
92   fT0OutTree->Branch("vertexT0", &fVertex, "fVertex/F");
93   fT0OutTree->Branch("vertexSPD", &fVertexSPD, "vertexSPD/F");
94   fT0OutTree->Branch("meanAC", &fMeanAC,"meanAC/F");
95   fT0OutTree->Branch("meanA", &fMeanA,"meanA/F");
96   fT0OutTree->Branch("meanC", &fMeanC,"meanC/F");
97   fT0OutTree->Branch("trackletSPD", &fTrackletSPD,"trackletSPD/I");
98   fT0OutTree->Branch("clustersSPD", &fClustersSPD,"clustersSPD/I");
99   fT0OutTree->Branch("TOFtracks", &fTOFtracks);
100   fT0OutTree->Branch("ESDtracks", &fESDtracks);
101   fT0OutTree->Branch("multV0A", &fMultV0A,"multV0A/F");
102   fT0OutTree->Branch("multV0C", &fMultV0C,"multV0C/F");
103   fT0OutTree->Branch("timeV0A", &fTimeV0A,"timeV0A/F");
104   fT0OutTree->Branch("timeV0C", &fTimeV0C,"timeV0C/F");
105   fT0OutTree->Branch("sumampA", &fSumampA);
106   fT0OutTree->Branch("sumampC", &fSumampC);
107   fT0OutTree->Branch("pileup", &fpileup);
108   fT0OutTree->Branch("trigger", &fTrigger);
109   fT0OutTree->Branch("triggerinput", &fTriggerinput);
110   fT0OutTree->Branch("T0trigger", &fT0Trigger);
111   fT0OutTree->Branch("t0tofTrack",  &fT0tofTrack);
112   fT0OutTree->Branch("ZDCcut", &fZDCcut);
113   fT0OutTree->Branch("ftimestamp", &ftimestamp);
114   fT0OutTree->Branch("centralityV0M", &fcentralityV0M);
115   fT0OutTree->Branch("centralityZDC", &fcentralityZDC);
116   fT0OutTree->Branch("centralityTRK", &fcentralityTRK);
117   fT0OutTree->Branch("centralityCLA", &fcentralityCLA);
118   for ( Int_t i=0; i<24; i++) {
119     fT0OutTree->Branch(Form("amp%i", i+1), &famp[i]);
120     fT0OutTree->Branch(Form("time%i", i+1), &ftime[i]);
121     fT0OutTree->Branch(Form("ftimefriend%i", i+1), &ftimefriend[i]);
122     fT0OutTree->Branch(Form("fampfriend%i", i+1), &fampfriend[i]);
123     fT0OutTree->Branch(Form("fampledfriend%i", i+1), &fampledfriend[i]);
124   } 
125   for ( Int_t i=0; i<5; i++) {
126     fT0OutTree->Branch(Form("fOrA%i", i+1), &fOrA[i]);
127     fT0OutTree->Branch(Form("fOrC%i", i+1), &fOrC[i]);
128     fT0OutTree->Branch(Form("fTVDC%i", i+1), &fTVDC[i]);  
129     for ( Int_t ii=0; ii<24; ii++) 
130       fT0OutTree->Branch(Form("fRawTime%i_%i", ii+1, i+1), &fRawTime[ii][i]);
131   }
132   TString pilename[3] = {"T0pileup", "T0background", "T0satellite"};
133   for ( Int_t i=0; i<3; i++) 
134     fT0OutTree->Branch(pilename[i].Data(), &fT0pileup[i]);
135   fT0OutTree->Branch("meanACbest", &fMeanACcalc,"meanACcalc/F");
136   fT0OutTree->Branch("meanAbest", &fMeanAcalc,"meanAcalc/F");
137   fT0OutTree->Branch("meanCbest", &fMeanCcalc,"meanCcalc/F");
138   fT0OutTree->Branch("multEstimator", &fMultiplicity);
139   //ZDC background
140  fT0OutTree->Branch("zbg", &fZDCbg);
141  fESDpid = new AliESDpid();
142
143   fOutputList->Add(fT0OutTree);
144  
145    PostData(1, fOutputList);
146   
147 }
148
149 //________________________________________________________________________
150 void AliT0CalibAnalysisTask::UserExec(Option_t *) 
151 {
152   // Main loop
153   // Called for each event
154    fVertex= fVertexSPD =  fMeanAC = fMeanA =  fMeanC =  fTrackletSPD=-99999;
155   fMeanACcalc = fMeanAcalc =  fMeanCcalc = fNcont = fNcontTPC = -99999;
156   fMultV0A=fMultV0C=fTimeV0A=fTimeV0C=fMultiplicity = -99999;
157   
158   for (Int_t i=0; i<24; i++) {
159     famp[i] = ftime[i] = -9999;
160     for ( Int_t i0=0; i0<5; i0++) {
161       fRawTime[i][i0] = -9999;
162       if(i0==0) { 
163         fTVDC[i0]=-9999;
164         fOrC[i0]=-9999;
165         fOrA[i0]=-9999;
166       }
167     }
168   }
169   Float_t orA , orC, tvdc;
170   orA = orC = tvdc = -9999;
171   fBC = -9999;
172   // Post output data.
173   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
174   if (!fESD) {
175     printf("ERROR: fESD not available\n");
176     return;
177   }
178
179   if(fESD)
180     {
181       //    if(fIsSelected) 
182       {
183         fZDCbg = kFALSE;
184         /*
185         AliTriggerAnalysis *fTriggerAnalysis;
186         Bool_t  fZDCbgA = fTriggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kZNABG);
187         Bool_t  fZDCbgC =  fTriggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kZNCBG);
188         fZDCbg = (fZDCbgA);
189         cout<<" ZDC back "<<fZDCbgA<<endl;
190         */
191         AliESDVZERO* esdV0 = fESD->GetVZEROData();
192         const AliMultiplicity *mult = fESD->GetMultiplicity();
193         // printf(" VZERO MULT \n");
194         AliESDTZERO* tz= (AliESDTZERO*) fESD->GetESDTZERO();
195
196         const Double32_t *amp, *time, *mean ;
197         const Double32_t *meanbest;
198         fVertex= fVertexSPD =  fMeanAC = fMeanA =  fMeanC =  fTrackletSPD=-99999;
199         fMeanACcalc = fMeanAcalc =  fMeanCcalc = fNcont = fNcontTPC -99999;
200         fMultV0A=fMultV0C=fTimeV0A=fTimeV0C=fMultiplicity = -99999;
201         
202         fSumampA=0;
203         fSumampC=0;
204         fNcont=0;
205         TString triggers = fESD->GetFiredTriggerClasses();
206         fTrigger.SetString(triggers.Data());
207         TString inputtriggers =fESD-> GetHeader()->GetFiredTriggerInputs();
208         // cout<<inputtriggers<<endl;
209         fTriggerinput.SetString(inputtriggers.Data());
210         
211         
212         fEvent=fESD->GetEventNumberInFile();
213         fT0Trigger=fESD->GetT0Trig();
214         
215         fESDtracks=fESD->GetNumberOfTracks();
216         AliTriggerAnalysis *trigAna = new AliTriggerAnalysis;
217         ftimestamp=fESD->GetTimeStamp(); // - 1301217882; 
218         //        printf (" timestamp %d \n", ftimestamp);
219         
220         fMultiplicity =  AliESDtrackCuts::GetReferenceMultiplicity(fESD);
221         //        cout<<" fMultiplicity "<<fMultiplicity<<endl;
222         fOrbit=fESD->GetOrbitNumber();
223         fBC=fESD->GetBunchCrossNumber();
224         fNcont = fESD->GetPrimaryVertexSPD()->GetNContributors();
225         fNcontTPC = fESD->GetPrimaryVertexTPC()->GetNContributors();
226         if(fNcontTPC>=2 )
227           fVertexTPC = fESD->GetPrimaryVertexTPC() ->GetZ();
228         if(fNcont>=0 ) {
229           fVertexSPD = fESD->GetPrimaryVertex()->GetZ();
230         }
231         fTrackletSPD = mult->GetNumberOfTracklets();
232         fClustersSPD = mult->GetNumberOfITSClusters(0); 
233         fMultV0A =  esdV0->GetMTotV0A();
234         fMultV0C =  esdV0->GetMTotV0C();
235         //TOF hit  
236         Int_t ntracksMatchedToTOF = 0; 
237         Int_t ntracks = fESD->GetNumberOfTracks();
238         for(Int_t itrk=0;itrk<ntracks;itrk++){
239           AliESDtrack* track = fESD->GetTrack(itrk);
240           if (!track) {
241             Printf("ERROR: Could not receive track %d", itrk);
242             continue;
243           }
244           //no track selection just TOF hit
245           if (track->IsOn(AliESDtrack::kTOFout)) ntracksMatchedToTOF++;
246         }
247         
248         fTOFtracks = ntracksMatchedToTOF;
249         if(fESDpid){ //get T0_TOF 
250           fESDpid->SetTOFResponse(fESD,AliESDpid::kTOF_T0);
251           fT0tofTrack =(Float_t) (fESDpid->GetTOFResponse().GetStartTime(10.0)); //Get start time "from all tracks 
252           //   printf("@@TOF T0 %f\n",fT0tofTrack);
253         }       
254         
255         
256         fpileup = fESD->IsPileupFromSPD();
257         fTimeV0C = esdV0->GetV0CTime();
258         fTimeV0A = esdV0->GetV0ATime();
259         //        printf(" Vo info \n");      
260         
261         
262         mean = fESD->GetT0TOF();
263         fVertex=fESD->GetT0zVertex();
264         fMeanA = mean[1];
265         fMeanC = mean[2];
266         if (TMath::Abs(fMeanA)<9999 && TMath::Abs(fMeanC)<9999) 
267             fMeanAC = mean[0] ;
268         fSumampC = tz->GetMultC();
269         fSumampA = tz->GetMultA();
270         //    printf("!!!!!!!!fSumampC %f fSumampA %f \n",fSumampC, fSumampA);
271         
272         amp=fESD->GetT0amplitude();
273         time=fESD->GetT0time();
274         for (Int_t i=0; i<24; i++){ 
275           if( time[i]>100 && amp[i]>0.1){
276             famp[i] = amp[i];
277             ftime[i] = time[i];
278           }
279         }
280         
281         //new raw OrA OrC TVDC
282         meanbest = tz->GetT0TOFbest();
283         fMeanAcalc = meanbest[1];
284         fMeanCcalc = meanbest[2];
285         if (TMath::Abs(fMeanAcalc)<9999 && TMath::Abs(fMeanCcalc)<9999) 
286           fMeanACcalc = meanbest[0] ;
287         for (Int_t ii=0; ii<5; ii++){ 
288           orA   = tz->GetOrA(ii);
289           orC   = tz->GetOrC(ii);
290           tvdc  = tz->GetTVDC(ii);
291           //    printf("@@@@@  new signasl %d  %f %f %f \n",ii, orA, orC, tvdc);
292           if(ii==0) {
293             fOrA[ii] = orA;
294             fOrC[ii] = orC;
295             fTVDC[ii] = tvdc;
296             //  cout<<tvdc<<" "<<fT0Trigger<<endl;
297             // if (tvdc>-5 && tvdc<5 && tvdc!=0) 
298             //  cout<<"@@@@@@@@@ "<< tvdc<<" "<<fT0Trigger<<endl;
299           } else {
300             if ( ii>0 && fOrA[ii-1] !=orA)      fOrA[ii] = orA;
301             else fOrA[ii]=-9999;
302             if ( ii>0 && fOrC[ii-1] !=orC)      fOrC[ii] = orC;
303             else fOrC[ii]=-9999;
304           if ( ii>0 && fTVDC[ii-1] != tvdc)     fTVDC[ii] = tvdc;
305           else fTVDC[ii]=-9999;
306           }
307           for (Int_t iii=0; iii<24; iii++)
308             fRawTime[iii][ii] =  tz->GetTimeFull(iii,ii);
309         }
310         fT0pileup[0] = tz->GetPileupFlag();
311         fT0pileup[1] = tz->GetBackgroundFlag();
312         fT0pileup[2] = tz->GetSatellite();
313         
314         AliCentrality *centrality = fESD->GetCentrality();
315         
316         //The centrality function you can use
317         if(centrality) {
318           fcentralityV0M = centrality->GetCentralityPercentile("V0M"); // returns the centrality 
319           fcentralityTRK = centrality->GetCentralityPercentile("TRK"); // returns the centrality 
320           fcentralityZDC = centrality->GetCentralityPercentile("ZEMvsZDC"); // returns the 
321         }
322         AliESDfriend* esdFriend = static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
323         if (esdFriend) {
324           //      cout<<" !!!!!esdFriend "<<endl;
325           AliESDTZEROfriend * t0friend =esdFriend-> GetTZEROfriend();
326           //      cout<<" !!!!t0friend "<<t0friend<<endl;
327           if(t0friend){
328             Double_t  *ftimefr =  t0friend->GetT0timeCorr();
329             Double_t  *fampqtcfr =  t0friend->GetT0ampQTC();
330             Double_t  *fampledfr =  t0friend->GetT0ampLEDminCFD();
331             for (Int_t iii=0; iii<24; iii++){
332               ftimefriend[iii]=ftimefr[iii];
333               fampfriend[iii]=fampqtcfr[iii];
334               fampledfriend[iii]=fampledfr[iii];
335               //              cout<<" friend "<<ftimefriend[iii]<<endl;
336             }
337           } //t0friend
338         } //ESDfriend
339         
340       
341       } //selected
342     } //ESD
343
344   fT0OutTree->Fill();
345     
346   PostData(1, fOutputList);
347   
348 }      
349
350 //________________________________________________________________________
351 void AliT0CalibAnalysisTask::Terminate(Option_t *) 
352 {
353   // Draw result to the screen
354   // Called once at the end of the query
355
356   /* fOutputList = dynamic_cast<TList*> (GetOutputData(1));
357   if (!fOutputList) {
358     printf("ERROR: Output list not available\n");
359     return;
360   }
361     */  
362
363 }