]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWG1/TRD/AliTRDcalibration.cxx
save MC pt for Kalman propagator test
[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDcalibration.cxx
1
2 /**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 *                                                                        *
5 * Author: The ALICE Off-line Project.                                    *
6 * Contributors are mentioned in the code where appropriate.              *
7 *                                                                        *
8 * Permission to use, copy, modify and distribute this software and its   *
9 * documentation strictly for non-commercial purposes is hereby granted   *
10 * without fee, provided that the above copyright notice appears in all   *
11 * copies and that both the copyright notice and this permission notice   *
12 * appear in the supporting documentation. The authors make no claims     *
13 * about the suitability of this software for any purpose. It is          *
14 * provided "as is" without express or implied warranty.                  *
15 **************************************************************************/
16
17 /////////////////////////////////////////////////////////////////////////////////
18 //                                                                             
19 // AliTRDcalibration                                                            
20 //                                                                             
21 // Task to run the calibration offline.
22 // Author:
23 //   R. Bailhache (rbailhache@ikf.uni-frankfurt.de, R.Bailhache@gsi.de)
24 //           
25 //////////////////////////////////////////////////////////////////////////////////
26
27
28 #include "Riostream.h"
29 #include "TChain.h"
30 #include "TTree.h"
31 #include "TProfile2D.h"
32 #include "TH2I.h"
33 #include "TH1F.h"
34 #include "TList.h"
35 #include "TMath.h"
36 #include "TCanvas.h"
37 #include "TObject.h"
38 #include "TFile.h"
39 #include "TObjArray.h"
40 #include "TGraph.h"
41 #include "TStyle.h"
42 #include "TLegend.h"
43 #include "TGraphErrors.h"
44
45 #include "AliTRDrecoTask.h"
46 #include "AliAnalysisManager.h"
47
48 #include "AliESDInputHandler.h"
49 #include "AliTRDtrackV1.h"
50 #include "AliTRDseedV1.h"
51 #include "AliTRDcluster.h"
52 #include "info/AliTRDtrackInfo.h"
53 #include "AliTRDcalibDB.h"
54
55 #include "AliTRDCalibraFillHisto.h"
56 #include "AliTRDCalibraFit.h"
57 #include "AliTRDCalibraVdriftLinearFit.h"
58 #include "AliTRDCalibraMode.h"
59 #include "AliTRDCalibraVector.h"
60 #include "./Cal/AliTRDCalPad.h"
61 #include "./Cal/AliTRDCalDet.h"
62
63 #include "AliLog.h"
64
65 #include "AliTRDcalibration.h"
66
67
68 ClassImp(AliTRDcalibration)
69
70 //________________________________________________________________________
71 AliTRDcalibration::AliTRDcalibration() 
72   :AliTRDrecoTask("Calibration", "Calibration on tracks")
73   ,fTrackInfo(0)
74   ,ftrdTrack(0)
75   ,fcl(0)
76   ,fTRDCalibraFillHisto(0)
77   ,fNbTRDTrack(0)
78   ,fNbTRDTrackOffline(0)
79   ,fNbTRDTrackStandalone(0)
80   ,fNbTRDTracklet(0)
81   ,fNbTRDTrackletOffline(0)
82   ,fNbTRDTrackletStandalone(0)
83   ,fNbTimeBin(0x0)
84   ,fNbTimeBinOffline(0x0)
85   ,fNbTimeBinStandalone(0x0)
86   ,fNbClusters(0)
87   ,fNbClustersOffline(0)
88   ,fNbClustersStandalone(0)
89   ,fPHSM(0)
90   ,fCHSM(0)
91   ,fPHSum(0)
92   ,fCHSum(0)
93   ,fDetSum(0)
94   ,fDetSumVector(0)
95   ,fHisto2d(kTRUE)
96   ,fVector2d(kFALSE)
97   ,fVdriftLinear(kTRUE)
98   ,flow(0)
99   ,fhigh(30)
100   ,fNbTimeBins(0)
101   ,ffillZero(kFALSE)
102   ,fnormalizeNbOfCluster(kFALSE)
103   ,fmaxCluster(0)
104   ,fOfflineTracks(kFALSE)
105   ,fStandaloneTracks(kFALSE)
106   ,fCompressPerDetector(kFALSE)
107   ,fGraph(0x0)
108   ,fArrayCalib(0x0)
109   ,fPostProcess(kFALSE)
110 {
111   // Constructor
112   
113   fNRefFigures = 17;
114
115   for(Int_t k = 0; k < 3; k++)
116     {
117       fNz[k]=0;
118       fNrphi[k]=0;
119     }
120
121 }  
122 //________________________________________________________________________
123 AliTRDcalibration::~AliTRDcalibration() 
124 {
125   // Default destructor
126
127   if(fNbTRDTrack) delete fNbTRDTrack;
128   if(fNbTRDTrackOffline) delete fNbTRDTrackOffline;
129   if(fNbTRDTrackStandalone) delete fNbTRDTrackStandalone;
130   if(fNbTRDTracklet) delete fNbTRDTracklet;
131   if(fNbTRDTrackletOffline) delete fNbTRDTrackletOffline;
132   if(fNbTRDTrackletStandalone) delete fNbTRDTrackletStandalone;
133   if(fNbTimeBin) delete fNbTimeBin;
134   if(fNbTimeBinOffline) delete fNbTimeBinOffline;
135   if(fNbTimeBinStandalone) delete fNbTimeBinStandalone;
136   if(fNbClusters) delete fNbClusters;
137   if(fNbClustersOffline) delete fNbClustersOffline;
138   if(fNbClustersStandalone) delete fNbClustersStandalone;
139   if(fPHSM) delete fPHSM;
140   if(fCHSM) delete fCHSM;
141   if(fPHSum) delete fPHSum;
142   if(fCHSum) delete fCHSum;
143   if(fDetSum) delete fDetSum;
144   if(fDetSumVector) delete fDetSumVector;
145   if(fGraph){fGraph->Delete(); delete fGraph;}
146   if(fArrayCalib){fArrayCalib->Delete(); delete fArrayCalib;}
147    
148 }
149 //________________________________________________________________________
150 void AliTRDcalibration::CreateOutputObjects() 
151 {
152   // Create output objects
153
154   OpenFile(0, "RECREATE");
155   
156   // Number of time bins
157   if(fNbTimeBins==0) {
158     AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
159     fNbTimeBins = cal->GetNumberOfTimeBinsDCS();
160     if(fNbTimeBins <= 0){ 
161       AliWarning(Form("No of TimeBins from DB [%d] use default [30]", fNbTimeBins));
162       fNbTimeBins = 30;
163     }
164   }
165   
166   // instance calibration: what to calibrate
167   fTRDCalibraFillHisto = AliTRDCalibraFillHisto::Instance();
168   fTRDCalibraFillHisto->SetHisto2d(fHisto2d); // choose to use histograms
169   fTRDCalibraFillHisto->SetVector2d(fVector2d); // choose to use vectors
170   fTRDCalibraFillHisto->SetCH2dOn();  // choose to calibrate the gain
171   fTRDCalibraFillHisto->SetPH2dOn();  // choose to calibrate the drift velocity
172   fTRDCalibraFillHisto->SetPRF2dOn(); // choose to look at the PRF
173   fTRDCalibraFillHisto->SetLinearFitterOn(fVdriftLinear); // Other possibility vdrift VDRIFT
174   fTRDCalibraFillHisto->SetLinearFitterDebugOn(fVdriftLinear); // Other possibility vdrift
175
176   // segmentation (should be per default the max and add at the end)
177   for(Int_t k = 0; k < 3; k++){
178     if(((fNz[k] != 10) && (fNrphi[k] != 10)) && ((fNz[k] != 100) && (fNrphi[k] != 100))) {
179       fTRDCalibraFillHisto->SetNz(k,fNz[k]);                                    // Mode calibration
180       fTRDCalibraFillHisto->SetNrphi(k,fNrphi[k]);                             // Mode calibration
181     }
182     else {
183       if((fNz[k] == 100) && (fNrphi[k] == 100))  {
184         if(fVector2d) AliInfo("The mode all together is not supported by the vector method");
185         fTRDCalibraFillHisto->SetAllTogether(k);
186       }
187       if((fNz[k] == 10) && (fNrphi[k] == 10))  {
188         if(fVector2d) AliInfo("The mode per supermodule is not supported by the vector method");
189         fTRDCalibraFillHisto->SetPerSuperModule(k);
190       }
191     }
192   }
193
194   // Debug level
195   fTRDCalibraFillHisto->SetDebugLevel(DebugLevel()); //debug stuff
196
197   // Init the stuff
198   fTRDCalibraFillHisto->Init2Dhistos(fNbTimeBins); // initialise the histos
199
200   // cuts
201   fTRDCalibraFillHisto->SetNumberClusters(flow); // At least flow clusters
202   fTRDCalibraFillHisto->SetNumberClustersf(fhigh); // The more fhigh clusters
203   fTRDCalibraFillHisto->SetFillWithZero(ffillZero); // Fill zeros
204   fTRDCalibraFillHisto->SetNormalizeNbOfCluster(fnormalizeNbOfCluster); // For iterations
205
206   // Add them to the container
207   fContainer = new TObjArray();
208   if(fHisto2d) {
209     fContainer->Add(fTRDCalibraFillHisto->GetCH2d()); //TH2I
210     fContainer->Add(fTRDCalibraFillHisto->GetPH2d()); //TProfile2D
211     fContainer->Add(fTRDCalibraFillHisto->GetPRF2d()); //TProfile2D
212   }
213   if(fVdriftLinear) fContainer->Add(fTRDCalibraFillHisto->GetVdriftLinearFit()); // Other drift velocity 
214   if(fVector2d) fContainer->Add(fTRDCalibraFillHisto->GetCalibraVector()); //calibra vector
215       
216   if(DebugLevel()) {
217     
218     // Init the debug histos
219     fNbTRDTrack = new TH1F("TRDTrack","TRDTrack",500,0,500);
220     fNbTRDTrack->Sumw2();
221     fNbTRDTrackOffline = new TH1F("TRDTrackOffline","TRDTrackOffline",500,0,500);
222     fNbTRDTrackOffline->Sumw2();
223     fNbTRDTrackStandalone = new TH1F("TRDTrackStandalone","TRDTrackStandalone",500,0,500);
224     fNbTRDTrackStandalone->Sumw2();
225     //
226     fNbTRDTracklet = new TH1F("TRDTracklet","TRDTracklet",540,0.,540.);
227     fNbTRDTracklet->Sumw2();
228     fNbTRDTrackletOffline = new TH1F("TRDTrackletOffline","TRDTrackletOffline",540,0.,540.);
229     fNbTRDTrackletOffline->Sumw2();
230     fNbTRDTrackletStandalone = new TH1F("TRDTrackletStandalone","TRDTrackletStandalone",540,0.,540.);
231     fNbTRDTrackletStandalone->Sumw2();
232     //
233     fNbTimeBin = new TH1F("TimeBin","TimeBin",35,0,35);
234     fNbTimeBin->Sumw2();
235     fNbTimeBinOffline = new TH1F("TimeBinOffline","TimeBinOffline",35,0,35);
236     fNbTimeBinOffline->Sumw2();
237     fNbTimeBinStandalone = new TH1F("TimeBinStandalone","TimeBinStandalone",35,0,35);
238     fNbTimeBinStandalone->Sumw2();
239     //
240     fNbClusters = new TH1F("NbClusters","",35,0,35);
241     fNbClusters->Sumw2();
242     fNbClustersOffline = new TH1F("NbClustersOffline","",35,0,35);
243     fNbClustersOffline->Sumw2();
244     fNbClustersStandalone = new TH1F("NbClustersStandalone","",35,0,35);
245     fNbClustersStandalone->Sumw2();
246     //
247     fPHSM = new TProfile2D("PH2dSM","Nz10Nrphi10"
248                             ,fNbTimeBins,-0.05,(Double_t)((fNbTimeBins-0.5)/10.0)
249                             ,18,0,18);
250     fPHSM->SetYTitle("Det/pad groups");
251     fPHSM->SetXTitle("time [#mus]");
252     fPHSM->SetZTitle("<PH> [a.u.]");
253     fPHSM->SetStats(0);
254     //
255     fCHSM = new TH2I("CH2dSM","Nz10Nrphi10",50,0,300,18,0,18);
256     fCHSM->SetYTitle("Det/pad groups");
257     fCHSM->SetXTitle("charge deposit [a.u]");
258     fCHSM->SetZTitle("counts");
259     fCHSM->SetStats(0);
260     fCHSM->Sumw2();
261     //
262     fPHSum = new TProfile2D("PH2dSum","Nz100Nrphi100"
263                             ,fNbTimeBins,-0.05,(Double_t)((fNbTimeBins-0.5)/10.0)
264                             ,1,0,1);
265     fPHSum->SetYTitle("Det/pad groups");
266     fPHSum->SetXTitle("time [#mus]");
267     fPHSum->SetZTitle("<PH> [a.u.]");
268     fPHSum->SetStats(0);
269     //
270     fCHSum = new TH2I("CH2dSum","Nz100Nrphi100",50,0,300,1,0,1);
271     fCHSum->SetYTitle("Det/pad groups");
272     fCHSum->SetXTitle("charge deposit [a.u]");
273     fCHSum->SetZTitle("counts");
274     fCHSum->SetStats(0);
275     fCHSum->Sumw2();
276     
277     // Add them
278     fContainer->Add(fNbTRDTrack);
279     fContainer->Add(fNbTRDTrackOffline);
280     fContainer->Add(fNbTRDTrackStandalone);
281     fContainer->Add(fNbTRDTracklet);
282     fContainer->Add(fNbTRDTrackletOffline);
283     fContainer->Add(fNbTRDTrackletStandalone);
284     fContainer->Add(fNbTimeBin);
285     fContainer->Add(fNbTimeBinOffline);
286     fContainer->Add(fNbTimeBinStandalone);
287     fContainer->Add(fNbClusters);
288     fContainer->Add(fNbClustersOffline);
289     fContainer->Add(fNbClustersStandalone);
290     fContainer->Add(fPHSM);
291     fContainer->Add(fCHSM);
292     fContainer->Add(fPHSum);
293     fContainer->Add(fCHSum);
294
295   }
296
297 }
298
299 //________________________________________________________________________
300 void AliTRDcalibration::Exec(Option_t *) 
301 {
302   //
303   // Execute function where the reference data are filled
304   //
305
306   if(!fTracks) return;
307   
308   // In total
309   Int_t nbTrdTracks = 0;
310   // standalone
311   Int_t nbTrdTracksStandalone = 0;
312   // offline
313   Int_t nbTrdTracksOffline = 0;
314   
315
316   //
317   // Loop on track in the event
318   //
319   //printf("Total of %d\n",fTracks->GetEntriesFast());
320   for(Int_t itrk=0; itrk < fTracks->GetEntriesFast(); itrk++){
321     
322     //printf("itrk %d\n",itrk);
323
324     fTrackInfo = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
325     ftrdTrack = fTrackInfo->GetTrack();
326     if(!ftrdTrack) continue;
327
328     nbTrdTracks++;
329   
330     fTRDCalibraFillHisto->UpdateHistogramsV1(ftrdTrack);
331
332     if(DebugLevel()) {
333       
334       Bool_t standalonetracklet = kFALSE;  
335       const AliTRDseedV1 *tracklet = 0x0;
336       //
337       // Loop on tracklet in the event
338       //
339       for(Int_t itr = 0; itr < 6; itr++){
340         //printf("itr %d\n",itr);
341         if(!(tracklet = ftrdTrack->GetTracklet(itr))) continue;
342         if(!tracklet->IsOK()) continue;
343         // standalone
344         if(tracklet->IsStandAlone()) standalonetracklet = kTRUE;
345         Int_t nbclusters = 0;
346         // For PH
347         Double_t phtb[AliTRDseedV1::kNtb];
348         memset(phtb, 0, AliTRDseedV1::kNtb*sizeof(Double_t));
349         // For CH
350         Double_t sum = 0.0;
351         // normalisation
352         Float_t normalisation = 6.67;
353         Int_t detector = 0;
354         Int_t sector = 0;
355         Int_t crossrow = 0;
356         // Check no shared clusters
357         for(int icc=AliTRDseedV1::kNtb; icc<AliTRDseedV1::kNclusters; icc++){
358           if((fcl = tracklet->GetClusters(icc)))  crossrow = 1;
359         }
360         // Loop on clusters
361         for(int ic=0; ic<AliTRDseedV1::kNtb; ic++){
362           if(!(fcl = tracklet->GetClusters(ic))) continue;
363           nbclusters++;
364           Int_t time = fcl->GetPadTime();
365           Float_t ch =  tracklet->GetdQdl(ic);
366           Float_t qcl = TMath::Abs(fcl->GetQ());
367           detector = fcl->GetDetector();          
368           if((time>-1) && (time<fNbTimeBins)) phtb[time]=qcl;
369           sum += ch/normalisation;
370           fNbTimeBin->Fill(time);
371           if(tracklet->IsStandAlone()) fNbTimeBinStandalone->Fill(time);
372           else fNbTimeBinOffline->Fill(time);
373         }
374         sector = AliTRDgeometry::GetSector(detector);
375
376         fNbTRDTracklet->Fill(detector);
377         if(tracklet->IsStandAlone()) fNbTRDTrackletStandalone->Fill(detector);
378         else fNbTRDTrackletOffline->Fill(detector);
379         
380         fNbClusters->Fill(nbclusters);
381         if(tracklet->IsStandAlone())  fNbClustersStandalone->Fill(nbclusters);
382         else  fNbClustersOffline->Fill(nbclusters);
383         
384         if((nbclusters > flow) && (nbclusters < fhigh)){
385           fCHSM->Fill(sum/20.0,sector+0.5);
386           fCHSum->Fill(sum/20.0,0.5);
387           for(int ic=0; ic<fNbTimeBins; ic++){
388             if(ffillZero) {
389               fPHSum->Fill((Double_t)(ic/10.0),0.5,(Double_t)phtb[ic]);
390               fPHSM->Fill((Double_t)(ic/10.0),sector+0.5,(Double_t)phtb[ic]);
391             }
392             else {
393               if(phtb[ic] > 0.0) {
394                 fPHSum->Fill((Double_t)(ic/10.0),0.0,(Double_t)phtb[ic]);
395                 fPHSM->Fill((Double_t)(ic/10.0),sector+0.5,(Double_t)phtb[ic]);
396               }
397             }
398           }
399         }
400       }
401     
402     if(standalonetracklet) nbTrdTracksStandalone++;
403     else nbTrdTracksOffline++;
404     
405     }
406     
407   }
408   
409   if(DebugLevel()) {
410     
411     //Fill Histos
412     fNbTRDTrack->Fill(nbTrdTracks);
413     fNbTRDTrackStandalone->Fill(nbTrdTracksStandalone);
414     fNbTRDTrackOffline->Fill(nbTrdTracksOffline);
415     
416   }
417
418   //printf("Nbof tracks %d\n",nbTrdTracks);
419   
420   // Post output data
421   PostData(0, fContainer);
422
423   //printf("post container\n");
424   
425 }      
426 //________________________________________________________________________
427 void AliTRDcalibration::Terminate(Option_t *) 
428 {
429   // Draw result to the screen
430   // Called once at the end of the query
431
432   //printf("terminate\n");
433
434   if(fTRDCalibraFillHisto) fTRDCalibraFillHisto->DestroyDebugStreamer();
435  
436 }
437 //________________________________________________________
438 Bool_t AliTRDcalibration::GetRefFigure(Int_t ifig)
439 {
440   //
441   // Draw filled histos
442   //
443   
444   gStyle->SetPalette(1);
445   gStyle->SetOptStat(1111);
446   gStyle->SetPadBorderMode(0);
447   gStyle->SetCanvasColor(10);
448   gStyle->SetPadLeftMargin(0.13);
449   gStyle->SetPadRightMargin(0.13);
450
451   if(!fContainer) return kFALSE;
452   
453   switch(ifig){
454   case kNbTrack:{
455     TCanvas *c0 = new TCanvas("c0","c0",10,10,510,510);
456     TLegend *legNbTrack = new TLegend(.7, .7, .98, .98);
457     legNbTrack->SetBorderSize(1);
458     TH1F *h  = 0x0;
459     TH1F *ho = 0x0;
460     TH1F *hs = 0x0;
461     if(!(h = (TH1F *)fContainer->FindObject("TRDTrack"))) break;
462     if(!(ho = (TH1F *)fContainer->FindObject("TRDTrackOffline"))) break;
463     if(!(hs = (TH1F *)fContainer->FindObject("TRDTrackStandalone"))) break;
464     c0->cd();
465     //gPad->SetLogy();
466     gPad->SetGridy();
467     gPad->SetGridx();
468     h->Draw();
469     ho->Draw("same");
470     hs->Draw("same");
471     legNbTrack->AddEntry(h, "all", "p");
472     legNbTrack->AddEntry(ho, "offline", "p");
473     legNbTrack->AddEntry(hs, "standalone", "p");
474     legNbTrack->Draw("same");
475     return kTRUE;
476   }
477   case kNbTracklet:{
478     TLegend *legNbTracklet = new TLegend(.7, .7, .98, .98);
479     legNbTracklet->SetBorderSize(1);
480     TH1F *h = 0x0;
481     TH1F *ho = 0x0;
482     TH1F *hs = 0x0;
483     if(!(h = (TH1F *)fContainer->FindObject("TRDTracklet"))) break;
484     if(!(ho = (TH1F *)fContainer->FindObject("TRDTrackletOffline"))) break;
485     if(!(hs = (TH1F *)fContainer->FindObject("TRDTrackletStandalone"))) break;
486     h->Draw();
487     ho->Draw("same");
488     hs->Draw("same");
489     legNbTracklet->AddEntry(h, "all", "p");
490     legNbTracklet->AddEntry(ho, "offline", "p");
491     legNbTracklet->AddEntry(hs, "standalone", "p");
492     legNbTracklet->Draw("same");
493     gPad->SetLogy();
494     gPad->SetGridy();
495     gPad->SetGridx();
496     return kTRUE;
497   }
498   case kNbTimeBin:{
499     TLegend *legNbTimeBin = new TLegend(.7, .7, .98, .98);
500     legNbTimeBin->SetBorderSize(1);
501     TH1F *h = 0x0;
502     TH1F *ho = 0x0;
503     TH1F *hs = 0x0;
504     if(!(h = (TH1F *)fContainer->FindObject("TimeBin"))) break;
505     if(!(ho = (TH1F *)fContainer->FindObject("TimeBinOffline"))) break;
506     if(!(hs = (TH1F *)fContainer->FindObject("TimeBinStandalone"))) break;
507     h->Draw();
508     ho->Draw("same");
509     hs->Draw("same");
510     legNbTimeBin->AddEntry(h, "all", "p");
511     legNbTimeBin->AddEntry(ho, "offline", "p");
512     legNbTimeBin->AddEntry(hs, "standalone", "p");
513     legNbTimeBin->Draw("same");
514     //gPad->SetLogy();
515     gPad->SetGridy();
516     gPad->SetGridx();
517     return kTRUE;
518   }
519   case kNbClusters:{
520     TLegend *legNbClusters = new TLegend(.7, .7, .98, .98);
521     legNbClusters->SetBorderSize(1);
522     TH1F *h = 0x0;
523     TH1F *ho = 0x0;
524     TH1F *hs = 0x0;
525     if(!(h = (TH1F *)fContainer->FindObject("NbClusters"))) break;
526     if(!(ho = (TH1F *)fContainer->FindObject("NbClustersOffline"))) break;
527     if(!(hs = (TH1F *)fContainer->FindObject("NbClustersStandalone"))) break;
528     h->Draw();
529     ho->Draw("same");
530     hs->Draw("same");
531     legNbClusters->AddEntry(h, "all", "p");
532     legNbClusters->AddEntry(ho, "offline", "p");
533     legNbClusters->AddEntry(hs, "standalone", "p");
534     legNbClusters->Draw("same");
535     gPad->SetLogy();
536     gPad->SetGridy();
537     gPad->SetGridx();
538     return kTRUE;
539   }
540   case kPHSum:{
541     TProfile2D *h = 0x0;
542     if(!(h = (TProfile2D *)fContainer->FindObject("PH2dSum"))) break;
543     TH1D *projh = h->ProjectionX("projh",1,1,"e");
544     projh->Draw();
545     gPad->SetGridy();
546     gPad->SetGridx();
547     return kTRUE;
548   }
549   case kCHSum:{
550     TH2I *h = 0x0;
551     if(!(h = (TH2I *)fContainer->FindObject("CH2dSum"))) break;
552     TH1D *projh = h->ProjectionX("projhh",1,1,"e");
553     projh->Draw();
554     gPad->SetGridy();
555     gPad->SetGridx();
556     return kTRUE;
557   }
558   case kPH2D:{
559     if(!fHisto2d) {
560       AliInfo("Histo was not filled!");
561       break;
562     }
563     TProfile2D *h = 0x0;
564     if(!(h = (TProfile2D *)fContainer->FindObject("PH2d"))) break;
565     h->Draw("lego");
566     return kTRUE;
567   }
568   case kCH2D:{
569     if(!fHisto2d) {
570       AliInfo("Histo was not filled!");
571       break;
572     }
573     TH2I *h = 0x0;
574     if(!(h = (TH2I *)fContainer->FindObject("CH2d"))) break;
575     h->Draw("lego");
576     return kTRUE;
577   }
578   case kPRF2D:{
579     if(!fHisto2d) {
580       AliInfo("Histo was not filled!");
581       break;
582     }
583     TProfile2D *h = 0x0;
584     if(!(h = (TProfile2D *)fContainer->FindObject("PRF2d"))) break;
585     h->Draw("lego");
586     return kTRUE;
587   }
588   case kPH2DVector:{
589     if(!fVector2d) {
590       AliInfo("vector was not filled!");
591       break;
592     }
593     AliTRDCalibraVector *v = 0x0;
594     TGraphErrors *vdet = 0x0; 
595     if(!(v = (AliTRDCalibraVector *)fContainer->FindObject("AliTRDCalibraVector"))) break;
596     Int_t detectormax = -1;
597     Int_t groupmax    = -1;
598     if(!v->FindTheMaxEntries(1,detectormax,groupmax)) break;
599     if(!(vdet = v->ConvertVectorPHTGraphErrors((Int_t)detectormax,groupmax,"plotPH2dVector"))) break;
600     Int_t nbeentries = 0;
601     TH1F *ko = v->CorrectTheError(vdet,nbeentries);
602     ko->Draw();
603     AliInfo(Form("There are %d entries in the detector %d and group %d",nbeentries,detectormax,groupmax));
604     return kTRUE;
605   }
606 case kCH2DVector:{
607     if(!fVector2d) {
608       AliInfo("vector was not filled!");
609       break;
610     }
611     AliTRDCalibraVector *v = 0x0;
612     TH1F *vdet = 0x0; 
613     if(!(v = (AliTRDCalibraVector *)fContainer->FindObject("AliTRDCalibraVector"))) break;
614     Int_t detectormax = -1;
615     Int_t groupmax    = -1;
616     if(!v->FindTheMaxEntries(0,detectormax,groupmax)) break;
617     if(!(vdet = v->ConvertVectorCHHisto((Int_t)detectormax,groupmax,"plotCH2dVector"))) break;
618     vdet->Draw();
619     AliInfo(Form("The detectormax and groupmax are %d and %d",detectormax,groupmax));
620     return kTRUE;
621   }
622   case kPRF2DVector:{
623     if(!fVector2d) {
624       AliInfo("vector was not filled!");
625       break;
626     }
627     AliTRDCalibraVector *v = 0x0;
628     TGraphErrors *vdet = 0x0; 
629     if(!(v = (AliTRDCalibraVector *)fContainer->FindObject("AliTRDCalibraVector"))) break;
630     Int_t detectormax  = -1;
631     Int_t groupmax     = -1;
632     Int_t nbeentries   = 0;
633     if(!v->FindTheMaxEntries(2,detectormax,groupmax)) break;
634     if(!(vdet = v->ConvertVectorPRFTGraphErrors((Int_t)detectormax,groupmax,"plotPRF2dVector"))) break;
635     TH1F *ko = v->CorrectTheError(vdet,nbeentries);
636     ko->Draw();
637     AliInfo(Form("The detectormax and groupmax are %d and %d",detectormax,groupmax));
638     return kTRUE;
639   }
640   case kLinearFitter:{
641     if(!fVdriftLinear) {
642       AliInfo("vdrift linear was not filled!");
643       break;
644     }
645     AliTRDCalibraVdriftLinearFit *h = 0x0;
646     TH2S *hdetector = 0x0; 
647     if(!(h = (AliTRDCalibraVdriftLinearFit *)fContainer->FindObject("AliTRDCalibraVdriftLinearFit"))) break;
648     Double_t entries[540];
649     for(Int_t k = 0; k < 540; k++){
650       entries[k] = 0.0;
651       hdetector = 0x0;
652       if(!(hdetector = (TH2S *)h->GetLinearFitterHisto(k,kFALSE))) continue;
653       entries[k] = hdetector->GetEntries();
654     }
655     Double_t max = -10.0;
656     Double_t detectormax = -1;
657     for(Int_t k = 0; k < 540; k++){
658       if(entries[k] > max) {
659         max = entries[k];
660         detectormax = k;
661       }
662     }
663     hdetector = 0x0;
664     if((TMath::Abs(max) <= 0.001) || (detectormax <0.0) || (detectormax >=540.0)) break;
665     if(!(hdetector = (TH2S *)h->GetLinearFitterHisto((Int_t)detectormax,kFALSE))) break;
666     AliInfo(Form("The detector with the maximum of entries is %d",detectormax));
667     hdetector->Draw();
668     return kTRUE;
669   }
670   case kGainFactor:{
671     if(!fPostProcess){
672       if(!PostProcess()) break;
673     }
674     TGraph *fgain = (TGraph *) fGraph->At(0);
675     if(!fgain) break;
676     fgain->Draw("ALP");
677     return kTRUE;
678   }
679   case kVdriftT0Factor:{
680     if(!fPostProcess){
681       if(!PostProcess()) break;
682     }
683     TCanvas *c = new TCanvas("c","c",10,10,510,510);
684     c->Divide(2,1);
685     TGraph *fvd = (TGraph *) fGraph->At(1);
686     if(fvd){
687       c->cd(1);
688       fvd->Draw("ALP");
689     } 
690     TGraph *ft0 = (TGraph *) fGraph->At(2);
691     if(ft0){
692       c->cd(2);
693       ft0->Draw("ALP");
694     } 
695     return kTRUE;
696   }
697   case kVdriftLorentzAngleFactor:{
698     if(!fVdriftLinear) {
699       AliInfo("vdrift linear was not filled!");
700       break;
701     }
702     if(!fPostProcess){
703       if(!PostProcess()) break;
704     }
705     TCanvas *c = new TCanvas("c","c",10,10,510,510);
706     c->Divide(2,1);
707     TGraph *fvdl = (TGraph *) fGraph->At(3);
708     if(fvdl){
709       c->cd(1);
710       fvdl->Draw("ALP");
711     } 
712     TGraph *flal = (TGraph *) fGraph->At(4);
713     if(flal){
714       c->cd(2);
715       flal->Draw("ALP");
716     } 
717     return kTRUE;
718   }
719   case kPRFFactor:{
720     if(!fPostProcess){
721       if(!PostProcess()) break;
722     }
723     TGraph *fprf = (TGraph *) fGraph->At(5);
724     if(!fprf) break;
725     fprf->Draw("ALP");
726     return kTRUE;
727   }
728   }
729   
730   return kFALSE;
731   
732 }
733 //________________________________________________________________________
734 Bool_t AliTRDcalibration::PostProcess()
735 {
736   // 
737   // Fit the filled histos
738   // Put the calibration object in fArrayCalib
739   // 0 and 1 AliTRDCalDet and AliTRDCalPad gain
740   // 2 and 3 AliTRDCalDet and AliTRDCalPad vdrift PH
741   // 4 and 5 AliTRDCalDet and AliTRDCalPad timeoffset PH
742   // 6 AliTRDCalPad PRF
743   // 7 and 8 AliTRDCalDet and AliTRDCalPad vdrift second
744   // 9 and 10 AliTRDCalDet and AliTRDCalPad lorentz angle second
745   //
746
747   if(!fArrayCalib){
748     fArrayCalib = new TObjArray(11);
749     fArrayCalib->SetOwner();
750   }
751   else {
752     delete fArrayCalib;
753     PostProcess();
754   }
755   
756   if(!fGraph){
757     fGraph = new TObjArray(6);
758     fGraph->SetOwner();
759   }
760   else {
761     delete fGraph;
762     PostProcess();
763   }
764
765   Bool_t storage[3] = {kFALSE,kFALSE,kFALSE};
766
767   // Objects for fitting
768   AliTRDCalibraFit *calibra = AliTRDCalibraFit::Instance();
769   calibra->SetDebugLevel(2); // 0 rien, 1 fitvoir, 2 debug files, 3 one detector  
770   
771   // Take the stuff
772   if (!fContainer) {
773     Printf("ERROR: list not available");
774     return kFALSE;
775   }
776
777   if(fHisto2d && fVector2d) AliInfo("We will only look at histos. Set fHisto2d off if you don't want");
778   AliTRDCalibraVector *calibraVector = 0x0;
779   if(fVector2d) calibraVector = (AliTRDCalibraVector *) fContainer->FindObject("CalibraVector");
780   //
781   // GAIN TH2I
782   //
783   Bool_t pass = kFALSE; 
784   AliTRDCalibraVector *vvect = 0x0;
785   if(fHisto2d) {
786     TH2I *histogain = (TH2I *) fContainer->FindObject("CH2d");  
787     if(histogain) {
788       histogain->SetDirectory(0);
789       calibra->SetMinEntries(20); 
790       if(fCompressPerDetector){
791         if(AddStatsPerDetector(histogain)) pass = calibra->AnalyseCH(fCHSum);
792       }
793       else pass = calibra->AnalyseCH(histogain);
794     }
795   }
796   else {
797     if(fVector2d && calibraVector) {
798       calibra->SetMinEntries(20); 
799       if(fCompressPerDetector){
800         if(!(vvect = calibraVector->AddStatsPerDetectorCH())) return kFALSE;
801         pass = calibra->AnalyseCH(vvect);
802       }
803       else pass = calibra->AnalyseCH(calibraVector);
804     }
805   }
806   
807   if(pass)
808     {
809       Int_t nbtg = 6*4*18*((Int_t) ((AliTRDCalibraMode *)calibra->GetCalibraMode())->GetDetChamb0(0))
810         + 6*  18*((Int_t) ((AliTRDCalibraMode *)calibra->GetCalibraMode())->GetDetChamb2(0));
811       Int_t nbfit = calibra->GetNumberFit();  //Number of fits
812       Int_t nbE   = calibra->GetNumberEnt();  //Number of detector mit entries
813       // enough statistics
814       if ((nbtg >                  0) && 
815           (nbfit        >= 0.001*nbE))
816         {
817           // create the cal objects
818           calibra->PutMeanValueOtherVectorFit(1,kTRUE); 
819           TObjArray object           = calibra->GetVectorFit();
820           AliTRDCalDet *objgaindet   = calibra->CreateDetObjectGain(&object);
821           TObject *objgainpad        = calibra->CreatePadObjectGain();
822           // store
823           fArrayCalib->AddAt(objgaindet,0);
824           fArrayCalib->AddAt(objgainpad,1);
825           storage[0] = kTRUE;
826           // Make graph
827           TGraph *graph = 0x0;
828           if(FillGraphIndex(&object,graph)){ 
829             fGraph->AddAt(graph,0);
830           }
831         }//if(enough statistics?)
832       calibra->ResetVectorFit();
833     }
834   else return kFALSE;
835   
836   //
837   // VDRIFT average pulse height
838   //
839   pass = kFALSE; 
840   if(fHisto2d) {
841     TProfile2D *histodriftvelocity = (TProfile2D *) fContainer->FindObject("PH2d");  
842     if(histodriftvelocity) {
843       histodriftvelocity->SetDirectory(0);  
844       calibra->SetMinEntries(20*20);  
845       if(fCompressPerDetector){
846         if(AddStatsPerDetector(histodriftvelocity,1)) {
847           pass = calibra->AnalysePH(fDetSumVector);
848         }
849       }
850       else pass = calibra->AnalysePH(histodriftvelocity);
851     }
852   }
853   else {
854     if(fVector2d && calibraVector) {
855       calibra->SetMinEntries(20*20);  
856       if(fCompressPerDetector){
857         if(!(vvect = calibraVector->AddStatsPerDetectorPH())) return kFALSE;
858         pass = calibra->AnalysePH(vvect);
859       }
860       else pass = calibra->AnalysePH(calibraVector);  
861     }
862   }
863
864   if(pass) {
865     Int_t nbtg = 6*4*18*((Int_t) ((AliTRDCalibraMode *)calibra->GetCalibraMode())->GetDetChamb0(1))
866       + 6*  18*((Int_t) ((AliTRDCalibraMode *)calibra->GetCalibraMode())->GetDetChamb2(1));
867     Int_t nbfit  = calibra->GetNumberFit();
868     Int_t nbE    = calibra->GetNumberEnt();
869     // enough statistics
870     if ((nbtg >                  0) && 
871         (nbfit        >= 0.001*nbE))
872       {
873         // create the cal objects
874         calibra->PutMeanValueOtherVectorFit(1,kTRUE);
875         calibra->PutMeanValueOtherVectorFit2(1,kTRUE); 
876         TObjArray object  = calibra->GetVectorFit();
877         AliTRDCalDet *objdriftvelocitydet = calibra->CreateDetObjectVdrift(&object,kTRUE);
878         TObject *objdriftvelocitypad      = calibra->CreatePadObjectVdrift();
879         TObjArray objectt          = calibra->GetVectorFit2();
880         AliTRDCalDet *objtime0det  = calibra->CreateDetObjectT0(&object,kTRUE);
881         TObject *objtime0pad       = calibra->CreatePadObjectT0();
882         // store
883         fArrayCalib->AddAt(objdriftvelocitydet,2);
884         fArrayCalib->AddAt(objdriftvelocitypad,3);
885         //
886         fArrayCalib->AddAt(objtime0det,4);
887         fArrayCalib->AddAt(objtime0pad,5);
888         // Make graph
889         TGraph *graph = 0x0;
890         if(FillGraphIndex(&object,graph)){ 
891           fGraph->AddAt(graph,1);
892         }
893         TGraph *graphh = 0x0;
894         if(FillGraphIndex(&objectt,graphh)){ 
895           fGraph->AddAt(graphh,2);
896         }
897       }//if(enough statistics)
898     calibra->ResetVectorFit();
899   }
900   else return kFALSE;
901   
902   //
903   // PRF
904   //
905   pass = kFALSE; 
906   if(fHisto2d) {
907     TProfile2D *histoprf = (TProfile2D *) fContainer->FindObject("PRF2d");
908     if (histoprf) {
909       histoprf->SetDirectory(0);  
910       calibra->SetMinEntries(600); 
911       if(fCompressPerDetector){
912         if(AddStatsPerDetector(histoprf,2)) pass = calibra->AnalysePRFMarianFit(fDetSumVector);
913       }
914       else pass = calibra->AnalysePRFMarianFit(histoprf);
915     }
916   }
917   else {
918     if(fVector2d && calibraVector) {
919       calibra->SetMinEntries(600);  
920       if(fCompressPerDetector){
921         if(!(vvect =calibraVector->AddStatsPerDetectorPRF())) return kFALSE;
922         pass = calibra->AnalysePRFMarianFit(vvect);
923       }
924       else pass = calibra->AnalysePRFMarianFit(calibraVector);  
925     }
926   }
927   
928   if(pass){
929     Int_t nbtg = 6*4*18*((Int_t) ((AliTRDCalibraMode *)calibra->GetCalibraMode())->GetDetChamb0(2))
930       + 6*  18*((Int_t) ((AliTRDCalibraMode *)calibra->GetCalibraMode())->GetDetChamb2(2));
931     Int_t nbfit        = calibra->GetNumberFit();
932     Int_t nbE          = calibra->GetNumberEnt();
933     // enough statistics
934     if ((nbtg >                  0) && 
935         (nbfit        >= 0.001*nbE)) {
936       TObjArray object            = calibra->GetVectorFit();
937       TObject *objPRFpad          = calibra->CreatePadObjectPRF(&object);
938       // store
939       fArrayCalib->AddAt(objPRFpad,6);
940       // Make graph
941       TGraph *graph = 0x0;
942       if(FillGraphIndex(&object,graph)){ 
943         fGraph->AddAt(graph,5);
944       }
945     }
946     calibra->ResetVectorFit();
947   }
948   else return kFALSE;
949   
950   //
951   // VDRIFT linear fit 
952   //
953   AliTRDCalibraVdriftLinearFit *vlinearfit = (AliTRDCalibraVdriftLinearFit *) fContainer->FindObject("LinearVdriftFit"); 
954   if (vlinearfit) {
955     calibra->SetMinEntries(20*20);     
956     if(calibra->AnalyseLinearFitters(vlinearfit)) {
957       
958       Int_t nbtg = 6*4*18*((Int_t) ((AliTRDCalibraMode *)calibra->GetCalibraMode())->GetDetChamb0(2))
959         + 6*  18*((Int_t) ((AliTRDCalibraMode *)calibra->GetCalibraMode())->GetDetChamb2(2));
960       Int_t nbfit  = calibra->GetNumberFit();
961       Int_t nbE    = calibra->GetNumberEnt();
962       // enough statistics
963       if ((nbtg >                  0) && 
964           (nbfit        >= 0.001*nbE))
965         {
966           // create the cal objects
967           calibra->PutMeanValueOtherVectorFit(1,kTRUE);
968           calibra->PutMeanValueOtherVectorFit2(1,kTRUE); 
969           TObjArray object  = calibra->GetVectorFit();
970           AliTRDCalDet *objdriftvelocitydet = calibra->CreateDetObjectVdrift(&object,kTRUE);
971           TObject *objdriftvelocitypad      = calibra->CreatePadObjectVdrift();
972           TObjArray objectt          = calibra->GetVectorFit2();
973           AliTRDCalDet *objtime0det  = calibra->CreateDetObjectT0(&object,kTRUE);
974           TObject *objtime0pad       = calibra->CreatePadObjectT0();
975           // store dummy
976           fArrayCalib->AddAt(objdriftvelocitydet,7);
977           fArrayCalib->AddAt(objdriftvelocitypad,8);
978           //
979           fArrayCalib->AddAt(objtime0det,9);
980           fArrayCalib->AddAt(objtime0pad,10);
981           // Make graph
982           TGraph *graph = 0x0;
983           if(FillGraphIndex(&object,graph)){ 
984             fGraph->AddAt(graph,3);
985           }
986           TGraph *graphh = 0x0;
987           if(FillGraphIndex(&objectt,graphh)){ 
988             fGraph->AddAt(graphh,4);
989           }
990         }//if(enough statistics)
991     }// if fit
992     calibra->ResetVectorFit();
993   }
994   else return kFALSE;
995   
996   fPostProcess = kTRUE;
997   
998   return kTRUE;
999   
1000 }
1001
1002 //________________________________________________________________________
1003 Bool_t AliTRDcalibration::FillGraphIndex(const TObjArray *vectora,TGraph *graph) const
1004 {
1005   //
1006   // Fill one value (the first one) per detector
1007   //
1008
1009   Int_t loop = (Int_t) vectora->GetEntriesFast();
1010   if(loop != 540) {
1011     AliInfo("The Vector Fit is not complete!");
1012     return kFALSE;
1013   }
1014   
1015   Double_t x[540];
1016   Double_t y[540];
1017   for (Int_t k = 0; k < loop; k++) {
1018     if(!vectora->At(k)){
1019       x[k] = -1.0;
1020       y[k] = -1.0;
1021       continue;
1022     }
1023     x[k]  = ((AliTRDCalibraFit::AliTRDFitInfo *) vectora->At(k))->GetDetector();
1024     y[k]  = ((Float_t *)((AliTRDCalibraFit::AliTRDFitInfo *) vectora->At(k))->GetCoef())[0];
1025   }
1026
1027   if(!graph) graph = new TGraph(540,&x[0],&y[0]);
1028   else{ 
1029     graph->~TGraph();
1030     new(graph) TGraph(540,&x[0],&y[0]);
1031   }
1032
1033   return kTRUE;
1034
1035 }
1036 //________________________________________________________________________
1037 Bool_t AliTRDcalibration::AddStatsPerDetector(const TH2I *ch) 
1038 {
1039   //
1040   // Add statistic per detector
1041   //
1042   
1043   AliTRDCalibraMode calibMode = AliTRDCalibraMode();
1044   const char *name = ch->GetTitle();
1045   if(!SetNzFromTObject(name,0,&calibMode)) return 0x0;
1046   if(!SetNrphiFromTObject(name,0,&calibMode)) return 0x0;
1047   if(((calibMode.GetNz(0) == 100) && (calibMode.GetNrphi(0) == 100)) || ((calibMode.GetNz(0) == 10) && (calibMode.GetNrphi(0) == 10))) return kFALSE;
1048
1049   Int_t    nybins  = ch->GetNbinsY();// groups number
1050   Int_t    nxbins  = ch->GetNbinsX();// number of bins X
1051   TAxis   *xaxis   = ch->GetXaxis();
1052   Double_t lowedge  = xaxis->GetBinLowEdge(1);
1053   Double_t upedge   = xaxis->GetBinUpEdge(nxbins);
1054
1055   // number per chamber 2
1056   calibMode.ModePadCalibration(2,0);
1057   calibMode.ModePadFragmentation(0,2,0,0);
1058   calibMode.SetDetChamb2(0);
1059   Int_t perChamber2 = (Int_t) calibMode.GetDetChamb2(0);
1060
1061   // number per other chamber
1062   calibMode.ModePadCalibration(0,0);
1063   calibMode.ModePadFragmentation(0,0,0,0);
1064   calibMode.SetDetChamb0(0);
1065   Int_t perChamber0 = (Int_t) calibMode.GetDetChamb0(0);
1066
1067   if(nybins != (6*18*perChamber2+6*4*18*perChamber0)) return kFALSE;
1068
1069   // Create Histo
1070   TString nname((const char *)ch->GetName());
1071   nname  += "PerDetector";
1072   TString title("Nz");
1073   title += 0;
1074   title += "Nrphi";
1075   title += 0;
1076   if(!fCHSum) fCHSum = new TH2I((const char *)nname,(const char *)title
1077                                 ,nxbins,lowedge,upedge,540,0,540);
1078   else{ 
1079     fCHSum->~TH2I();
1080     new(fCHSum) TH2I((const Char_t *) nname,(const char *)title
1081                      ,nxbins,lowedge,upedge,540,0,540);
1082   }
1083   fCHSum->SetYTitle("Detector number");
1084   fCHSum->SetXTitle("charge deposit [a.u]");
1085   fCHSum->SetZTitle("counts");
1086   fCHSum->SetStats(0);
1087   fCHSum->Sumw2();
1088
1089   Int_t counter = 0;
1090   
1091   for(Int_t det = 0; det < 540; det++){
1092
1093     Int_t numberofgroup = 0;
1094     if(AliTRDgeometry::GetStack(det) == 2) numberofgroup = perChamber2;
1095     else numberofgroup = perChamber0;
1096     TH1I *projch = (TH1I *) ch->ProjectionX("projch",counter+1,counter+numberofgroup,(Option_t *)"e");
1097     projch->SetDirectory(0);
1098        
1099     for(Int_t nx = 0; nx <= nxbins; nx++) {
1100       fCHSum->SetBinContent(nx,det+1,projch->GetBinContent(nx));
1101       fCHSum->SetBinError(nx,det+1,projch->GetBinError(nx));
1102     }
1103
1104     counter += numberofgroup;
1105     
1106     delete projch;
1107
1108   }
1109
1110   return kTRUE;
1111
1112 }
1113 //_____________________________________________________________________________________________________________________
1114 Bool_t AliTRDcalibration::AddStatsPerDetector(const TProfile2D *ph,Int_t i)
1115 {
1116   //
1117   // Add statistic per detector
1118   //
1119
1120   AliTRDCalibraMode calibMode = AliTRDCalibraMode();
1121   const char *name = ph->GetTitle();
1122   //printf("name %s\n",name);
1123   if(!SetNzFromTObject(name,0,&calibMode)) return kFALSE;
1124   if(!SetNrphiFromTObject(name,0,&calibMode)) return kFALSE;
1125   if(((calibMode.GetNz(0) == 100) && (calibMode.GetNrphi(0) == 100)) || ((calibMode.GetNz(0) == 10) && (calibMode.GetNrphi(0) == 10))) return kFALSE;
1126   //printf("Found mode Mz %d, Nrphi %d\n",calibMode.GetNz(0),calibMode.GetNrphi(0));  
1127
1128
1129   Int_t    nybins  = ph->GetNbinsY();// groups number
1130   Int_t    nxbins  = ph->GetNbinsX();// number of bins X
1131   TAxis   *xaxis = ph->GetXaxis();
1132   Double_t lowedge  = xaxis->GetBinLowEdge(1);
1133   Double_t upedge   = xaxis->GetBinUpEdge(nxbins);
1134
1135   // number per chamber 2
1136   calibMode.ModePadCalibration(2,0);
1137   calibMode.ModePadFragmentation(0,2,0,0);
1138   calibMode.SetDetChamb2(0);
1139   Int_t perChamber2 = (Int_t) calibMode.GetDetChamb2(0);
1140
1141   // number per other chamber
1142   calibMode.ModePadCalibration(0,0);
1143   calibMode.ModePadFragmentation(0,0,0,0);
1144   calibMode.SetDetChamb0(0);
1145   Int_t perChamber0 = (Int_t) calibMode.GetDetChamb0(0);
1146
1147   if(nybins != (6*18*perChamber2+6*4*18*perChamber0)) return kFALSE;
1148   
1149   // Create calvector 
1150   TString nbname((const char *)ph->GetName());
1151   nbname  += "PerDetectorVector";
1152   if(!fDetSumVector) fDetSumVector = new AliTRDCalibraVector();
1153   else{ 
1154     fDetSumVector->~AliTRDCalibraVector();
1155     new(fDetSumVector) AliTRDCalibraVector();
1156   }
1157   if(i==1){
1158     fDetSumVector->SetTimeMax(nxbins);
1159   }
1160   if(i==2){
1161     fDetSumVector->SetNumberBinPRF(nxbins);
1162     fDetSumVector->SetPRFRange(TMath::Abs(lowedge));
1163   }
1164   fDetSumVector->SetDetCha0(i,1);
1165   fDetSumVector->SetDetCha2(i,1);
1166   fDetSumVector->SetNzNrphi(i,0,0);
1167   if(i==2) {
1168     Int_t nbg = GetNumberOfGroupsPRF((const char *)name);
1169     fDetSumVector->SetNbGroupPRF(nbg);
1170   }
1171
1172   // Create Histo
1173   TString nname((const char *)ph->GetName());
1174   nname  += "PerDetector";
1175   TString title("Nz");
1176   title += 0;
1177   title += "Nrphi";
1178   title += 0;
1179   if(!fDetSum) fDetSum = new TH2D((const char *)nname,(const Char_t *) title
1180                                 ,nxbins,lowedge,upedge,540,0,540);
1181   else{ 
1182     fDetSum->~TH2D();
1183     new(fDetSum) TH2D((const Char_t *) nname,(const Char_t *) title
1184                      ,nxbins,lowedge,upedge,540,0,540);
1185   }
1186   fDetSum->SetYTitle("Detector number");
1187   fDetSum->SetXTitle(xaxis->GetTitle());
1188   fDetSum->SetStats(0);
1189   
1190   Int_t counter = 0;
1191
1192   for(Int_t det = 0; det < 540; det++){
1193
1194     Int_t numberofgroup = 0;
1195     if(AliTRDgeometry::GetStack(det) == 2) numberofgroup = perChamber2;
1196     else numberofgroup = perChamber0;
1197     
1198     for(Int_t nx = 1; nx <= nxbins; nx++) {
1199       
1200       Double_t entries = 0.0;
1201       Double_t sumw2 = 0.0;
1202       Double_t sumw = 0.0;
1203
1204       for(Int_t k = counter+1; k <= (counter+numberofgroup); k++){
1205         Int_t  binnumber = ph->GetBin(nx,k);
1206         entries += ph->GetBinEntries(binnumber);
1207         sumw2 += (ph->GetBinError(binnumber)*ph->GetBinError(binnumber)+ph->GetBinContent(binnumber)*ph->GetBinContent(binnumber))*ph->GetBinEntries(binnumber);
1208         sumw += ph->GetBinContent(binnumber)*ph->GetBinEntries(binnumber);
1209       }
1210
1211       Double_t mean = 0.0;
1212       if(entries > 0.0) mean = sumw/entries;
1213       Double_t squaremean = 0.0;
1214       if(entries > 0.0) squaremean = sumw2/entries;
1215       Double_t errorf = squaremean - mean*mean;
1216       Double_t error = 0.0;
1217       if(entries > 0.0) error = TMath::Sqrt(TMath::Abs(errorf)/entries);
1218       
1219       fDetSum->SetBinContent(nx,det+1,mean);
1220       fDetSum->SetBinError(nx,det+1,error);
1221
1222       if(i==1) fDetSumVector->FillVectorPH(det,0,nx-1,(Int_t)entries,(Float_t)mean,(Float_t)squaremean);
1223       if(i==2) fDetSumVector->FillVectorPRF(det,0,nx-1,(Int_t)entries,(Float_t)mean,(Float_t)squaremean);
1224       
1225     }
1226     
1227     counter += numberofgroup;
1228
1229   }
1230
1231   return kTRUE;
1232
1233   
1234 }
1235 //_____________________________________________________________________________
1236 Bool_t AliTRDcalibration::SetNrphiFromTObject(const char *name, Int_t i, AliTRDCalibraMode *calibMode) const
1237 {
1238   //
1239   // Set the granularity from object
1240   //  
1241   
1242   const Char_t *patternrphi0 = "Nrphi0";
1243   const Char_t *patternrphi1 = "Nrphi1";
1244   const Char_t *patternrphi2 = "Nrphi2";
1245   const Char_t *patternrphi3 = "Nrphi3";
1246   const Char_t *patternrphi4 = "Nrphi4";
1247   const Char_t *patternrphi5 = "Nrphi5";
1248   const Char_t *patternrphi6 = "Nrphi6";
1249
1250   
1251   const Char_t *patternrphi10 = "Nrphi10";
1252   const Char_t *patternrphi100 = "Nrphi100";
1253   const Char_t *patternz10 = "Nz10";
1254   const Char_t *patternz100 = "Nz100";
1255
1256   // Nrphi mode
1257   if ((strstr(name,patternrphi100)) && (strstr(name,patternz100))) {
1258     calibMode->SetAllTogether(i);
1259     return kTRUE;
1260   }
1261   if ((strstr(name,patternrphi10)) && (strstr(name,patternz10))) {
1262     calibMode->SetPerSuperModule(i);
1263     return kTRUE;
1264   }
1265   
1266   if (strstr(name,patternrphi0)) {
1267     calibMode->SetNrphi(i ,0);
1268     return kTRUE;
1269   }
1270   if (strstr(name,patternrphi1)) {
1271     calibMode->SetNrphi(i, 1);
1272     return kTRUE;
1273   }
1274   if (strstr(name,patternrphi2)) {
1275     calibMode->SetNrphi(i, 2);
1276     return kTRUE;
1277   }
1278   if (strstr(name,patternrphi3)) {
1279     calibMode->SetNrphi(i, 3);
1280     return kTRUE;
1281   }
1282   if (strstr(name,patternrphi4)) {
1283     calibMode->SetNrphi(i, 4);
1284     return kTRUE;
1285   }
1286   if (strstr(name,patternrphi5)) {
1287     calibMode->SetNrphi(i, 5);
1288     return kTRUE;
1289   }
1290   if (strstr(name,patternrphi6)) {
1291     calibMode->SetNrphi(i, 6);
1292     return kTRUE;
1293   }
1294   
1295   calibMode->SetNrphi(i ,0);
1296   return kFALSE;
1297   
1298 }
1299 //_____________________________________________________________________________
1300 Bool_t AliTRDcalibration::SetNzFromTObject(const char *name, Int_t i, AliTRDCalibraMode *calibMode) const
1301 {
1302   //
1303   // Set fNz[i] of the AliTRDCalibraFit::Instance()
1304   // corresponding to the given TObject
1305   //
1306
1307   // Some patterns
1308   const Char_t *patternz0    = "Nz0";
1309   const Char_t *patternz1    = "Nz1";
1310   const Char_t *patternz2    = "Nz2";
1311   const Char_t *patternz3    = "Nz3";
1312   const Char_t *patternz4    = "Nz4";
1313
1314   const Char_t *patternrphi10 = "Nrphi10";
1315   const Char_t *patternrphi100 = "Nrphi100";
1316   const Char_t *patternz10 = "Nz10";
1317   const Char_t *patternz100 = "Nz100";
1318
1319   if ((strstr(name,patternrphi100)) && (strstr(name,patternz100))) {
1320     calibMode->SetAllTogether(i);
1321     return kTRUE;
1322   }
1323   if ((strstr(name,patternrphi10)) && (strstr(name,patternz10))) {
1324     calibMode->SetPerSuperModule(i);
1325     return kTRUE;
1326   }
1327   if (strstr(name,patternz0)) {
1328     calibMode->SetNz(i, 0);
1329     return kTRUE;
1330   }
1331   if (strstr(name,patternz1)) {
1332     calibMode->SetNz(i ,1);
1333     return kTRUE;
1334   }
1335   if (strstr(name,patternz2)) {
1336     calibMode->SetNz(i ,2);
1337     return kTRUE;
1338   }
1339   if (strstr(name,patternz3)) {
1340     calibMode->SetNz(i ,3);
1341     return kTRUE;  
1342   }
1343   if (strstr(name,patternz4)) {
1344     calibMode->SetNz(i ,4);
1345     return kTRUE;
1346   }
1347  
1348   calibMode->SetNz(i ,0);
1349   return kFALSE;
1350 }
1351 //____________________________________________________________________________________________________
1352 Int_t AliTRDcalibration::GetNumberOfGroupsPRF(const char* nametitle) const
1353 {
1354   //
1355   // Get numberofgroupsprf
1356   //
1357   
1358   // Some patterns
1359   const Char_t *pattern0 = "Ngp0";
1360   const Char_t *pattern1 = "Ngp1";
1361   const Char_t *pattern2 = "Ngp2";
1362   const Char_t *pattern3 = "Ngp3";
1363   const Char_t *pattern4 = "Ngp4";
1364   const Char_t *pattern5 = "Ngp5";
1365   const Char_t *pattern6 = "Ngp6";
1366
1367   // Nrphi mode
1368   if (strstr(nametitle,pattern0)) {
1369     return 0;
1370   }
1371   if (strstr(nametitle,pattern1)) {
1372     return 1;
1373   }
1374   if (strstr(nametitle,pattern2)) {
1375     return 2;
1376   }
1377   if (strstr(nametitle,pattern3)) {
1378     return 3;
1379   }
1380   if (strstr(nametitle,pattern4)) {
1381     return 4;
1382   }
1383   if (strstr(nametitle,pattern5)) {
1384     return 5;
1385   }
1386   if (strstr(nametitle,pattern6)){
1387     return 6;
1388   }
1389   else return -1;
1390 }