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