/[PAMELA software]/PamelaLevel2/src/PamLevel2.cpp
ViewVC logotype

Contents of /PamelaLevel2/src/PamLevel2.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.42 - (show annotations) (download)
Tue Aug 7 15:37:12 2007 UTC (17 years, 6 months ago) by pam-fi
Branch: MAIN
Changes since 1.41: +6 -4 lines
sorttracks modified again

1 #include <PamLevel2.h>
2 //--------------------------------------
3 //
4 //
5 //--------------------------------------
6 /**
7 * Default constructor
8 */
9 PamTrack::PamTrack(){
10 trk_track = 0;
11 calo_track = 0;
12 tof_track = 0;
13 candeleteobj = 0;
14 };
15 //--------------------------------------
16 //
17 //
18 //--------------------------------------
19 /**
20 * Constructor
21 */
22 PamTrack::PamTrack(TrkTrack* t, CaloTrkVar* c, ToFTrkVar* o){
23
24 trk_track = 0;
25 calo_track = 0;
26 tof_track = 0;
27 // if(t)trk_track = new TrkTrack(*t);
28 // if(c)calo_track = new CaloTrkVar(*c);
29 // if(o)tof_track = new ToFTrkVar(*o);
30 if(t)trk_track = t;
31 if(c)calo_track = c;
32 if(o)tof_track = o;
33 candeleteobj = 0;
34
35 };
36 PamTrack::PamTrack(const PamTrack& track){
37
38 TrkTrack *t = track.trk_track;
39 CaloTrkVar *c = track.calo_track;
40 ToFTrkVar *o = track.tof_track;
41
42 trk_track = 0;
43 calo_track = 0;
44 tof_track = 0;
45 if(t)trk_track = new TrkTrack(*t);
46 if(c)calo_track = new CaloTrkVar(*c);
47 if(o)tof_track = new ToFTrkVar(*o);
48 candeleteobj = 1;
49
50 }
51 void PamTrack::Clear(){
52
53 // cout << "PamTrack::Clear() "<<candeleteobj<<endl;
54 if(candeleteobj){
55 if(trk_track) trk_track->TrkTrack::Clear();
56 if(calo_track) calo_track->CaloTrkVar::Clear();//???
57 if(tof_track) tof_track->ToFTrkVar::Clear();//???
58 }else{
59 trk_track = 0;
60 calo_track = 0;
61 tof_track = 0;
62 }
63 }
64 void PamTrack::Delete(){
65 // cout << "PamTrack::Delete() "<<candeleteobj<<endl;
66 if(candeleteobj){
67 if(trk_track) {
68 trk_track->TrkTrack::Clear();
69 delete trk_track;
70 }
71 if(calo_track){
72 calo_track->CaloTrkVar::Clear();//???
73 delete calo_track;
74 }
75 if(tof_track){
76 tof_track->ToFTrkVar::Clear();//???
77 delete tof_track;
78 }
79 }else{
80 Clear();
81 }
82 }
83 //--------------------------------------
84 //
85 //
86 //--------------------------------------
87 /**
88 * Default Constructor
89 */
90 PamLevel2::PamLevel2(){
91 Initialize();
92 };
93
94 /**
95 * Constructor with given dir, list and detectors
96 */
97 PamLevel2::PamLevel2(TString ddir,TString list,TString detlist){
98 Initialize();
99 GetPamTree(GetListOfLevel2Files(ddir,list),detlist);
100 GetRunTree(GetListOfLevel2Files(ddir,list));
101 };
102
103 /**
104 * Constructor with given dir and list
105 */
106 PamLevel2::PamLevel2(TString ddir,TString list){
107 Initialize();
108 GetPamTree(GetListOfLevel2Files(ddir,list),"");
109 GetRunTree(GetListOfLevel2Files(ddir,list));
110 };
111
112
113 void PamLevel2::Initialize(){
114
115 // trk2_obj = TrkLevel2::GetTrkLevel2();
116 // trk1_obj = TrkLevel1::GetTrkLevel1();
117 // trkh_obj = TrkHough::GetTrkHough();
118 // calo1_obj = CaloLevel1::GetCaloLevel1();
119 // calo2_obj = CaloLevel2::GetCaloLevel2();
120 // tof_obj = ToFLevel2::GetToFLevel2();
121 // trig_obj = TrigLevel2::GetTrigLevel2();
122 // s4_obj = S4Level2::GetS4Level2();
123 // nd_obj = NDLevel2::GetNDLevel2();
124 // ac_obj = AcLevel2::GetAcLevel2();
125 // orb_obj = OrbitalInfo::GetOrbitalInfo();
126
127 h0_obj = 0;
128 trk0_obj = 0;
129 calo0_obj = 0;
130
131 trk2_obj = 0;
132 trk1_obj = 0;
133 trkh_obj = 0;
134 calo1_obj = 0;
135 calo2_obj = 0;
136 tof_obj = 0;
137 trig_obj = 0;
138 s4_obj = 0;
139 nd_obj = 0;
140 ac_obj = 0;
141 orb_obj = 0;
142
143 run_obj = 0;//new GL_RUN();
144 soft_obj = 0;// Emiliano
145 irun = -1LL;
146 runfirstentry = 0LL;
147 runlastentry = 0LL;
148
149 l0_file = NULL;
150 l0_tree = NULL;
151 iroot = -1;
152 dbc = 0;
153
154 prevshift = 0;
155
156 run_tree = NULL;
157 run_tree_clone = NULL;
158 sel_tree = NULL;
159 sel_tree_clone = NULL;
160
161 irunentry = -1LL;
162 pam_tree = NULL;
163 for(Int_t i=0; i<NCLONES; i++ )pam_tree_clone[i]=NULL;
164
165 host = "mysql://localhost/pamelaprod";
166 user = "anonymous";
167 psw = "";
168 const char *pamdbhost=gSystem->Getenv("PAM_DBHOST");
169 const char *pamdbuser=gSystem->Getenv("PAM_DBUSER");
170 const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW");
171 if ( !pamdbhost ) pamdbhost = "";
172 if ( !pamdbuser ) pamdbuser = "";
173 if ( !pamdbpsw ) pamdbpsw = "";
174 if ( strcmp(pamdbhost,"") ) host = pamdbhost;
175 if ( strcmp(pamdbuser,"") ) user = pamdbuser;
176 if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw;
177
178
179
180 // sorted_tracks = 0;//new TRefArray();
181
182 CAL0 = false;
183 CAL1 = true;
184 CAL2 = true;
185 TRK2 = true;
186 TRK1 = false;
187 TRKh = false;
188 TRKh = false;
189 TRG = true;
190 TOF = true;
191 TOF0 = false;
192 S4 = true;
193 ND = true;
194 AC = true;
195 ORB = true;
196
197 RUN = true;
198
199 SELLI = -1;
200
201 tsorted=0;
202 timage=0;
203
204 howtosort = "+CAL+TOF";
205 //howtosort = "+TOF";
206 sortthr = 100.;
207
208 };
209 /**
210 * Delete the event (act as Dtor)
211 */
212 void PamLevel2::Delete(){
213
214 if(run_obj)delete run_obj;
215 if(soft_obj)delete soft_obj; //Emiliano
216
217 // cout << "void PamLevel2::Clear()"<<endl;
218 if(h0_obj) delete h0_obj;
219 if(trk0_obj) delete trk0_obj;
220 if(calo0_obj) delete calo0_obj;
221 if(trk1_obj) delete trk1_obj;
222 if(trk2_obj) delete trk2_obj;
223 if(trkh_obj) delete trkh_obj;
224 if(calo1_obj)delete calo1_obj;
225 if(calo2_obj)delete calo2_obj;
226 if(tof_obj) delete tof_obj;
227 if(trig_obj) delete trig_obj;
228 if(s4_obj) delete s4_obj;
229 if(nd_obj) delete nd_obj;
230 if(ac_obj) delete ac_obj;
231 if(orb_obj) delete orb_obj;
232
233 if(tsorted){
234 tsorted->Delete();
235 delete tsorted;
236 }
237 if(timage){
238 timage->Delete();
239 delete timage;
240 }
241
242 if(dbc){
243 dbc->Close();
244 delete dbc;
245 }
246
247 if(l0_file)l0_file->Close();
248 // if(pam_tree)pam_tree->Delete();;
249
250 if ( pam_tree ){
251 //
252 // we have first to find which chains we have to delete, then delete the main chain and only after delete the friend chains. Any other order causes a segfault...
253 //
254 TList *temp = pam_tree->GetListOfFriends();
255 TList *contents = new TList; // create chain friend list
256 contents->SetOwner();
257 TIter next(temp);
258 TChain *questo = 0;
259 while ( (questo = (TChain*) next()) ){
260 TString name = questo->GetName();
261 contents->Add((TChain*)gROOT->FindObject(name.Data()));// add object to the list
262 };
263 //
264 // deleting the main chain
265 //
266 pam_tree->Delete();
267 //
268 // deleting the friends...
269 //
270 TIter next2(contents);
271 TChain *questa = 0;
272 while ( questa = (TChain*)next2() ){
273 TString name = questa->GetName();
274 questa->Delete();
275 questa=NULL;
276 };
277 //
278 };
279 pam_tree = NULL;
280
281 if(run_tree)run_tree->Delete();;
282 if(sel_tree)sel_tree->Delete();;
283 for(Int_t i=0; i<NCLONES; i++ )if(pam_tree_clone[i])pam_tree_clone[i]->Delete();;
284 if(run_tree_clone)run_tree_clone->Delete();;
285 if(sel_tree_clone)sel_tree_clone->Delete();;
286
287
288 };
289
290 /**
291 * Clear the event (NB! does not deallocate objects)
292 */
293 void PamLevel2::Clear(){
294
295 // cout << "void PamLevel2::Clear()"<<endl;
296
297 //
298 // This method is called once for every entry but RunInfo and SoftInfo do not change until the next run so we cannot clear them here unless we don't
299 // want to load them for each event even if they are the same...
300 //
301 // if(run_obj)delete run_obj;
302 // if(run_obj) run_obj->Clear(); // Emiliano: Do not deallocate run_obj here, it will give segmentation fault! call clear instead
303 // if(soft_obj) soft_obj->Clear();
304
305 if(h0_obj) h0_obj->Clear();
306 // if(trk0_obj) trk0_obj->Clear();
307 if(trk1_obj) trk1_obj->Clear();
308 if(trk2_obj) trk2_obj->Clear();
309 if(trkh_obj) trkh_obj->Clear();
310 if(calo0_obj) calo0_obj->Clear();
311 if(calo1_obj)calo1_obj->Clear();
312 if(calo2_obj)calo2_obj->Clear();
313 if(tof_obj) tof_obj->Clear();
314 if(trig_obj) trig_obj->Clear();
315 if(s4_obj) s4_obj->Clear();
316 if(nd_obj) nd_obj->Clear();
317 if(ac_obj) ac_obj->Clear();
318 if(orb_obj) orb_obj->Clear();
319
320 // if(sorted_tracks)sorted_tracks->Clear();
321 // sorted_tracks.Clear();
322
323 if(tsorted){
324 tsorted->Delete();
325 }
326 if(timage){
327 timage->Delete();
328 }
329 };
330
331 void PamLevel2::Reset(){
332 //
333 // First of all clear everything
334 //
335 Clear();
336 //
337 // close and reset chains and pointers
338 //
339 if ( pam_tree ){
340 //
341 // we have first to find which chains we have to delete, then delete the main chain and only after delete the friend chains. Any other order causes a segfault...
342 //
343 TList *temp = pam_tree->GetListOfFriends();
344 TList *contents = new TList; // create chain friend list
345 contents->SetOwner();
346 TIter next(temp);
347 TChain *questo = 0;
348 while ( (questo = (TChain*) next()) ){
349 TString name = questo->GetName();
350 contents->Add((TChain*)gROOT->FindObject(name.Data()));// add object to the list
351 };
352 //
353 // deleting the main chain
354 //
355 pam_tree->Delete();
356 //
357 // deleting the friends...
358 //
359 TIter next2(contents);
360 TChain *questa = 0;
361 while ( questa = (TChain*)next2() ){
362 TString name = questa->GetName();
363 questa->Delete();
364 questa=NULL;
365 };
366 //
367 };
368 pam_tree = NULL;
369 //
370 if(run_tree)run_tree->Delete();;
371 run_tree = NULL;
372 if(sel_tree)sel_tree->Delete();;
373 sel_tree = NULL;
374 //
375 // Close file
376 //
377 if(l0_file)l0_file->Close("R");
378 l0_file = NULL;
379 //
380 h0_obj = 0;
381 trk0_obj = 0;
382 calo0_obj = 0;
383 //
384 trk2_obj = 0;
385 trk1_obj = 0;
386 trkh_obj = 0;
387 calo1_obj = 0;
388 calo2_obj = 0;
389 tof_obj = 0;
390 trig_obj = 0;
391 s4_obj = 0;
392 nd_obj = 0;
393 ac_obj = 0;
394 orb_obj = 0;
395 //
396 // Reset run pointers
397 //
398 run_obj = 0;//new GL_RUN();
399 soft_obj = 0;// Emiliano
400 irun = -1;
401 runfirstentry = 0ULL;
402 runlastentry = 0ULL;
403 //
404 };
405
406 Bool_t PamLevel2::IsGood(){
407 Bool_t goodev=true;
408 // if(trk2_obj && trk2_obj->UnpackError() != 0 ) goodev = false;
409 if(calo2_obj && calo2_obj->good != 1) goodev = false;
410 if(tof_obj && tof_obj->unpackError != 0) goodev = false;
411 if(trig_obj && trig_obj->unpackError != 0) goodev = false;
412 if(s4_obj && s4_obj->unpackError != 0) goodev = false;
413 if(nd_obj && nd_obj->unpackError != 0) goodev = false;
414 if(ac_obj && ac_obj->unpackError != 255) goodev = false;
415 // if(orb_obj)
416 return goodev;
417 };
418
419 //--------------------------------------
420 //
421 //
422 //--------------------------------------
423 void *PamLevel2::GetPointerTo(const char* c ){
424
425 TString objname = c;
426
427 if(!objname.CompareTo("TrkLevel1")) {
428 if(!trk1_obj){
429 trk1_obj = new TrkLevel1();
430 trk1_obj->Set();
431 }
432 return &trk1_obj;
433 };
434 if(!objname.CompareTo("TrkLevel2")) {
435 if(!trk2_obj){
436 trk2_obj = new TrkLevel2();
437 trk2_obj->Set();
438 }
439 return &trk2_obj;
440 };
441 if(!objname.CompareTo("TrkHough")) {
442 if(!trkh_obj) trkh_obj = new TrkHough();
443 return &trkh_obj;
444 };
445 if(!objname.CompareTo("CaloLevel1")) {
446 if(!calo1_obj) calo1_obj = new CaloLevel1();
447 return &calo1_obj;
448 };
449 if(!objname.CompareTo("CaloLevel2")) {
450 if(!calo2_obj){
451 calo2_obj = new CaloLevel2();
452 calo2_obj->Set();
453 };
454 return &calo2_obj;
455 };
456 if(!objname.CompareTo("ToFLevel2")) {
457 if(!tof_obj){
458 tof_obj = new ToFLevel2();
459 tof_obj->Set();
460 }
461 return &tof_obj;
462 };
463 if(!objname.CompareTo("TrigLevel2")) {
464 if(!trig_obj) trig_obj = new TrigLevel2();
465 return &trig_obj;
466 };
467 if(!objname.CompareTo("S4Level2")) {
468 if(!s4_obj) s4_obj = new S4Level2();
469 return &s4_obj;
470 };
471 if(!objname.CompareTo("NDLevel2")) {
472 if(!nd_obj) nd_obj = new NDLevel2();
473 return &nd_obj;
474 };
475 if(!objname.CompareTo("AcLevel2")) {
476 if(!ac_obj) ac_obj = new AcLevel2();
477 return &ac_obj;
478 };
479 if(!objname.CompareTo("OrbitalInfo")){
480 if(!orb_obj) orb_obj = new OrbitalInfo();
481 return &orb_obj;
482 };
483
484 if(!objname.CompareTo("RunInfo"))return &run_obj;
485
486 if(!objname.CompareTo("SoftInfo"))return &soft_obj; // Emiliano
487
488 return NULL;
489 };
490 //--------------------------------------
491 //
492 //
493 //--------------------------------------
494 /**
495 * Retrieves the calorimeter track matching the seqno-th tracker stored track.
496 * (If seqno = -1 retrieves the self-trigger calorimeter track)
497 */
498 CaloTrkVar *PamLevel2::GetCaloStoredTrack(int seqno){
499
500 if( !calo2_obj )return 0;
501
502 if( calo2_obj->CaloLevel2::ntrk()==0 ){
503 cout << "PamLevel2::GetCaloStoredTrack(int) : requested tracker SeqNo "<< seqno <<" but no Calorimeter tracks are stored"<<endl;
504 return NULL;
505 };
506
507 CaloTrkVar *c = 0;
508 Int_t it_calo=0;
509
510 do{
511 c = calo2_obj->CaloLevel2::GetCaloTrkVar(it_calo);
512 it_calo++;
513 } while( c && seqno != c->trkseqno && it_calo < calo2_obj->CaloLevel2::ntrk());
514
515 if(!c || seqno != c->trkseqno){
516 c = 0;
517 if(seqno!=-1)cout << "PamLevel2::GetCaloStoredTrack(int) : requested tracker SeqNo "<< seqno <<" does not match Calorimeter stored tracks"<<endl;
518 };
519 return c;
520
521 };
522 //--------------------------------------
523 //
524 //
525 //--------------------------------------
526 /**
527 * Retrieves the ToF track matching the seqno-th tracker stored track.
528 * (If seqno = -1 retrieves the tracker-independent tof track)
529 */
530 ToFTrkVar *PamLevel2::GetToFStoredTrack(int seqno){
531
532 if( !tof_obj )return 0;
533
534 if( tof_obj->ToFLevel2::ntrk()==0 ){
535 cout << "PamLevel2::GetToFStoredTrack(int) : requested tracker SeqNo "<< seqno <<" but no ToF tracks are stored"<<endl;
536 return NULL;
537 };
538
539 ToFTrkVar *c = 0;
540 Int_t it_tof=0;
541
542 do{
543 c = tof_obj->ToFLevel2::GetToFTrkVar(it_tof);
544 it_tof++;
545 } while( c && seqno != c->trkseqno && it_tof < tof_obj->ToFLevel2::ntrk());
546
547 if(!c || seqno != c->trkseqno){
548 c = 0;
549 if(seqno!=-1)cout << "PamLevel2::GetToFStoredTrack(int) : requested tracker SeqNo "<< seqno <<" does not match ToF stored tracks"<<endl;
550 };
551 return c;
552
553 };
554
555 //--------------------------------------
556 //
557 //
558 //--------------------------------------
559 /**
560 * Give the pamela track associated to a tracker track, retrieving related calorimeter and tof track information.
561 */
562 PamTrack* PamLevel2::GetPamTrackAlong(TrkTrack* t){
563
564 cout <<"PamLevel2::GetPamTrackAlong(TrkTrack* t) **obsolete** "<<endl;
565 cout <<"(if you use it, remember to delete the PamTrack object)"<<endl;
566
567 CaloTrkVar *c = 0;
568 ToFTrkVar *o = 0;
569
570 if(CAL2) c = GetCaloStoredTrack(t->GetSeqNo());
571 if(TOF) o = GetToFStoredTrack(t->GetSeqNo());
572
573 // if(t && c && o)track = new PamTrack(t,c,o);
574 PamTrack *track = new PamTrack(t,c,o);
575
576 return track;
577
578 };
579 //--------------------------------------
580 //
581 //
582 //--------------------------------------
583 /**
584 * Retrieves the it-th stored track.
585 * It override TrkLevel2::GetTrack(int it).
586 * @param itrk Track number, ranging from 0 to GetNTracks().
587 */
588
589 PamTrack* PamLevel2::GetStoredTrack(Int_t itrk){
590
591 cout <<"PamLevel2::GetStoredTrack(Int_t itrk) **to-be-updated** "<<endl;
592 cout <<"for the moment, better use separately the methods: TrkLevel2::GetStoredTrack(seqno) CaloLevel2::GetCaloTrkVar(Int_t notrack) ToFLevel2::GetToFTrkVar(Int_t notrack)"<<endl;
593 cout <<"(if you use it, remember to delete the PamTrack object)"<<endl;
594 PamTrack *track = 0;
595
596 if( itrk >=0 && itrk < trk2_obj->TrkLevel2::ntrk() ){
597
598 TrkTrack *t = trk2_obj->TrkLevel2::GetStoredTrack(itrk);
599 track = GetPamTrackAlong(t);
600
601
602
603 }else{
604 cout << "PamLevel2::GetStoredTrack(int) : tracker track SeqNo "<< itrk <<" does not exist (GetNTracks() = "<<trk2_obj->TrkLevel2::GetNTracks()<<")"<<endl;
605 };
606
607 return track;
608
609 }
610 //--------------------------------------
611 //
612
613 /**
614 * Sort physical (tracker) tracks and stores them in the TRefArray (of TrkTrack objects) which pointer is PamLevel2::sorted_tracks. Left here as backward compatibility method.
615 **/
616 void PamLevel2::SortTracks(TString how){
617 printf(" WARNING! obsolete, use SortTracks() and SetSortingMethod(TString) instead! \n Setting sorting method to %s \n",how.Data());
618 howtosort = how;
619 SortTracks();
620 };
621
622 //
623 //--------------------------------------
624 /**
625 * Sort physical (tracker) tracks and stores them in the TRefArray (of TrkTrack objects) which pointer is PamLevel2::sorted_tracks.
626 * @param how String to set the sorting cryterium (es: "CAL" or "TRK+CAL+TOF" ecc...).
627 * Sorting cryteria:
628 * TRK: lower chi**2
629 * CAL: lower Y spatial residual on the first calorimeter plane
630 * TOF: bigger numebr of hit PMTs along the track, on S12 S21 S32 (where paddles are along the Y axis).
631 * The default sorting cryterium is "TOF+CAL".
632 *
633 * The total number of physical tracks is always given by GetNTracks() and the it-th physical track can be retrieved by means of the methods GetTrack(int it) and GetTrack(int it, TString how).
634 */
635 void PamLevel2::SortTracks(){
636 TString how = howtosort;
637
638 // cout <<" PamLevel2::SortTracks(TString how) "<<endl;
639 if( !trk2_obj ){
640 cout << "void PamLevel2::SortTracks(): TrkLevel2 not loaded !!!";
641 return;
642 };
643 // cout << "call SortTracs() "<<endl;
644 //Save current Object count
645 Int_t ObjectNumber = TProcessID::GetObjectCount();
646
647 // cout << "ObjectNumber "<<ObjectNumber <<endl;
648
649 // if(!sorted_tracks)sorted_tracks = new TRefArray();
650 // sorted_tracks->Clear();
651 // sorted_tracks.Clear();
652
653 if(!tsorted)tsorted = new TClonesArray("PamTrack",trk2_obj->GetNTracks());
654 tsorted->Delete();
655 TClonesArray &ttsorted = *tsorted;
656 if(!timage)timage = new TClonesArray("PamTrack",trk2_obj->GetNTracks());
657 timage->Delete();
658 TClonesArray &ttimage = *timage;
659
660
661
662 //--------------------------------------------------
663 // retrieve sorting method
664 //--------------------------------------------------
665 Bool_t use_TRK = how.Contains("TRK", TString::kIgnoreCase);
666 Bool_t use_CAL = how.Contains("CAL", TString::kIgnoreCase);
667 Bool_t use_TOF = how.Contains("TOF", TString::kIgnoreCase);
668
669 if( !CAL2 && use_CAL) use_CAL = false;
670 if( !TOF && use_TOF) use_TOF = false;
671
672 if( !TRK2 ){
673 // cout << "SortTracks() : without tracker does not work!!! (not yet)" << endl;
674 return;
675 };
676
677 // cout << "use_CAL "<<use_CAL<<" use_TOF "<<use_TOF<<" use_TRK "<<use_TRK <<endl;
678
679
680 //--------------------------------------------------
681 // loop over "physical" tracks sorted by the tracker
682 //--------------------------------------------------
683 for(Int_t i=0; i < trk2_obj->TrkLevel2::GetNTracks(); i++){
684
685 TrkTrack *ts = 0;
686 CaloTrkVar *cs = 0;
687 ToFTrkVar *os = 0;
688
689 // get tracker tracks
690 TrkTrack *tp = trk2_obj->TrkLevel2::GetTrack(i); //tracker
691 CaloTrkVar *cp = GetCaloStoredTrack(tp->GetSeqNo());
692 ToFTrkVar *op = GetToFStoredTrack(tp->GetSeqNo());
693
694 TrkTrack *ti = 0; //tracker (image)
695 CaloTrkVar *ci = 0;
696 ToFTrkVar *oi = 0;
697 // cout << "trk track n. "<<i << " "<<hex<< tp <<dec<< endl;
698 // if track has an image, check image selection
699 if(tp->HasImage()){
700
701 ti = trk2_obj->TrkLevel2::GetTrackImage(i); //tracker (image)
702 ci = GetCaloStoredTrack(ti->GetSeqNo());
703 oi = GetToFStoredTrack(ti->GetSeqNo());
704
705 // cout << "its image "<<i << " "<<hex<< ti <<dec<< endl;
706
707 //assign starting scores
708 Int_t tp_score = 0; //main track sorted by the tracker
709 Int_t ti_score = 0; //image track
710
711 // -----------------------------------------------------------------------------------------
712 // calorimeter check
713 // -----------------------------------------------------------------------------------------
714 // check the Y spatial residual on the first calorimeter plane
715 // (cut on calorimeter variables from Emiliano)
716 if( use_CAL && !calo2_obj ){
717 cout << "void PamLevel2::SortTracks(): howtosort= "<<how<<" but CaloLevel2 not loaded !!!";
718 return;
719 };
720 if(
721 use_CAL &&
722 // calo2_obj->npcfit[1] > 5 && //no. of fit planes on Y view
723 // calo2_obj->varcfit[1] < 1000. && //fit variance on Y view
724 cp && ci &&
725 true){
726
727
728 // Float_t resy_p = cp->tbar[0][1] - calo2_obj->cbar[0][1]; if(resy_p < 0)resy_p= - resy_p;
729 // Float_t resy_i = ci->tbar[0][1] - calo2_obj->cbar[0][1]; if(resy_i < 0)resy_i= - resy_i;
730
731 // if(resy_p <= resy_i) tp_score++;
732 // else ti_score++;
733
734
735 if( cp->npresh > ci->npresh) tp_score++;
736 if( cp->npresh < ci->npresh) ti_score++;
737 else ;//niente
738
739 // cout << "CALO "<<tp_score<<ti_score<<endl;
740
741 };
742 // -----------------------------------------------------------------------------------------
743 // TOF check
744 // -----------------------------------------------------------------------------------------
745 // check the number of hit pmts along the track
746 // on S12 S21 and S32, where paddles are parallel to Y axis
747 if( use_TOF && !tof_obj ){
748 cout << "void PamLevel2::SortTracks(): howtosort= "<<how<<" but ToFLevel2 not loaded !!!";
749 return;
750 };
751 if( use_TOF && op && oi ){
752
753 //
754 Float_t sen = 0.;
755 for (Int_t ih=0; ih < op->npmtadc; ih++){
756 Int_t pl = tof_obj->GetPlaneIndex( (op->pmtadc).At(ih) );
757 if ( pl == 2 || pl == 3 || pl == 4 || pl == 5 ) sen += (op->dedx).At(ih);
758 };
759 for (Int_t ih=0; ih < oi->npmtadc; ih++){
760 Int_t pl = tof_obj->GetPlaneIndex( (oi->pmtadc).At(ih) );
761 if ( pl == 2 || pl == 3 || pl == 4 || pl == 5 ) sen += (oi->dedx).At(ih);
762 };
763 //
764 if ( sen >= sortthr ){
765 //printf(" IS A NUCLEUS! en = %f \n",sen);
766 //
767 // is a nucleus use a different algorithm
768 //
769 Int_t nz = 6; Float_t zin[6]; // << define TOF z-coordinates
770 for(Int_t ip=0; ip<nz; ip++)
771 zin[ip] = tof_obj->ToFLevel2::GetZTOF(tof_obj->ToFLevel2::GetToFPlaneID(ip)); // << read ToF plane z-coordinates
772 Trajectory *tr = new Trajectory(nz,zin);
773 //
774 Int_t nphit_p =0;
775 Int_t nphit_i =0;
776 Float_t enhit_p = 0.;
777 Float_t enhit_i = 0.;
778 //
779 for (Int_t ih=0; ih < op->npmtadc; ih++){
780 Int_t pl = tof_obj->GetPlaneIndex( (op->pmtadc).At(ih) );
781 if(pl == 1 || pl == 2 || pl == 5){
782 nphit_p++;
783 enhit_p += (op->dedx).At(ih);
784 };
785 };
786 //
787 tp->DoTrack2(tr);
788 //
789 if ( fabs(tr->y[0]-oi->ytofpos[0]) < 2. ){
790 for (Int_t ih=0; ih < op->npmtadc; ih++){
791 Int_t pl = tof_obj->GetPlaneIndex( (op->pmtadc).At(ih) );
792 if(pl == 0){
793 nphit_p++;
794 enhit_p += (op->dedx).At(ih);
795 };
796 };
797 };
798 if ( fabs(tr->y[3]-oi->ytofpos[1]) < 2. ){
799 for (Int_t ih=0; ih < op->npmtadc; ih++){
800 Int_t pl = tof_obj->GetPlaneIndex( (op->pmtadc).At(ih) );
801 if(pl == 3){
802 nphit_p++;
803 enhit_p += (op->dedx).At(ih);
804 };
805 };
806 };
807 if ( fabs(tr->y[4]-oi->ytofpos[2]) < 2. ){
808 for (Int_t ih=0; ih < op->npmtadc; ih++){
809 Int_t pl = tof_obj->GetPlaneIndex( (op->pmtadc).At(ih) );
810 if(pl == 4){
811 nphit_p++;
812 enhit_p += (op->dedx).At(ih);
813 };
814 };
815 };
816
817 for (Int_t ih=0; ih < oi->npmtadc; ih++){
818 Int_t pl = tof_obj->GetPlaneIndex( (oi->pmtadc).At(ih) );
819 if(pl == 1 || pl == 2 || pl == 5){
820 nphit_i++;
821 enhit_i += (op->dedx).At(ih);
822 };
823 };
824 //
825 ti->DoTrack2(tr);
826 //
827 if ( fabs(tr->y[0]-oi->ytofpos[0]) < 2. ){
828 for (Int_t ih=0; ih < oi->npmtadc; ih++){
829 Int_t pl = tof_obj->GetPlaneIndex( (oi->pmtadc).At(ih) );
830 if(pl == 0){
831 nphit_i++;
832 enhit_i += (op->dedx).At(ih);
833 };
834 };
835 };
836 if ( fabs(tr->y[3]-oi->ytofpos[1]) < 2. ){
837 for (Int_t ih=0; ih < oi->npmtadc; ih++){
838 Int_t pl = tof_obj->GetPlaneIndex( (oi->pmtadc).At(ih) );
839 if(pl == 3){
840 nphit_i++;
841 enhit_i += (op->dedx).At(ih);
842 };
843 };
844 };
845 if ( fabs(tr->y[4]-oi->ytofpos[2]) < 2. ){
846 for (Int_t ih=0; ih < oi->npmtadc; ih++){
847 Int_t pl = tof_obj->GetPlaneIndex( (oi->pmtadc).At(ih) );
848 if(pl == 4){
849 nphit_i++;
850 enhit_i += (op->dedx).At(ih);
851 };
852 };
853 };
854
855
856 if(
857 use_TOF &&
858 (nphit_p+nphit_i) !=0 &&
859 true){
860
861 // printf(" seqno %i nphit_p %i nphit_i %i enhit_p %f enhit_i %f \n",trk2_obj->TrkLevel2::GetSeqNo(i),nphit_p,nphit_i,enhit_p,enhit_i);
862 // printf(" score p %i score i %i \n",tp_score,ti_score);
863 // if( enhit_p > enhit_i ) tp_score++;
864 // if( nphit_p >= nphit_i && enhit_p > enhit_i ) tp_score++;
865 if ( nphit_p > nphit_i ) tp_score++;
866 if ( nphit_p < nphit_i ) ti_score++;
867 if ( nphit_p == nphit_i ){
868 if ( enhit_p > enhit_i ) tp_score++;
869 else ti_score++;
870 };
871 // printf(" dopo score p %i score i %i \n",tp_score,ti_score);
872 };
873 delete tr;
874 //
875 } else {
876 //
877 // NOT a NUCLEUS
878 //
879 //printf(" NOT a NUCLEUS! en = %f \n",sen);
880
881 Int_t nphit_p =0;
882 Int_t nphit_i =0;
883
884
885 /* cout << "track: npmtadc "<< op->npmtadc << endl;
886 cout << "track: npmttdc "<< op->npmttdc << endl;
887 cout << "image: npmtadc "<< oi->npmtadc << endl;
888 cout << "image: npmttdc "<< oi->npmttdc << endl;*/
889
890 // for (Int_t ih=0; ih < op->npmtadc; ih++){
891 // Int_t pl = tof_obj->GetPlaneIndex( (op->pmtadc).At(ih) );
892 // if(pl == 1 || pl == 2 || pl == 5)nphit_p++;
893 // };
894
895 // for (Int_t ih=0; ih < oi->npmtadc; ih++){
896 // Int_t pl = tof_obj->GetPlaneIndex( (oi->pmtadc).At(ih) );
897 // if(pl == 1 || pl == 2 || pl == 5)nphit_i++;
898 // };
899 // --- modified to count tdc signals (more efficient?)
900 // --- and to implement check on tdcflag
901 for (Int_t ih=0; ih < op->npmttdc; ih++){
902 Int_t pl = tof_obj->GetPlaneIndex( (op->pmttdc).At(ih) );
903 // if( (op->tdcflag).At(ih)==0 && (pl == 1 || pl == 2 || pl == 5) )nphit_p++;
904 if( (op->tdcflag).At(ih)==0 )nphit_p++;
905 };
906
907 for (Int_t ih=0; ih < oi->npmttdc; ih++){
908 Int_t pl = tof_obj->GetPlaneIndex( (oi->pmttdc).At(ih) );
909 // if( (oi->tdcflag).At(ih)==0 && (pl == 1 || pl == 2 || pl == 5) )nphit_i++;
910 if( (oi->tdcflag).At(ih)==0 )nphit_i++;
911 };
912
913 if(
914 (nphit_p+nphit_i) !=0 &&
915 true){
916
917 if ( nphit_p > nphit_i) tp_score++;
918 else if( nphit_p < nphit_i) ti_score++;
919 else ;//niente
920 };
921 };
922 // cout << "TOF "<<tp_score<<ti_score<<endl;
923 };
924
925
926 // if(tp_score == ti_score) use_TRK = true;
927
928
929 // -----------------------------------------------------------------------------------------
930 // tracker check
931 // -----------------------------------------------------------------------------------------
932 // chi**2 difference is not always large enough to distinguish among
933 // the real track and its image.
934 // Tracker check will be applied always when calorimeter and tof information is ambiguous.
935 // *** MODIFIED ON AUGUST 2007 ***
936 if(use_TRK){
937 // if( tp->chi2 > 0 && tp->chi2 < ti->chi2 ) tp_score++ ;
938 // else if( ti->chi2 > 0 && ti->chi2 < tp->chi2 ) ti_score++ ;
939
940 // CHECK 1 : number of points along X
941 if ( tp->GetNX() >= ti->GetNX() )tp_score++ ;
942 if ( tp->GetNX() <= ti->GetNX() )ti_score++ ;
943 // CHECK 2 : number of points along Y
944 if ( tp->GetNY() >= ti->GetNY() )tp_score++ ;
945 if ( tp->GetNY() <= ti->GetNY() )ti_score++ ;
946 // CHECK 3 : chi**2 along X
947 // if( tp->GetChi2X() > 0 && (tp->GetChi2X() < ti->GetChi2X() || ti->GetChi2X() <=0) ) tp_score++ ;
948 // if( ti->GetChi2X() > 0 && (ti->GetChi2X() < tp->GetChi2X() || tp->GetChi2X() <=0) ) ti_score++ ;
949 // CHECK 4 : chi**2 along Y
950 // if( tp->GetChi2Y() > 0 && (tp->GetChi2Y() < ti->GetChi2Y() || ti->GetChi2Y() <=0) ) tp_score++ ;
951 // if( ti->GetChi2Y() > 0 && (ti->GetChi2Y() < tp->GetChi2Y() || tp->GetChi2Y() <=0) ) ti_score++ ;
952 // CHECK 5 : total chi**2
953 // if( tp->chi2 > 0 && (tp->chi2 < ti->chi2 || ti->chi2 <=0) ) tp_score++ ;
954 // if( ti->chi2 > 0 && (ti->chi2 < tp->chi2 || tp->chi2 <=0) ) ti_score++ ;
955
956 // cout << "TRK "<<tp_score<<ti_score<<endl;
957 };
958
959
960
961 // *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
962 // the winner is....
963 // *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
964 if (tp_score > ti_score) {
965
966
967 }else if (tp_score < ti_score) {
968
969
970 ts = ti;//its image!!
971 cs = ci;
972 os = oi;
973
974 ti = tp;//its image!!
975 ci = cp;
976 oi = op;
977
978 tp = ts;//its image!!
979 cp = cs;
980 op = os;
981
982
983 }else {
984
985 // cout << "Warning - track image ambiguity not solved" << endl;
986
987 };
988
989 }else{
990 // ts = tp;
991 // cs = cp;
992 // os = op;
993 };
994
995 // cout <<" SortTracks() "<<i<<" -- "<<ts<<endl;
996 // sorted_tracks->Add(ts);//save the track in the sorted array
997 // sorted_tracks.Add(ts);//save the track in the sorted array
998 // sorted_tracks.Add(tp);//save the track in the sorted array
999 // cout << "SortTracks:: sorted_tracks->Add(it) "<<i<<" "<<ts<<endl;
1000 // cout<<"o "<<tp<<endl;
1001 // cout<<"o "<<cp<<endl;
1002 // cout<<"o "<<op<<endl;
1003 new(ttsorted[i]) PamTrack(tp,cp,op);
1004 new(ttimage[i]) PamTrack(ti,ci,oi);
1005 };
1006
1007 if( tsorted->GetEntries() != trk2_obj->GetNTracks() ){
1008 cout << "void PamLevel2::SortTracks(): tsorted->GetEntries() "<<tsorted->GetEntries()<<" != trk2_obj->GetNTracks() = "<<trk2_obj->GetNTracks() <<endl;
1009 tsorted->Delete(); tsorted=0;
1010 timage->Delete(); timage=0;
1011 }
1012
1013 //Restore Object count
1014 //To save space in the table keeping track of all referenced objects
1015 //We reset the object count to what it was at the beginning of the event.
1016 TProcessID::SetObjectCount(ObjectNumber);
1017
1018 };
1019 //--------------------------------------
1020 //
1021 //
1022 //--------------------------------------
1023 /**
1024 * This method overrides TrkLevel2::GetTracks(), where sorting is done by decreasing number of fit points and increasing chi^2.
1025 * PamLevel2::GetTracks() keeps the same track order given by TrkLevel2::GetTracks(), but checks image selection by using calorimeter and ToF tracking information.
1026 */
1027 // TRefArray *PamLevel2::GetTracks(){
1028
1029 // // *-*-*-*-*-*-*-*-*-*-*-*-*
1030 // SortTracks("+CAL+TOF");
1031 // // *-*-*-*-*-*-*-*-*-*-*-*-*
1032
1033 // // return sorted_tracks;
1034 // return &sorted_tracks;
1035
1036 // };
1037 TClonesArray *PamLevel2::GetTracks(){
1038
1039 // *-*-*-*-*-*-*-*-*-*-*-*-*
1040 SortTracks();
1041 // *-*-*-*-*-*-*-*-*-*-*-*-*
1042
1043 return tsorted;
1044
1045 };
1046 //--------------------------------------
1047 //
1048 //
1049 //--------------------------------------
1050 /**
1051 * Retrieves the it-th Pamela "physical" track.
1052 * It override TrkLevel2::GetTrack(int it).
1053 * @param it Track number, ranging from 0 to GetNTracks().
1054 */
1055 PamTrack *PamLevel2::GetTrack(int it){
1056
1057 // if(!trk2_obj) return 0;
1058
1059 // // *-*-*-*+-*-*-*-*-*-*-*-*-*
1060 // SortTracks("+CAL+TOF");
1061 // // *-*-*-*-*-*-*-*-*-*-*-*-*
1062 // // if(!sorted_tracks)return 0;
1063 // if(sorted_tracks.IsEmpty())return 0;
1064
1065 // PamTrack *track = 0;
1066
1067 // if( it >=0 && it < trk2_obj->TrkLevel2::GetNTracks() && it<sorted_tracks.GetEntriesFast() ){
1068 // TrkTrack *t = (TrkTrack*)sorted_tracks.At(it);
1069 // track = GetPamTrackAlong(t);
1070 // }else{
1071 // cout << "PamLevel2::GetTrack(int) : tracker track SeqNo "<< it <<" does not exist (GetNTracks() = "<<trk2_obj->TrkLevel2::GetNTracks()<<")"<<endl;
1072 // };
1073
1074 // return track;
1075
1076 // cout << "PamLevel2::GetTrack(int it) "<<endl;
1077 // *-*-*-*-*-*-*-*-*-*-*-*-*
1078 SortTracks();
1079 // *-*-*-*-*-*-*-*-*-*-*-*-*
1080 if(!tsorted)return 0;
1081 if(!tsorted->GetEntries())return 0;
1082
1083 PamTrack *track = 0;
1084
1085 if( it >=0 && it < trk2_obj->TrkLevel2::GetNTracks()){
1086 // TrkTrack *t = (TrkTrack*)sorted_tracks.At(it);
1087 // track = GetPamTrackAlong(t);
1088 TClonesArray &t = *(tsorted);
1089 track = (PamTrack*)t[it];
1090
1091 }else{
1092 cout << "PamLevel2::GetTrack(int) : tracker track SeqNo "<< it <<" does not exist (GetNTracks() = "<<trk2_obj->TrkLevel2::GetNTracks()<<")"<<endl;
1093 };
1094
1095 return track;
1096
1097 };
1098
1099 //--------------------------------------
1100 //
1101 //
1102 //--------------------------------------
1103 /**
1104 * Retrieves (if present) the image of the it-th Pamela "physical" track, sorted by the method PamLevel2::SortTracks().
1105 * @param it Track number, ranging from 0 to GetNTracks().
1106 */
1107 PamTrack *PamLevel2::GetTrackImage(int it){
1108
1109 // if(!trk2_obj) return 0;
1110
1111 // // *-*-*-*-*-*-*-*-*-*-*-*-*
1112 // SortTracks("+CAL+TOF");
1113 // // *-*-*-*-*-*-*-*-*-*-*-*-*
1114 // // if(!sorted_tracks)return 0;
1115 // if(sorted_tracks.IsEmpty())return 0;
1116
1117 // PamTrack *image = 0;
1118
1119 // if( it >=0 && it < trk2_obj->TrkLevel2::GetNTracks() ){
1120 // TrkTrack *temp = (TrkTrack*)sorted_tracks.At(it);
1121 // if( temp->HasImage() ){
1122 // TrkTrack *t = trk2_obj->TrkLevel2::GetStoredTrack(temp->GetImageSeqNo());
1123 // image = GetPamTrackAlong(t);
1124 // }else{
1125 // cout <<"PamLevel2::GetTrackImage(int) : Track SeqNo "<<it<<" does not have image"<<endl;
1126 // };
1127 // }else{
1128 // cout << "PamLevel2::GetTrackImage(int) : Tracker track SeqNo "<< it <<" does not exist (GetNTracks() = "<<trk2_obj->TrkLevel2::GetNTracks()<<")"<<endl;
1129 // };
1130
1131 // return image;
1132
1133
1134 // *-*-*-*-*-*-*-*-*-*-*-*-*
1135 SortTracks();
1136 // *-*-*-*-*-*-*-*-*-*-*-*-*
1137 if(!timage)return 0;
1138 if(!timage->GetEntries())return 0;
1139
1140 PamTrack *image = 0;
1141
1142 if( it >=0 && it < trk2_obj->TrkLevel2::GetNTracks() ){
1143 TClonesArray &t = *(tsorted);
1144 PamTrack *temp = (PamTrack*)t[it];
1145 if( temp->GetTrkTrack()->HasImage() ){
1146 TClonesArray &t = *(timage);
1147 image = (PamTrack*)t[it];
1148 }else{
1149 // cout <<"PamLevel2::GetTrackImage(int) : Track SeqNo "<<it<<" does not have image"<<endl;
1150 };
1151 }else{
1152 cout << "PamLevel2::GetTrackImage(int) : Tracker track SeqNo "<< it <<" does not exist (GetNTracks() = "<<trk2_obj->TrkLevel2::GetNTracks()<<")"<<endl;
1153 };
1154
1155 return image;
1156 }
1157
1158 //--------------------------------------
1159 //
1160 //
1161 //--------------------------------------
1162 /**
1163 * Get the Pamela detector trees in a single file and make them friends.
1164 * @param f TFile pointer
1165 * @param detlist String to select trees to be included
1166 * @return Pointer to a TTree
1167 */
1168 TTree *PamLevel2::GetPamTree(TFile *f, TString detlist ){
1169
1170 // if( !detlist.IsNull() ) SetWhichTrees(detlist);
1171 // else GetWhichTrees(f);
1172 if ( pam_tree ){
1173 printf("WARNING: TTree *PamLevel2::GetPamTree(TFile *fl, TString detlist) -- pam_tree already exists!\n ");
1174 return pam_tree;
1175 };
1176 //
1177
1178 cout << "TTree *PamLevel2::GetPamTree(TFile *f, TString detlist ) -- obsolte "<<endl;
1179
1180 SetWhichTrees(detlist);
1181
1182 TTree *Trout =0;
1183
1184 TString fname = f->GetName();
1185 if(!CheckLevel2File(fname))return NULL;
1186
1187 // UInt_t *found=0;
1188
1189 cout<< "GetPamTree(TFile*,TString): detector list --> ";
1190 if(TRK1)cout<<"TRK1 ";
1191 if(TRK2)cout<<"TRK2 ";
1192 if(TRKh)cout<<"TRKH ";
1193 if(CAL1)cout<<"CAL1 ";
1194 if(CAL2)cout<<"CAL2 ";
1195 if(TOF)cout<<"TOF ";
1196 if(TRG)cout<<"TRG ";
1197 if(AC)cout<<"AC ";
1198 if(ND)cout<<"ND ";
1199 if(S4)cout<<"S4 ";
1200 if(ORB)cout<<"ORB ";
1201 cout << endl;
1202 if(SELLI)cout<<">>> Found selection-list <<<"<<endl;
1203
1204 f->cd();
1205
1206 // Tracker
1207 TTree *T = (TTree*)f->Get("Tracker");
1208 if(T && (TRK2||TRK1||TRKh)) {
1209 if(TRK2)T->SetBranchAddress("TrkLevel2", GetPointerTo("TrkLevel2"));
1210 // else T->SetBranchStatus("TrkLevel2",0,found);
1211 if(TRK2)cout << "Tracker : set branch address TrkLevel2"<<endl;
1212 if(TRK1)T->SetBranchAddress("TrkLevel1", GetPointerTo("TrkLevel1"));
1213 // else T->SetBranchStatus("TrkLevel1",0,found);
1214 if(TRK1)cout << "Tracker : set branch address TrkLevel1"<<endl;
1215 if(TRKh)T->SetBranchAddress("TrkHough", GetPointerTo("TrkHough"));
1216 // else T->SetBranchStatus("TrkHough",0,found);
1217 if(TRKh)cout << "Tracker : set branch address TrkHough"<<endl;
1218 if(!Trout)Trout=T;
1219 else Trout->AddFriend(T);
1220 }else{
1221 cout << "Tracker : missing tree"<<endl;
1222 };
1223 // Calorimeter
1224 TTree *C = (TTree*)f->Get("Calorimeter");
1225 if(C && (CAL2||CAL1)) {
1226 if(CAL2)C->SetBranchAddress("CaloLevel2", GetPointerTo("CaloLevel2"));
1227 // else C->SetBranchStatus("CaloLevel2",0,found);
1228 if(CAL2)cout << "Calorimeter : set branch address CaloLevel2"<<endl;
1229 if(CAL1)C->SetBranchAddress("CaloLevel1", GetPointerTo("CaloLevel1"));
1230 // else C->SetBranchStatus("CaloLevel1",0,found);
1231 if(CAL1)cout << "Calorimeter : set branch address CaloLevel1"<<endl;
1232 if(!Trout)Trout=C;
1233 else Trout->AddFriend(C);
1234 }else{
1235 cout << "Calorimeter : missing tree"<<endl;
1236 };
1237
1238 // ToF
1239 TTree *O = (TTree*)f->Get("ToF");
1240 if(O && TOF) {
1241 O->SetBranchAddress("ToFLevel2", GetPointerTo("ToFLevel2"));
1242 cout << "ToF : set branch address ToFLevel2"<<endl;
1243 if(!Trout)Trout=O;
1244 else Trout->AddFriend(O);
1245 }else{
1246 cout << "ToF : missing tree"<<endl;
1247 };
1248 // Trigger
1249 TTree *R = (TTree*)f->Get("Trigger");
1250 if(R && TRG) {
1251 R->SetBranchAddress("TrigLevel2", GetPointerTo("TrigLevel2"));
1252 cout << "Trigger : set branch address TrigLevel2"<<endl;
1253 if(!Trout)Trout=R;
1254 else Trout->AddFriend(R);
1255 }else{
1256 cout << "Trigger : missing tree"<<endl;
1257 };
1258 // S4
1259 TTree *S = (TTree*)f->Get("S4");
1260 if(S && S4) {
1261 S->SetBranchAddress("S4Level2", GetPointerTo("S4Level2"));
1262 cout << "S4 : set branch address S4Level2"<<endl;
1263 if(!Trout)Trout=S;
1264 else Trout->AddFriend(S);
1265 }else{
1266 cout << "S4 : missing tree"<<endl;
1267 };
1268 // Neutron Detector
1269 TTree *N = (TTree*)f->Get("NeutronD");
1270 if(N && ND) {
1271 N->SetBranchAddress("NDLevel2", GetPointerTo("NDLevel2"));
1272 cout << "NeutronD : set branch address NDLevel2"<<endl;
1273 if(!Trout)Trout=N;
1274 else Trout->AddFriend(N);
1275 }else{
1276 cout << "NeutronD : missing tree"<<endl;
1277 };
1278 // Anticounters
1279 TTree *A = (TTree*)f->Get("Anticounter");
1280 if(A && AC) {
1281 A->SetBranchAddress("AcLevel2", GetPointerTo("AcLevel2"));
1282 cout << "Anticounter : set branch address AcLevel2"<<endl;
1283 if(!Trout)Trout=A;
1284 else Trout->AddFriend(A);
1285 }else{
1286 cout << "Anticounter : missing tree"<<endl;
1287 };
1288 // Orbital Info
1289 TTree *B = (TTree*)f->Get("OrbitalInfo");
1290 if(B && ORB) {
1291 B->SetBranchAddress("OrbitalInfo", GetPointerTo("OrbitalInfo"));
1292 cout << "OrbitalInfo : set branch address OrbitalInfo"<<endl;
1293 if(!Trout)Trout=B;
1294 else Trout->AddFriend(B);
1295 }else{
1296 cout << "OrbitalInfo : missing tree"<<endl;
1297 };
1298
1299 TTree *L = (TTree*)f->Get("SelectionList");
1300 if(L && SELLI==1) {
1301 // L->SetBranchAddress("RunEntry",&irun);
1302 // cout << "SelectionList: set branch address RunEntry"<<endl;
1303 // L->SetBranchAddress("EventEntry",&irunentry);
1304 // cout << "SelectionList: set branch address EventEntry"<<endl;
1305 // if(!Trout)Trout=O;
1306 // else Trout->AddFriend("SelectionList");
1307 cout << " TTree *PamLevel2::GetPamTree(TFile, TString) >>> SelectionList not implemente!!!!"<<endl;
1308 sel_tree = 0;
1309 }else{
1310 cout << "SelectionList : missing tree"<<endl;
1311 };
1312
1313 cout<<endl<<" Number of entries: "<<Trout->GetEntries()<<endl<<endl;
1314
1315 pam_tree = (TChain*)Trout;
1316
1317 return Trout;
1318
1319 }
1320 //--------------------------------------
1321 //
1322 //
1323 //--------------------------------------
1324 /**
1325 * Get list of Level2 files.
1326 * @param ddir Level2 data directory.
1327 * @param flisttxt Name of txt file containing file list.
1328 * @return Pointer to a TList of TSystemFiles
1329 * If no input file list is given , all the Level2 files inside the directory are processed.
1330 */
1331 TList* PamLevel2::GetListOfLevel2Files(TString ddir, TString flisttxt = ""){
1332
1333 TString wdir = gSystem->WorkingDirectory();
1334
1335 // if(ddir=="")ddir = wdir;
1336 // TSystemDirectory *datadir = new TSystemDirectory(gSystem->BaseName(ddir),ddir);
1337 if ( ddir != ""){
1338 cout << "Level2 data directory : "<< ddir << endl;
1339 } else {
1340 cout << "Level2 data directory not given as input: trying to evaluate from list or taking working directory " << endl;
1341 };
1342 TList *contents = new TList; // create output list
1343 contents->SetOwner();
1344
1345 // char *fullpath;
1346 // const char *fullpath;
1347
1348 // if no input file list is given:
1349 if ( flisttxt != "" ){
1350
1351 // if( !gSystem->IsFileInIncludePath(flisttxt,&fullpath) ){
1352 // if( !(fullpath = gSystem->FindFile("./",flisttxt)) ){
1353 // cout <<"File "<<flisttxt<<" not found"<<endl;
1354 // return 0;
1355 // }
1356 // flisttxt = fullpath;
1357 if ( !flisttxt.EndsWith(".root") ){
1358
1359 flisttxt = gSystem->ConcatFileName(gSystem->DirName(flisttxt),gSystem->BaseName(flisttxt));
1360
1361 if( !gSystem->ChangeDirectory(ddir) ){
1362 cout << "Cannot change directory : "<<ddir<<endl;
1363 return 0;
1364 }
1365
1366 cout <<"Input file list : " << flisttxt <<endl;
1367 ifstream in;
1368 in.open(flisttxt, ios::in); //open input file list
1369 if(!in.good()){
1370 cout <<" ERROR opening the file "<<endl;
1371 gSystem->ChangeDirectory(wdir); // back to the working directory
1372 return 0;
1373 }
1374 int line=0;
1375 while (1) {
1376 TString file;
1377 in >> file;
1378 if (!in.good()) break;
1379 line++;
1380 // cout <<"(1) " << file << endl;
1381 if(file.IsNull()){
1382 cout << "-- list interrupted at line "<<line <<endl;
1383 break;
1384 }
1385 if(file.Contains("#"))file = file(0,file.First("#"));
1386 // cout <<"(2) " << file << endl;
1387 // if( gSystem->IsFileInIncludePath(file,&fullpath) ){
1388 // if( (fullpath = gSystem->FindFile(ddir,file)) ){
1389 if( file.EndsWith(".root") ){
1390 TString filedir;
1391 if (ddir != ""){
1392 filedir = ddir; // take the input dir
1393 } else {
1394 gSystem->ChangeDirectory(wdir); // back to the working directory
1395 filedir = gSystem->DirName(file); // this will take the path if exist in the list otherwise it will return automatically the working dir
1396 };
1397 char *fullpath = gSystem->ConcatFileName(gSystem->DirName(filedir),gSystem->BaseName(file));
1398 contents->Add(new TSystemFile(fullpath,gSystem->DirName(fullpath)));// add file to the list
1399 delete fullpath;
1400 }
1401 // }else{
1402 // if(file.Data()!="")cout << "File: "<<file<<" ---> missing "<< endl;
1403 // };
1404 };
1405 in.close();
1406 } else {
1407 if(flisttxt.Contains("#"))flisttxt = flisttxt(0,flisttxt.First("#"));
1408 char *fullpath = gSystem->ConcatFileName(gSystem->DirName(flisttxt),gSystem->BaseName(flisttxt));
1409 contents->Add(new TSystemFile(fullpath,gSystem->DirName(fullpath)));// add file to the list
1410 delete fullpath;
1411 };
1412 }else{
1413
1414 cout << "No input file list given."<<endl;
1415 cout << "Check for existing root files."<<endl;
1416 // cout << "Warking directory: "<< gSystem->WorkingDirectory()<< endl;
1417 if ( ddir == "" ){
1418 ddir = wdir;
1419 cout << "Level2 data directory : "<< ddir << endl;
1420 };
1421
1422 TSystemDirectory *datadir = new TSystemDirectory(gSystem->BaseName(ddir),ddir);
1423 TList *temp = datadir->GetListOfFiles();
1424 // temp->Print();
1425 // cout << "*************************************" << endl;
1426
1427 TIter next(temp);
1428 TSystemFile *questo = 0;
1429
1430
1431 while ( (questo = (TSystemFile*) next()) ) {
1432 TString name = questo-> GetName();
1433 if( name.EndsWith(".root") ){
1434 // const char *fullpath = gSystem->FindFile(ddir,name);
1435 // char *fullpath;
1436 // gSystem->IsFileInIncludePath(name,&fullpath);
1437 char *fullpath = gSystem->ConcatFileName(gSystem->DirName(ddir),gSystem->BaseName(name));
1438 contents->Add(new TSystemFile(fullpath,gSystem->DirName(fullpath)));
1439 delete fullpath;
1440 };
1441 }
1442 delete temp;
1443 delete datadir;
1444
1445 };
1446 gSystem->ChangeDirectory(wdir); // back to the working directory
1447 // cout << endl << "Selected files:" << endl;
1448 // contents->Print();
1449 cout << contents->GetEntries()<<" files \n";
1450 // cout << endl;
1451 // cout << "Working directory: "<< gSystem->WorkingDirectory()<< endl;
1452 return contents;
1453 };
1454 //--------------------------------------
1455 //
1456 //
1457 //--------------------------------------
1458 /**
1459 * Get the Pamela detector chains from a list of files and make them friends.
1460 * @param fl Pointer to a TList of TSystemFiles
1461 * @param detlist String to select trees to be included
1462 * @return Pointer to a TChain
1463 */
1464 TChain *PamLevel2::GetPamTree(TList *fl, TString detlist ){
1465 //
1466 //
1467 //
1468 if ( pam_tree ){
1469 printf("WARNING: TChain *PamLevel2::GetPamTree(TList *fl, TString detlist) -- pam_tree already exists!\n ");
1470 return pam_tree;
1471 };
1472 //
1473
1474 TChain *Trout =0;
1475
1476 // if( !detlist.IsNull() )SetWhichTrees(detlist);
1477 SetWhichTrees(detlist);
1478
1479 cout<< "GetPamTree(TList*,TString): input detector list --> ";
1480 if(TRK1)cout<<"TRK1 ";
1481 if(TRK2)cout<<"TRK2 ";
1482 if(TRKh)cout<<"TRKH ";
1483 if(CAL1)cout<<"CAL1 ";
1484 if(CAL2)cout<<"CAL2 ";
1485 if(TOF)cout<<"TOF ";
1486 if(TRG)cout<<"TRG ";
1487 if(AC)cout<<"AC ";
1488 if(ND)cout<<"ND ";
1489 if(S4)cout<<"S4 ";
1490 if(ORB)cout<<"ORB ";
1491 cout << endl;
1492
1493 TChain *T = 0;
1494 TChain *C = 0;
1495 TChain *O = 0;
1496 TChain *R = 0;
1497 TChain *S = 0;
1498 TChain *N = 0;
1499 TChain *A = 0;
1500 TChain *B = 0;
1501
1502 TChain *L = 0;
1503
1504 if(TRK2||TRK1||TRKh) T = new TChain("Tracker");
1505 if(CAL2||CAL1) C = new TChain("Calorimeter");
1506 if(TOF) O = new TChain("ToF");
1507 if(TRG) R = new TChain("Trigger");
1508 if(S4) S = new TChain("S4");
1509 if(ND) N = new TChain("NeutronD");
1510 if(AC) A = new TChain("Anticounter");
1511 if(ORB) B = new TChain("OrbitalInfo");
1512 L = new TChain("SelectionList");
1513
1514 // loop over files and create chains
1515 TIter next(fl);
1516 TSystemFile *questo = 0;
1517 while ( (questo = (TSystemFile*) next()) ) {
1518 TString name = questo->GetName();
1519 cout << "File: "<< name << endl;
1520 if( CheckLevel2File(name) ){
1521 if(TRK2||TRK1||TRKh) T->Add(name);
1522 if(CAL1||CAL2) C->Add(name);
1523 if(TOF) O->Add(name);
1524 if(TRG) R->Add(name);
1525 if(S4) S->Add(name);
1526 if(ND) N->Add(name);
1527 if(AC) A->Add(name);
1528 if(ORB) B->Add(name);
1529 if(SELLI==1) L->Add(name);
1530 };
1531 };
1532
1533 cout << "done chain \n";
1534 cout << "+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+" <<endl;
1535
1536 // UInt_t *found=0;
1537 // Tracker
1538 if(T && (TRK2||TRK1||TRKh)) {
1539 if(TRK2)T->SetBranchAddress("TrkLevel2", GetPointerTo("TrkLevel2"));
1540 // else T->SetBranchStatus("TrkLevel2",0,found);
1541 if(TRK2)cout << "Tracker : set branch address TrkLevel2"<<endl;
1542 if(TRK1)T->SetBranchAddress("TrkLevel1", GetPointerTo("TrkLevel1"));
1543 // else T->SetBranchStatus("TrkLevel1",0,found);
1544 if(TRK1)cout << "Tracker : set branch address TrkLevel1"<<endl;
1545 if(TRKh)T->SetBranchAddress("TrkHough", GetPointerTo("TrkHough"));
1546 // else T->SetBranchStatus("TrkHough",0,found);
1547 if(TRKh)cout << "Tracker : set branch address TrkHough"<<endl;
1548 if(!Trout)Trout=T;
1549 else Trout->AddFriend("Tracker");
1550 }else{
1551 cout << "Tracker : missing tree"<<endl;
1552 };
1553 // Calorimeter
1554 if(C && (CAL2||CAL1)) {
1555 if(CAL2)C->SetBranchAddress("CaloLevel2", GetPointerTo("CaloLevel2"));
1556 // else C->SetBranchStatus("CaloLevel2",0,found);
1557 if(CAL2)cout << "Calorimeter : set branch address CaloLevel2"<<endl;
1558 if(CAL1)C->SetBranchAddress("CaloLevel1", GetPointerTo("CaloLevel1"));
1559 // else C->SetBranchStatus("CaloLevel1",0,found);
1560 if(CAL1)cout << "Calorimeter : set branch address CaloLevel1"<<endl;
1561 if(!Trout)Trout=C;
1562 else Trout->AddFriend("Calorimeter");
1563 }else{
1564 cout << "Calorimeter : missing tree"<<endl;
1565 };
1566 // ToF
1567 if(O && TOF) {
1568 O->SetBranchAddress("ToFLevel2", GetPointerTo("ToFLevel2"));
1569 cout << "ToF : set branch address ToFLevel2"<<endl;
1570 if(!Trout)Trout=O;
1571 else Trout->AddFriend("ToF");
1572 }else{
1573 cout << "ToF : missing tree"<<endl;
1574 };
1575 // Trigger
1576 if(R && TRG) {
1577 R->SetBranchAddress("TrigLevel2", GetPointerTo("TrigLevel2"));
1578 cout << "Trigger : set branch address TrigLevel2"<<endl;
1579 if(!Trout)Trout=O;
1580 else Trout->AddFriend("Trigger");
1581 }else{
1582 cout << "Trigger : missing tree"<<endl;
1583 };
1584 // S4
1585 if(S && S4) {
1586 S->SetBranchAddress("S4Level2", GetPointerTo("S4Level2"));
1587 cout << "S4 : set branch address S4Level2"<<endl;
1588 if(!Trout)Trout=O;
1589 else Trout->AddFriend("S4");
1590 }else{
1591 cout << "S4 : missing tree"<<endl;
1592 };
1593 // Neutron Detector
1594 if(N && ND) {
1595 N->SetBranchAddress("NDLevel2", GetPointerTo("NDLevel2"));
1596 cout << "NeutronD : set branch address NDLevel2"<<endl;
1597 if(!Trout)Trout=O;
1598 else Trout->AddFriend("NeutronD");
1599 }else{
1600 cout << "NeutronD : missing tree"<<endl;
1601 };
1602 // Anticounters
1603 if(A && AC) {
1604 A->SetBranchAddress("AcLevel2", GetPointerTo("AcLevel2"));
1605 cout << "Anticounter : set branch address AcLevel2"<<endl;
1606 if(!Trout)Trout=O;
1607 else Trout->AddFriend("Anticounter");
1608 }else{
1609 cout << "Anticounter : missing tree"<<endl;
1610 };
1611 // Orbital Info
1612 if(B && ORB) {
1613 B->SetBranchAddress("OrbitalInfo", GetPointerTo("OrbitalInfo"));
1614 cout << "OrbitalInfo : set branch address OrbitalInfo"<<endl;
1615 if(!Trout)Trout=O;
1616 else Trout->AddFriend("OrbitalInfo");
1617 }else{
1618 cout << "OrbitalInfo : missing tree"<<endl;
1619 };
1620 cout << "+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+" <<endl;
1621
1622 // Selection List
1623 if(L && SELLI==1) {
1624 cout<<">>> Found selection-list <<<"<<endl;
1625 L->SetBranchAddress("RunEntry",&irun);
1626 cout << "SelectionList: set branch address RunEntry"<<endl;
1627 L->SetBranchAddress("EventEntry",&irunentry);
1628 cout << "SelectionList: set branch address EventEntry"<<endl;
1629 sel_tree = L;
1630 // if(!Trout)Trout=O;
1631 // else Trout->AddFriend("SelectionList");
1632 cout << "----------------------------------------------------" <<endl;
1633 }else{
1634 // cout << "SelectionList : missing tree"<<endl;
1635 if(L)L->Delete();
1636 };
1637
1638 // cout<<endl<<" Number of entries: "<<Trout->GetEntries()<<endl<<endl;
1639
1640 pam_tree = Trout;
1641
1642 return Trout;
1643 }
1644
1645
1646
1647 //--------------------------------------
1648 //
1649 //
1650 //--------------------------------------
1651 /**
1652 * Set branch addresses for Pamela friend trees
1653 */
1654 void PamLevel2::SetBranchAddress(TTree *t){
1655
1656 TRK2 = TRK2 & t->GetBranchStatus("TrkLevel2");
1657 TRK1 = TRK1 & t->GetBranchStatus("TrkLevel1");
1658 TRKh = TRKh & t->GetBranchStatus("TrkHough");
1659 CAL2 = CAL2 & t->GetBranchStatus("CaloLevel2");
1660 CAL1 = CAL1 & t->GetBranchStatus("CaloLevel1");
1661 TOF = TOF & t->GetBranchStatus("ToFLevel2");
1662 TRG = TRG & t->GetBranchStatus("TrigLevel2");
1663 S4 = S4 & t->GetBranchStatus("S4Level2");
1664 ND = ND & t->GetBranchStatus("NDLevel2");
1665 AC = AC & t->GetBranchStatus("AcLevel2");
1666 ORB = ORB & t->GetBranchStatus("OrbitalInfo");
1667
1668
1669 // Tracker
1670 if(TRK1) {
1671 t->SetBranchAddress("TrkLevel1", GetPointerTo("TrkLevel2"));
1672 cout << "Tracker : set branch address TrkLevel1"<<endl;
1673 };
1674 if(TRK2) {
1675 t->SetBranchAddress("TrkLevel2", GetPointerTo("TrkLevel1"));
1676 cout << "Tracker : set branch address TrkLevel2"<<endl;
1677 };
1678 if(TRKh) {
1679 t->SetBranchAddress("TrkHough", GetPointerTo("TrkHough"));
1680 cout << "Tracker : set branch address TrkHough"<<endl;
1681 };
1682
1683 // Calorimeter
1684 if(CAL1) {
1685 t->SetBranchAddress("CaloLevel1", GetPointerTo("CaloLevel1"));
1686 cout << "Calorimeter : set branch address CaloLevel1"<<endl;
1687 };
1688 if(CAL2) {
1689 t->SetBranchAddress("CaloLevel2", GetPointerTo("CaloLevel2"));
1690 cout << "Calorimeter : set branch address CaloLevel2"<<endl;
1691 };
1692
1693 // ToF
1694 if(TOF) {
1695 t->SetBranchAddress("ToFLevel2", GetPointerTo("ToFLevel2"));
1696 cout << "ToF : set branch address ToFLevel2"<<endl;
1697 };
1698 // Trigger
1699 if(TRG) {
1700 t->SetBranchAddress("TrigLevel2", GetPointerTo("TrigLevel2"));
1701 cout << "Trigger : set branch address TrigLevel2"<<endl;
1702 };
1703 // S4
1704 if(S4) {
1705 t->SetBranchAddress("S4Level2", GetPointerTo("S4Level2"));
1706 cout << "S4 : set branch address S4Level2"<<endl;
1707 };
1708 // Neutron Detector
1709 if(ND) {
1710 t->SetBranchAddress("NDLevel2", GetPointerTo("NDLevel2"));
1711 cout << "NeutronD : set branch address NDLevel2"<<endl;
1712 };
1713 // Anticounters
1714 if(AC) {
1715 t->SetBranchAddress("AcLevel2", GetPointerTo("AcLevel2"));
1716 cout << "Anticounter : set branch address AcLevel2"<<endl;
1717 };
1718 // OrbitalInfo
1719 if(ORB) {
1720 t->SetBranchAddress("OrbitalInfo", GetPointerTo("OrbitalInfo"));
1721 cout << "OrbitalInfo : set branch address OrbitalInfo"<<endl;
1722 };
1723 // SelectionList
1724 if(SELLI==1) {
1725 t->SetBranchAddress("RunEntry", &irun);
1726 cout << "SelectionList: set branch address RunEntry"<<endl;
1727 t->SetBranchAddress("EventEntry", &irunentry);
1728 cout << "SelectionList: set branch address EventEntry"<<endl;
1729 };
1730
1731 }
1732 /**
1733 * Set branch addresses for Pamela friend trees
1734 */
1735 void PamLevel2::SetBranchAddress(TChain *t){
1736
1737 TRK2 = TRK2 & t->GetBranchStatus("TrkLevel2");
1738 TRK1 = TRK1 & t->GetBranchStatus("TrkLevel1");
1739 TRKh = TRKh & t->GetBranchStatus("TrkHough");
1740 CAL1 = CAL1 & t->GetBranchStatus("CaloLevel1");
1741 CAL2 = CAL2 & t->GetBranchStatus("CaloLevel2");
1742 TOF = TOF & t->GetBranchStatus("ToFLevel2");
1743 TRG = TRG & t->GetBranchStatus("TrigLevel2");
1744 S4 = S4 & t->GetBranchStatus("S4Level2");
1745 ND = ND & t->GetBranchStatus("NDLevel2");
1746 AC = AC & t->GetBranchStatus("AcLevel2");
1747 ORB = ORB & t->GetBranchStatus("OrbitalInfo");
1748
1749 // Tracker
1750 if(TRK2) {
1751 t->SetBranchAddress("TrkLevel2", GetPointerTo("TrkLevel2"));
1752 cout << "Tracker : set branch address TrkLevel2"<<endl;
1753 };
1754 if(TRK1) {
1755 t->SetBranchAddress("TrkLevel1", GetPointerTo("TrkLevel1"));
1756 cout << "Tracker : set branch address TrkLevel1"<<endl;
1757 };
1758 if(TRKh) {
1759 t->SetBranchAddress("TrkHough", GetPointerTo("TrkHough"));
1760 cout << "Tracker : set branch address TrkHough"<<endl;
1761 };
1762
1763 // Calorimeter
1764 if(CAL2) {
1765 t->SetBranchAddress("CaloLevel2", GetPointerTo("CaloLevel2"));
1766 cout << "Calorimeter : set branch address CaloLevel2"<<endl;
1767 };
1768 if(CAL1) {
1769 t->SetBranchAddress("CaloLevel1", GetPointerTo("CaloLevel1"));
1770 cout << "Calorimeter : set branch address CaloLevel1"<<endl;
1771 };
1772
1773 // ToF
1774 if(TOF) {
1775 t->SetBranchAddress("ToFLevel2", GetPointerTo("ToFLevel2"));
1776 cout << "ToF : set branch address ToFLevel2"<<endl;
1777 };
1778 // Trigger
1779 if(TRG) {
1780 t->SetBranchAddress("TrigLevel2", GetPointerTo("TrigLevel2"));
1781 cout << "Trigger : set branch address TrigLevel2"<<endl;
1782 };
1783 // S4
1784 if(S4) {
1785 t->SetBranchAddress("S4Level2", GetPointerTo("S4Level2"));
1786 cout << "S4 : set branch address S4Level2"<<endl;
1787 };
1788 // Neutron Detector
1789 if(ND) {
1790 t->SetBranchAddress("NDLevel2", GetPointerTo("NDLevel2"));
1791 cout << "NeutronD : set branch address NDLevel2"<<endl;
1792 };
1793 // Anticounters
1794 if(AC) {
1795 t->SetBranchAddress("AcLevel2", GetPointerTo("AcLevel2"));
1796 cout << "Anticounter : set branch address AcLevel2"<<endl;
1797 };
1798 // OrbitalInfo
1799 if(ORB) {
1800 t->SetBranchAddress("OrbitalInfo", GetPointerTo("OrbitalInfo"));
1801 cout << "OrbitalInfo : set branch address OrbitalInfo"<<endl;
1802 };
1803 // SelectionList
1804 if(SELLI==1) {
1805 t->SetBranchAddress("RunEntry", &irun);
1806 cout << "SelectionList: set branch address RunEntry"<<endl;
1807 t->SetBranchAddress("EventEntry", &irunentry);
1808 cout << "SelectionList: set branch address EventEntry"<<endl;
1809 };
1810
1811 }
1812
1813
1814 //--------------------------------------
1815 //
1816 //
1817 //--------------------------------------
1818 /**
1819 * Get the Run tree chain from a list of files.
1820 * @param fl Pointer to a TList of TSystemFiles
1821 * @return Pointer to a TChain
1822 */
1823 TChain *PamLevel2::GetRunTree(TList *fl){
1824 //
1825 //
1826 //
1827 if ( run_tree ){
1828 printf("WARNING: TChain *PamLevel2::GetRunTree(TList *fl) -- run_tree already exists!\n ");
1829 return run_tree;
1830 };
1831 //
1832 TChain *R = new TChain("Run");
1833
1834 // loop over files and create chains
1835 TIter next(fl);
1836 TSystemFile *questo = 0;
1837 while ( (questo = (TSystemFile*) next()) ) {
1838 TString name = questo->GetName();
1839 // cout << "File: "<< name << endl;
1840 if( CheckLevel2File(name) ){
1841 R->Add(name);
1842 };
1843 }
1844
1845 if(R->GetNtrees()){
1846
1847 R->SetBranchAddress("RunInfo", GetPointerTo("RunInfo"));
1848 cout << "Run : set branch address RunInfo"<<endl;
1849 R->SetBranchAddress("SoftInfo", GetPointerTo("SoftInfo")); // Emiliano
1850 cout << "Software : set branch address SoftInfo"<<endl; // Emiliano
1851 }else{
1852 delete R;
1853 R=0;
1854 }
1855
1856 run_tree = R;
1857
1858 return R;
1859
1860 }
1861 //--------------------------------------
1862 //
1863 //
1864 //--------------------------------------
1865 /**
1866 * Get the Run tree from a file.
1867 * @param f Pointer to a TFile
1868 * @return Pointer to a TTree
1869 */
1870 TTree *PamLevel2::GetRunTree(TFile *f){
1871 if ( run_tree ){
1872 printf("WARNING: TTree *PamLevel2::GetRunTree(TFile *f) -- run_tree already exists!\n ");
1873 return run_tree;
1874 };
1875
1876
1877 cout << "TTree *PamLevel2::GetRunTree(TFile *f) -- obsolte "<<endl;
1878
1879 TTree *T = (TTree*)f->Get("Run");
1880
1881 if(T){
1882 T->SetBranchAddress("RunInfo", GetPointerTo("RunInfo"));
1883 cout << "Run : set branch address RunInfo"<<endl;
1884 T->SetBranchAddress("SoftInfo", GetPointerTo("SoftInfo")); // Emiliano
1885 cout << "Software : set branch address SoftInfo"<<endl; // Emiliano
1886 }
1887
1888 run_tree = (TChain*)T;
1889
1890 return T;
1891
1892 }
1893 /**
1894 * Update the runinfo informations (to be used to have Run infos event by event basis)
1895 * @param run Pointer to the chain/tree which contains run infos
1896 * @return true if a new run has been read, false if it is still the same run
1897 */
1898 Bool_t PamLevel2::UpdateRunInfo(TChain *run, Long64_t iev){
1899 //
1900 // check if we have already called once GetEntry, if not call it
1901 //
1902 cout << "Bool_t PamLevel2::UpdateRunInfo(TChain *run, Long64_t iev) --- NON FUNZIONA BENISSIMO.... "<<endl;
1903 if(!run){
1904 cout << "Bool_t PamLevel2::UpdateRunInfo(TChain *run, ULong64_t iev) -- ERROR -- missing RunInfo tree "<<endl;
1905 return(false);
1906 }
1907 if ( run->GetEntries() <= 0 ) return(false);
1908 //
1909
1910 // Int_t oldrun = irun;
1911 Long64_t oldrun = irun;
1912
1913 // --------------------------------------
1914 // if it is a full file (not preselected)
1915 // --------------------------------------
1916 if(SELLI==0){
1917
1918 //
1919 // the absolute time is necessary to relate the event with the run
1920 //
1921 if( !GetOrbitalInfo() ){
1922 cout << "Bool_t PamLevel2::UpdateRunInfo(TChain *run, ULong64_t iev) -- ERROR -- missing OrbitalInfo "<<endl;
1923 return(false);
1924 }
1925
1926 //
1927 // the first time the routine is called, set run search from the beginning
1928 //
1929 if ( irun < 0LL ){
1930 irun = 0LL;
1931 run->GetEntry(irun);
1932 runfirstentry = 0LL;
1933 runlastentry += (Long64_t)(this->GetRunInfo()->NEVENTS);
1934 if ( (Long64_t)(this->GetRunInfo()->NEVENTS) > 0LL ) runlastentry -= 1LL;
1935 };
1936 //
1937 if ( irun == run->GetEntries()-1LL &&
1938 !(GetOrbitalInfo()->absTime >= GetRunInfo()->RUNHEADER_TIME &&
1939 GetOrbitalInfo()->absTime <= GetRunInfo()->RUNTRAILER_TIME)
1940 ){
1941 irun = -1LL;
1942 runfirstentry = 0LL;
1943 runlastentry = -1LL;
1944 };
1945 // modificato il controllo sull'aggiornamento del run, per evitare problemi
1946 // dovuti agli eventi annidati (NB! NEVENTS conta anche questi!!)
1947 //
1948 bool fromfirst = true;
1949 //
1950 while ( !(GetOrbitalInfo()->absTime >= GetRunInfo()->RUNHEADER_TIME && GetOrbitalInfo()->absTime <= GetRunInfo()->RUNTRAILER_TIME) && irun < run->GetEntries()-1LL ){
1951 // while ( iev > (runfirstentry+(ULong64_t)(this->GetRunInfo()->NEVENTS-1)) && irun < run->GetEntries() ){
1952 irun++;
1953 run->GetEntry(irun);
1954 runfirstentry = runlastentry;
1955 if ( (Long64_t)(this->GetRunInfo()->NEVENTS) > 0LL ) runfirstentry += 1LL;
1956 runlastentry += (Long64_t)(this->GetRunInfo()->NEVENTS);
1957 // cout << " ))))) UPDATE RUN INFO ((((( @iev "<<iev<<" run "<<GetRunInfo()->ID<<" irun "<<irun<<endl;
1958 // cout << "runfirstentry "<<runfirstentry<<endl;
1959 // printf(" iev %llu %u %llu \n",iev,this->GetRunInfo()->NEVENTS,(ULong64_t)(runfirstentry+(ULong64_t)(this->GetRunInfo()->NEVENTS)));
1960 // printf(" abstime %u trailertime %u \n",GetOrbitalInfo()->absTime,GetRunInfo()->RUNTRAILER_TIME);
1961 // printf(" IDRUN %u \n",GetRunInfo()->ID);
1962 //
1963 // prevshift = 0;
1964 //
1965 if ( irun == (Long64_t)(run->GetEntries()-1LL) && fromfirst && !(GetOrbitalInfo()->absTime >= GetRunInfo()->RUNHEADER_TIME && GetOrbitalInfo()->absTime <= GetRunInfo()->RUNTRAILER_TIME)){
1966 printf(" resetting irun (it should NOT happen!!!)\n");
1967 fromfirst = false;
1968 irun = 0;
1969 run->GetEntry(irun);
1970 runfirstentry = 0ULL;
1971 runlastentry += (Long64_t)(this->GetRunInfo()->NEVENTS);
1972 if ( (Long64_t)(this->GetRunInfo()->NEVENTS) > 0LL ) runlastentry -= 1LL;
1973 };
1974 //
1975 };
1976 //
1977 if (
1978 !(GetOrbitalInfo()->absTime >= GetRunInfo()->RUNHEADER_TIME &&
1979 GetOrbitalInfo()->absTime <= GetRunInfo()->RUNTRAILER_TIME)
1980 ) {
1981 printf(" Something very wrong here: cannot find RUN containing absolute time %u \n",GetOrbitalInfo()->absTime);
1982 return false;
1983 }
1984 //
1985 if ( irun == oldrun || irun >= run->GetEntries() ) return(false);
1986 //
1987 // printf(" iev %llu irun %i nevents %u 1st %llu last %llu \n",iev,irun,this->GetRunInfo()->NEVENTS,(ULong64_t)runfirstentry,(ULong64_t)runlastentry);
1988 //
1989 prevshift = 0;
1990 cout << " ))))) UPDATE RUN INFO ((((( @iev "<<iev<<" run "<<GetRunInfo()->ID<<" irun "<<irun<<endl;
1991 // cout << "runfirstentry "<<runfirstentry<<endl;
1992 return(true);
1993 };
1994 // ----------------------------------------------------
1995 // if it is a preselected file (there is SelectionList)
1996 // NBNB - the event tree MUST be read first
1997 // ----------------------------------------------------
1998 if(SELLI==1){
1999 sel_tree->GetEntry(iev);
2000 // cout << irun << " "<< irunentry << endl;
2001 if(irun != oldrun){
2002 run->GetEntry(irun);
2003 cout << " ))))) UPDATE RUN INFO ((((( @iev "<<iev<<" run "<<GetRunInfo()->ID<<" irun "<<irun<<endl;
2004 prevshift = 0;
2005 return true;
2006 }
2007 return false;
2008 }
2009
2010 return false;
2011 //
2012 };
2013
2014 Bool_t PamLevel2::UpdateRunInfo(Long64_t iev){
2015
2016 if(!run_tree){
2017 cout << " Bool_t PamLevel2::UpdateRunInfo(ULong64_t iev) -- ERROR -- run tree not loaded"<<endl;
2018 return false;
2019 }
2020 if ( run_tree->GetEntries() <= 0 ) {
2021 cout << " Bool_t PamLevel2::UpdateRunInfo(ULong64_t iev) -- ERROR -- run tree is empty"<<endl;
2022 return(false);
2023 }
2024
2025 Int_t oldrun = irun; // store current run index
2026
2027 // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
2028 // if it is a full file (not preselected)
2029 // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
2030 if(SELLI==0){
2031
2032 // ---------------------------------------------------------------
2033 // the absolute time is necessary to relate the event with the run
2034 // ---------------------------------------------------------------
2035 if( !GetOrbitalInfo() ){
2036 cout << "Bool_t PamLevel2::UpdateRunInfo(Long64_t iev) -- ERROR -- missing OrbitalInfo "<<endl;
2037 return(false);
2038 }
2039
2040 // -----------------------------------------------------------------------
2041 // the first time the routine is called, set run search from the beginning
2042 // -----------------------------------------------------------------------
2043 if ( irun < 0 ){
2044 irun = 0LL;
2045 irunentry = 0;
2046 prevshift = 0;
2047 run_tree->GetEntry(irun);
2048 };
2049 //
2050 bool fromfirst = true; // first loop over runs
2051
2052 // Bool_t hasfrag = false;
2053 // if( GetRunInfo()->ID_RUN_FRAG!=0 && GetRunInfo()->ID_RUN_FRAG != GetRunInfo()->ID)hasfrag=true;
2054 // ULong64_t fragid = GetRunInfo()->ID_RUN_FRAG;
2055
2056 // ------------------------------------------------------
2057 // loop over runs to find the one that contains the event
2058 // ------------------------------------------------------
2059 while (
2060 (
2061 (
2062 !(GetOrbitalInfo()->absTime >= GetRunInfo()->RUNHEADER_TIME && // check on absolute time (s)
2063 GetOrbitalInfo()->absTime <= GetRunInfo()->RUNTRAILER_TIME) &&
2064 !(GetOrbitalInfo()->OBT >= GetRunInfo()->RUNHEADER_OBT && // additional check on OBT (ms)
2065 GetOrbitalInfo()->OBT <= GetRunInfo()->RUNTRAILER_OBT)
2066 )
2067 || GetRunInfo()->NEVENTS==0
2068 || !(irunentry < GetRunInfo()->NEVENTS-1-prevshift)
2069 )
2070 && irun < run_tree->GetEntries() ){
2071
2072 irun++;
2073 // ------------------------------------
2074 // if the end of run tree is reached...
2075 // ------------------------------------
2076 if( irun == run_tree->GetEntries() ){
2077 if(!fromfirst){
2078 // -----------------------------------------------------
2079 // if it happened already once and the run was not found
2080 // ---> exit with error
2081 // -----------------------------------------------------
2082 cout << "Bool_t PamLevel2::UpdateRunInfo(Long64_t iev) -- ERROR -- event entry #"<<iev<<" does not belong to any run (should not happen)" <<endl;
2083 return false;
2084 }
2085 // -----------------------------------------
2086 // ...otherwise repeat search from beginning
2087 // -----------------------------------------
2088 cout << "Bool_t PamLevel2::UpdateRunInfo(Long64_t iev) -- WARNING -- reached end of run tree. searchin again from beginning " <<endl;
2089 fromfirst = false;
2090 irun = 0LL;
2091 runfirstentry = 0LL;
2092 }
2093 // -------------------------------------------------------------------
2094 // save the index of the first entry of the run, relative to pam_tree,
2095 // and read a new run
2096 // -------------------------------------------------------------------
2097 if(irun>0)runfirstentry += (GetRunInfo()->NEVENTS)-prevshift;
2098 irunentry = 0;
2099 prevshift = 0;
2100 run_tree->GetEntry(irun);
2101 if(GetRunInfo()->RUNHEADER_OBT>GetRunInfo()->RUNTRAILER_OBT ){
2102 cout << "Bool_t PamLevel2::UpdateRunInfo(Long64_t iev) -- WARNING -- irun "<<irun<<" has RUNHEADER_OBT>=RUNTRAILER_OBT " <<endl;
2103 cout << " (NB!! in this case some events are assigned to a wrong run)"<<endl;
2104 }
2105 // if(hasfrag && fragid != GetRunInfo()->ID){
2106 // cout << "... where is the next fragment ??"<<endl;
2107 // }
2108 };
2109
2110
2111 // --------------------------------------
2112 // if there was no need to update the run
2113 // ---> exit with FALSE
2114 // --------------------------------------
2115 if ( irun == oldrun ) return(false);
2116
2117 // --------------------------------------
2118 // ... otherwise
2119 // ---> exit with TRUE
2120 // --------------------------------------
2121 cout << endl << " ))))) UPDATE RUN INFO ((((( @iev "<<iev<<" run "<<GetRunInfo()->ID<<" irun "<<irun<<endl;
2122 // ----------------------------------------------------
2123 // then check if the run has a fragment
2124 // in this case we have to switch to the next fragment
2125 // when the end of the first fragment is reached
2126 // ----------------------------------------------------
2127 if(
2128 GetRunInfo()->ID_RUN_FRAG != 0 &&
2129 // GetRunInfo()->ID_RUN_FRAG != GetRunInfo()->ID &&
2130 true ){
2131 cout << "* fragment *"<<endl;
2132 }
2133
2134 return(true);
2135 };
2136 // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
2137 // if it is a preselected file (there is SelectionList)
2138 // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
2139 if(SELLI==1){
2140 sel_tree->GetEntry(iev);
2141 if(irun != oldrun){
2142 run_tree->GetEntry(irun);
2143 cout << endl << " ))))) UPDATE RUN INFO ((((( @iev "<<iev<<" run "<<GetRunInfo()->ID<<" irun "<<irun<<endl;
2144 prevshift = 0;
2145 return true;
2146 }
2147 return false;
2148 }
2149
2150 return false;
2151 //
2152 };
2153 /**
2154 * Update the runinfo informations (to be used to have Run infos event by event basis)
2155 * @param run Pointer to the chain/tree which contains run infos
2156 * @return true if a new run has been read, false if it is still the same run
2157 */
2158 Bool_t PamLevel2::UpdateRunInfo(TTree *run, Long64_t iev){
2159 return(UpdateRunInfo((TChain*)run,iev));
2160 };
2161
2162 //--------------------------------------
2163 //
2164 //
2165 //--------------------------------------
2166 /**
2167 * Set which trees shoul be analysed
2168 * @param detlist TString containing the sequence of trees required
2169 */
2170 void PamLevel2::SetWhichTrees(TString detlist){
2171
2172 if(detlist.IsNull() || detlist.Contains("+ALL", TString::kIgnoreCase)){
2173 CAL0 = false;
2174 CAL1 = true;
2175 CAL2 = true;
2176 TRK2 = true;
2177 TRK1 = false;
2178 TRKh = false;
2179 TRK0 = false;
2180 TRG = true;
2181 TOF = true;
2182 TOF0 = false;
2183 S4 = true;
2184 ND = true;
2185 AC = true;
2186 ORB = true;
2187 }else if( detlist.Contains("-ALL", TString::kIgnoreCase) ){
2188 CAL0 = false;
2189 CAL1 = false;
2190 CAL2 = false;
2191 TRK2 = false;
2192 TRK1 = false;
2193 TRKh = false;
2194 TRK0 = false;
2195 TRG = false;
2196 TOF = false;
2197 TOF0 = false;
2198 S4 = false;
2199 ND = false;
2200 AC = false;
2201 ORB = false;
2202 };
2203
2204 // -------------------------------------------------------------------------
2205 if( detlist.Contains("CAL1", TString::kIgnoreCase) ){
2206 if ( detlist.Contains("-CAL1", TString::kIgnoreCase) )CAL1=false;
2207 if ( detlist.Contains("+CAL1", TString::kIgnoreCase) )CAL1=true;
2208 };
2209
2210 if( detlist.Contains("CAL0", TString::kIgnoreCase) ){
2211 if ( detlist.Contains("-CAL0", TString::kIgnoreCase) )CAL0=false;
2212 if ( detlist.Contains("+CAL0", TString::kIgnoreCase) )CAL0=true;
2213 };
2214
2215 if( detlist.Contains("CAL2", TString::kIgnoreCase)){
2216 if ( detlist.Contains("-CAL2", TString::kIgnoreCase) )CAL2=false;
2217 if ( detlist.Contains("+CAL2", TString::kIgnoreCase) )CAL2=true;
2218 };
2219
2220 if( detlist.Contains("+CAL", TString::kIgnoreCase) && !CAL1 && !CAL2 )CAL2=true;
2221 if( detlist.Contains("-CAL", TString::kIgnoreCase) && CAL1 && CAL2 ){
2222 CAL2=false;
2223 CAL1=false;
2224 }
2225 // -------------------------------------------------------------------------
2226 if( detlist.Contains("TRK0", TString::kIgnoreCase) ){
2227 if ( detlist.Contains("-TRK0", TString::kIgnoreCase) )TRK0=false;
2228 if ( detlist.Contains("+TRK0", TString::kIgnoreCase) )TRK0=true;
2229 };
2230
2231 if( detlist.Contains("TRK1", TString::kIgnoreCase) ){
2232 if ( detlist.Contains("-TRK1", TString::kIgnoreCase) )TRK1=false;
2233 if ( detlist.Contains("+TRK1", TString::kIgnoreCase) )TRK1=true;
2234 };
2235
2236 if( detlist.Contains("TRK2", TString::kIgnoreCase)){
2237 if ( detlist.Contains("-TRK2", TString::kIgnoreCase) )TRK2=false;
2238 if ( detlist.Contains("+TRK2", TString::kIgnoreCase) )TRK2=true;
2239 };
2240
2241 if( detlist.Contains("TRKh", TString::kIgnoreCase)){
2242 if ( detlist.Contains("-TRKh", TString::kIgnoreCase) )TRKh=false;
2243 if ( detlist.Contains("+TRKh", TString::kIgnoreCase) )TRKh=true;
2244 };
2245
2246 if( detlist.Contains("+TRK", TString::kIgnoreCase) && !TRK1 && !TRK2 && !TRKh )TRK2=true;
2247 if( detlist.Contains("-TRK", TString::kIgnoreCase) && TRK1 && TRK2 && TRKh){
2248 TRK2=false;
2249 TRK1=false;
2250 TRKh=false;
2251 }
2252 // -------------------------------------------------------------------------
2253
2254 if( detlist.Contains("-TRG", TString::kIgnoreCase) )TRG = false;
2255 else if( detlist.Contains("+TRG", TString::kIgnoreCase) )TRG = true;
2256
2257 if( detlist.Contains("-TOF", TString::kIgnoreCase) )TOF = false;
2258 else if( detlist.Contains("+TOF", TString::kIgnoreCase) )TOF = true;
2259
2260 if( detlist.Contains("-TOF0", TString::kIgnoreCase) )TOF0 = false;
2261 else if( detlist.Contains("+TOF0", TString::kIgnoreCase) )TOF0 = true;
2262
2263 if( detlist.Contains("-S4", TString::kIgnoreCase) )S4 = false;
2264 else if( detlist.Contains("+S4", TString::kIgnoreCase) )S4 = true;
2265
2266 if( detlist.Contains("-ND", TString::kIgnoreCase) )ND = false;
2267 else if( detlist.Contains("+ND", TString::kIgnoreCase) )ND = true;
2268
2269 if( detlist.Contains("-AC", TString::kIgnoreCase) )AC = false;
2270 else if( detlist.Contains("+AC", TString::kIgnoreCase) )AC = true;
2271
2272 if( detlist.Contains("-ORB", TString::kIgnoreCase) )ORB = false;
2273 else if( detlist.Contains("+ORB", TString::kIgnoreCase) )ORB = true;
2274
2275 // cout<< "Set detector list --> ";
2276 // if(TRK1)cout<<"TRK1 ";
2277 // if(TRK2)cout<<"TRK2 ";
2278 // if(TRKh)cout<<"TRKH ";
2279 // if(CAL1)cout<<"CAL1 ";
2280 // if(CAL2)cout<<"CAL2 ";
2281 // if(TOF)cout<<"TOF ";
2282 // if(TRG)cout<<"TRG ";
2283 // if(AC)cout<<"AC ";
2284 // if(ND)cout<<"ND ";
2285 // if(S4)cout<<"S4 ";
2286 // if(ORB)cout<<"ORB ";
2287 // cout << endl;
2288
2289 };
2290
2291
2292 /**
2293 * Set tree/branch detector flags from the content of a tree
2294 */
2295 void PamLevel2::GetWhichTrees(TFile* f){
2296
2297
2298
2299 cout << "void PamLevel2::GetWhichTrees(TFile* f) --- WARNING!! --- ...potrebbe non funzionare "<<endl;
2300 // -----------
2301 // reset flags
2302 // -----------
2303 CAL1 = false;
2304 CAL2 = false;
2305 TRK2 = false;
2306 TRK1 = false;
2307 TRKh = false;
2308 TRG = false;
2309 TOF = false;
2310 S4 = false;
2311 ND = false;
2312 AC = false;
2313 ORB = false;
2314
2315 RUN = false;
2316
2317 cout << "Checking file: "<<f->GetName()<<endl;
2318 if( !f || f->IsZombie() ){
2319 cout << "File: "<< f->GetName() <<" Non valid root file"<< endl;
2320 return;
2321 }
2322
2323 TList *lk = f->GetListOfKeys();
2324 TIter next(lk);
2325 TKey *key =0;
2326
2327 Int_t nev = 0;
2328
2329 while( (key = (TKey*)next()) ){
2330
2331 if( !strcmp(key->GetName(),"Run" ) )RUN = true;
2332
2333 //=========================================================
2334 if( !strcmp(key->GetName(),"Trigger" ) ){
2335 TRG = true;
2336 Int_t nevt = ((TTree*)f->Get("Trigger"))->GetEntries();
2337 if( nev && nevt!=nev){
2338 cout << "File: "<< f->GetName() <<" Trigger tree has "<<nevt<<" events instead of "<<nev<< endl;
2339 TRG = false;
2340 }else nev=nevt;
2341 }
2342 //=========================================================
2343 if( !strcmp(key->GetName(),"ToF" ) ){
2344 TOF = true;
2345 Int_t nevt = ((TTree*)f->Get("ToF"))->GetEntries();
2346 if( nev && nevt!=nev){
2347 cout << "File: "<< f->GetName() <<" ToF tree has "<<nevt<<" events instead of "<<nev<< endl;
2348 TOF = false;
2349 }else nev=nevt;
2350 }
2351 //=========================================================
2352 if( !strcmp(key->GetName(),"S4" ) ){
2353 S4 = true;
2354 Int_t nevt = ((TTree*)f->Get("S4"))->GetEntries();
2355 if( nev && nevt!=nev){
2356 cout << "File: "<< f->GetName() <<" S4 tree has "<<nevt<<" events instead of "<<nev<< endl;
2357 S4 = false;
2358 }else nev=nevt;
2359 }
2360 //=========================================================
2361
2362 if( !strcmp(key->GetName(),"NeutronD" ) ){
2363 ND = true;
2364 Int_t nevt = ((TTree*)f->Get("NeutronD"))->GetEntries();
2365 if( nev && nevt!=nev){
2366 cout << "File: "<< f->GetName() <<"NeutronD tree has "<<nevt<<" events instead of "<<nev<< endl;
2367 ND =false;
2368 }else nev=nevt;
2369 }
2370 //=========================================================
2371 if( !strcmp(key->GetName(),"Anticounter") ){
2372 AC = true;
2373 Int_t nevt = ((TTree*)f->Get("Anticounter"))->GetEntries();
2374 if( nev && nevt!=nev){
2375 cout << "File: "<< f->GetName() <<" Anticounter tree has "<<nevt<<" events instead of "<<nev<< endl;
2376 AC =false;
2377 }else nev=nevt;
2378 }
2379 //=========================================================
2380 if( !strcmp(key->GetName(),"OrbitalInfo") ){
2381 ORB = true;
2382 Int_t nevt = ((TTree*)f->Get("OrbitalInfo"))->GetEntries();
2383 if( nev && nevt!=nev){
2384 cout << "File: "<< f->GetName() <<" OrbitalInfo tree has "<<nevt<<" events instead of "<<nev<< endl;
2385 ORB = false;
2386 }else nev=nevt;
2387 }
2388 //=========================================================
2389 if( !strcmp(key->GetName(),"Tracker" ) ){
2390 TTree *T = (TTree*)f->Get("Tracker");
2391 for(Int_t i=0; i<T->GetListOfBranches()->GetEntries(); i++){
2392 TString name = T->GetListOfBranches()->At(i)->GetName();
2393 if( !name.CompareTo("TrkLevel1") )TRK1=true;
2394 if( !name.CompareTo("TrkLevel2") )TRK2=true;
2395 if( !name.CompareTo("TrkHough") )TRKh=true;
2396 };
2397 Int_t nevt = T->GetEntries();
2398 if( nev && nevt!=nev){
2399 cout << "File: "<< f->GetName() <<" Tracker tree has "<<nevt<<" events instead of "<<nev<< endl;
2400 TRK1 = false;
2401 TRK2 = false;
2402 TRKh = false;
2403 }else nev=nevt;
2404 T->Delete();
2405 };
2406 //=========================================================
2407 if( !strcmp(key->GetName(),"Calorimeter" ) ){
2408 TTree *T = (TTree*)f->Get("Calorimeter");
2409 for(Int_t i=0; i<T->GetListOfBranches()->GetEntries(); i++){
2410 TString name = T->GetListOfBranches()->At(i)->GetName();
2411 if( !name.CompareTo("CaloLevel1") )CAL1=true;
2412 if( !name.CompareTo("CaloLevel2") )CAL2=true;
2413 };
2414 Int_t nevt = T->GetEntries();
2415 if( nev && nevt!=nev){
2416 cout << "File: "<< f->GetName() <<" Calorimeter tree has "<<nevt<<" events instead of "<<nev<< endl;
2417 CAL1 = false;
2418 CAL2 = false;
2419 }else nev=nevt;
2420 T->Delete();
2421 };
2422
2423 };
2424
2425 delete lk;
2426
2427 // cout<< "Get detector list from input file --> ";
2428 // if(TRK1)cout<<"TRK1 ";
2429 // if(TRK2)cout<<"TRK2 ";
2430 // if(TRKh)cout<<"TRKH ";
2431 // if(CAL1)cout<<"CAL1 ";
2432 // if(CAL2)cout<<"CAL2 ";
2433 // if(TOF)cout<<"TOF ";
2434 // if(TRG)cout<<"TRG ";
2435 // if(AC)cout<<"AC ";
2436 // if(ND)cout<<"ND ";
2437 // if(S4)cout<<"S4 ";
2438 // if(ORB)cout<<"ORB ";
2439 // cout << endl;
2440
2441 return ;
2442
2443 };
2444
2445
2446 //--------------------------------------
2447 //
2448 //
2449 //--------------------------------------
2450 /**
2451 * Check if a file contains selected Pamela Level2 trees.
2452 * @param name File name
2453 * @return true if the file is ok.
2454 */
2455 Bool_t PamLevel2::CheckLevel2File(TString name){
2456
2457 Bool_t CAL1__ok = false;
2458 Bool_t CAL2__ok = false;
2459 Bool_t TRK2__ok = false;
2460 Bool_t TRK1__ok = false;
2461 Bool_t TRKh__ok = false;
2462 Bool_t TRG__ok = false;
2463 Bool_t TOF__ok = false;
2464 Bool_t S4__ok = false;
2465 Bool_t ND__ok = false;
2466 Bool_t AC__ok = false;
2467 Bool_t ORB__ok = false;
2468
2469 Bool_t RUN__ok = false;
2470
2471 Bool_t SELLI__ok = false;
2472
2473 cout << "Checking file: "<<name<<endl;
2474 TFile *f = new TFile(name.Data());
2475 if( !f || f->IsZombie() ){
2476 cout << "File: "<< f->GetName() <<" discarded ---- Non valid root file"<< endl; return false;
2477 }
2478 // cout << "Get list of keys: "<<f<<endl;
2479 TList *lk = f->GetListOfKeys();
2480 // lk->Print();
2481 TIter next(lk);
2482 TKey *key =0;
2483
2484 Int_t nev = 0;
2485
2486 while( (key = (TKey*)next()) ){
2487
2488 // cout << key->GetName() << endl;
2489 // cout << key->GetName() << ""<<key->GetClassName()<<endl;
2490 // cout << " Get tree: " << f->Get(key->GetName())<<endl;
2491 // nev_previous = nev;
2492 // cout << " n.entries "<< nev <<endl;
2493 // if( key->GetClassName()=="TTree" && nev_previous && nev != nev_previous ){
2494 // nev = ((TTree*)f->Get(key->GetName()))->GetEntries();
2495 // cout << "File: "<< f->GetName() <<" discarded ---- "<< key->GetName() << " tree: n.entries does not match "<<nev<<" "<<nev_previous<< endl;
2496 // return false;
2497 // };
2498
2499
2500 if( !strcmp(key->GetName(),"Run" ) )RUN__ok = true;
2501
2502 //=========================================================
2503 if( !strcmp(key->GetName(),"SelectionList" ) ){
2504 SELLI__ok = true;
2505 if(SELLI==1){
2506 Int_t nevt = ((TTree*)f->Get("SelectionList"))->GetEntries();
2507 if( nev && nevt!=nev){
2508 cout << "File: "<< f->GetName() <<" discarded ---- SelectionList tree has "<<nevt<<" events instead of "<<nev<< endl;
2509 return false;
2510 }
2511 nev=nevt;
2512 }
2513 }
2514
2515 //=========================================================
2516 if( !strcmp(key->GetName(),"Trigger" ) ){
2517 TRG__ok = true;
2518 if(TRG){
2519 Int_t nevt = ((TTree*)f->Get("Trigger"))->GetEntries();
2520 if( nev && nevt!=nev){
2521 cout << "File: "<< f->GetName() <<" discarded ---- Trigger tree has "<<nevt<<" events instead of "<<nev<< endl;
2522 return false;
2523 }
2524 nev=nevt;
2525 }
2526 }
2527 //=========================================================
2528 if( !strcmp(key->GetName(),"ToF" ) ){
2529 TOF__ok = true;
2530 if(TOF){
2531 Int_t nevt = ((TTree*)f->Get("ToF"))->GetEntries();
2532 if( nev && nevt!=nev){
2533 cout << "File: "<< f->GetName() <<" discarded ---- ToF tree has "<<nevt<<" events instead of "<<nev<< endl;
2534 return false;
2535 }
2536 nev=nevt;
2537 }
2538 }
2539 //=========================================================
2540 if( !strcmp(key->GetName(),"S4" ) ){
2541 S4__ok = true;
2542 if(S4){
2543 Int_t nevt = ((TTree*)f->Get("S4"))->GetEntries();
2544 if( nev && nevt!=nev){
2545 cout << "File: "<< f->GetName() <<" discarded ---- S4 tree has "<<nevt<<" events instead of "<<nev<< endl;
2546 return false;
2547 }
2548 nev=nevt;
2549 }
2550 }
2551 //=========================================================
2552
2553 if( !strcmp(key->GetName(),"NeutronD" ) ){
2554 ND__ok = true;
2555 if(ND){
2556 Int_t nevt = ((TTree*)f->Get("NeutronD"))->GetEntries();
2557 if( nev && nevt!=nev){
2558 cout << "File: "<< f->GetName() <<" discarded ---- NeutronD tree has "<<nevt<<" events instead of "<<nev<< endl;
2559 return false;
2560 }
2561 nev=nevt;
2562 }
2563 }
2564 //=========================================================
2565 if( !strcmp(key->GetName(),"Anticounter") ){
2566 AC__ok = true;
2567 if(AC){
2568 Int_t nevt = ((TTree*)f->Get("Anticounter"))->GetEntries();
2569 if( nev && nevt!=nev){
2570 cout << "File: "<< f->GetName() <<" discarded ---- Anticounter tree has "<<nevt<<" events instead of "<<nev<< endl;
2571 return false;
2572 }
2573 nev=nevt;
2574 }
2575 }
2576 //=========================================================
2577 if( !strcmp(key->GetName(),"OrbitalInfo") ){
2578 ORB__ok = true;
2579 if(ORB){
2580 Int_t nevt = ((TTree*)f->Get("OrbitalInfo"))->GetEntries();
2581 if( nev && nevt!=nev){
2582 cout << "File: "<< f->GetName() <<" discarded ---- OrbitalInfo tree has "<<nevt<<" events instead of "<<nev<< endl;
2583 return false;
2584 }
2585 nev=nevt;
2586 }
2587 }
2588 //=========================================================
2589 if( !strcmp(key->GetName(),"Tracker" ) ){
2590 TTree *T = (TTree*)f->Get("Tracker");
2591 if(TRK1||TRK2||TRKh){
2592 Int_t nevt = T->GetEntries();
2593 if( nev && nevt!=nev){
2594 cout << "File: "<< f->GetName() <<" discarded ---- Tracker tree has "<<nevt<<" events instead of "<<nev<< endl;
2595 return false;
2596 }
2597 nev=nevt;
2598 }
2599 for(Int_t i=0; i<T->GetListOfBranches()->GetEntries(); i++){
2600 TString name = T->GetListOfBranches()->At(i)->GetName();
2601 if( !name.CompareTo("TrkLevel1") )TRK1__ok=true;
2602 if( !name.CompareTo("TrkLevel2") )TRK2__ok=true;
2603 if( !name.CompareTo("TrkHough") )TRKh__ok=true;
2604 };
2605 T->Delete();
2606 };
2607 //=========================================================
2608 if( !strcmp(key->GetName(),"Calorimeter" ) ){
2609 TTree *T = (TTree*)f->Get("Calorimeter");
2610 if(CAL1||CAL2){
2611 Int_t nevt = T->GetEntries();
2612 if( nev && nevt!=nev){
2613 cout << "File: "<< f->GetName() <<" discarded ---- Calorimeter tree has "<<nevt<<" events instead of "<<nev<< endl;
2614 return false;
2615 }
2616 nev=nevt;
2617 }
2618 for(Int_t i=0; i<T->GetListOfBranches()->GetEntries(); i++){
2619 TString name = T->GetListOfBranches()->At(i)->GetName();
2620 if( !name.CompareTo("CaloLevel1") )CAL1__ok=true;
2621 if( !name.CompareTo("CaloLevel2") )CAL2__ok=true;
2622 };
2623 T->Delete();
2624 };
2625
2626 };
2627
2628 if( SELLI==-1 )SELLI = (Int_t)SELLI__ok;
2629 if( SELLI==0 && SELLI__ok ){
2630 cout << "File: "<< f->GetName() <<" discarded ---- found SelectionList (it is not a full-event file)" << endl;
2631 return false;
2632 }
2633 if( SELLI==1 && !SELLI__ok ){
2634 cout << "File: "<< f->GetName() <<" discarded ---- SelectionList missing" << endl;
2635 return false;
2636 }
2637
2638 // cout << "SELLI "<<SELLI<<endl;
2639
2640 // cout<< "CheckLevel2File(TString): detector list --> ";
2641 // if(TRK1__ok)cout<<"TRK1 ";
2642 // if(TRK2__ok)cout<<"TRK2 ";
2643 // if(TRKh__ok)cout<<"TRKH ";
2644 // if(CAL1__ok)cout<<"CAL1 ";
2645 // if(CAL2__ok)cout<<"CAL2 ";
2646 // if(TOF__ok)cout<<"TOF ";
2647 // if(TRG__ok)cout<<"TRG ";
2648 // if(AC__ok)cout<<"AC ";
2649 // if(ND__ok)cout<<"ND ";
2650 // if(S4__ok)cout<<"S4 ";
2651 // if(ORB__ok)cout<<"ORB ";
2652 // cout << endl;
2653
2654
2655 if(TRK2 && TRK1__ok)TRK1=1;
2656 // ----------------------------------------------------------------------------
2657 // NOTA
2658 // se c'e` il level1, lo devo necessarimente leggere.
2659 // infatti (non ho capito perche`) i cluster vengono letti e allocati in memoria
2660 // comunque, ma non vengono disallocati da PamLevel2::Clear()
2661 // ----------------------------------------------------------------------------
2662
2663
2664 if(!RUN__ok) {
2665 cout << "File: "<< f->GetName() <<" *WARNING* ---- Missing RunInfo tree"<< endl;
2666 // return false;
2667 };
2668
2669 if(CAL1 && !CAL1__ok){
2670 cout << "File: "<< f->GetName() <<" discarded ---- Missing CaloLevel1 branch"<< endl;
2671 return false;
2672 };
2673 if(CAL2 && !CAL2__ok){
2674 cout << "File: "<< f->GetName() <<" discarded ---- Missing CaloLevel2 branch"<< endl;
2675 return false;
2676 };
2677 if(TRK2 && !TRK2__ok){
2678 cout << "File: "<< f->GetName() <<" discarded ---- Missing TrkLevel2 branch"<< endl;
2679 return false;
2680 };
2681 if(TRK1 && !TRK1__ok){
2682 cout << "File: "<< f->GetName() <<" discarded ---- Missing TrkLevel1 branch"<< endl;
2683 return false;
2684 };
2685 if(TRKh && !TRKh__ok){
2686 cout << "File: "<< f->GetName() <<" discarded ---- Missing TrkHough branch"<< endl;
2687 return false;
2688 };
2689 if(ORB && !ORB__ok){
2690 cout << "File: "<< f->GetName() <<" discarded ---- Missing ORB tree"<< endl;
2691 return false;
2692 };
2693 if(AC && !AC__ok){
2694 cout << "File: "<< f->GetName() <<" discarded ---- Missing AC tree"<< endl;
2695 return false;
2696 };
2697 if(S4 && !S4__ok){
2698 cout << "File: "<< f->GetName() <<" discarded ---- Missing S4 tree"<< endl;
2699 return false;
2700 };
2701 if(TOF && !TOF__ok){
2702 cout << "File: "<< f->GetName() <<" discarded ---- Missing ToF tree"<< endl;
2703 return false;
2704 };
2705
2706 if(ND && !ND__ok){
2707 cout << "File: "<< f->GetName() <<" discarded ---- Missing ND tree"<< endl;
2708 return false;
2709 };
2710 if(TRG && !TRG__ok){
2711 cout << "File: "<< f->GetName() <<" discarded ---- Missing Trigger tree"<< endl;
2712 return false;
2713 };
2714
2715
2716 // lk->Delete();
2717 // delete lk;
2718 f->Close();
2719
2720 // cout<< "CheckLevel2File(TString): detector list --> ";
2721 // if(TRK1)cout<<"TRK1 ";
2722 // if(TRK2)cout<<"TRK2 ";
2723 // if(TRKh)cout<<"TRKH ";
2724 // if(CAL1)cout<<"CAL1 ";
2725 // if(CAL2)cout<<"CAL2 ";
2726 // if(TOF)cout<<"TOF ";
2727 // if(TRG)cout<<"TRG ";
2728 // if(AC)cout<<"AC ";
2729 // if(ND)cout<<"ND ";
2730 // if(S4)cout<<"S4 ";
2731 // if(ORB)cout<<"ORB ";
2732 // cout << endl;
2733
2734 return true;
2735
2736 };
2737
2738
2739 /**
2740 * Create clone-trees
2741 */
2742 void PamLevel2::CreateCloneTrees0( TChain *fChain, TFile *ofile ){
2743
2744 cout << "+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+" <<endl;
2745 cout << "Create clones of PAMELA trees "<<endl;
2746
2747 Int_t i=0;
2748 pam_tree_clone[i] = fChain->GetTree()->CloneTree(0);
2749 TString name = pam_tree_clone[i]->GetName();
2750 name.Append("_clone");
2751 // pam_tree_clone[i]->SetName(name.Data());
2752 cout << pam_tree_clone[i]->GetName() <<endl;
2753 i++;
2754
2755 TList *li = fChain->GetListOfFriends();
2756 TIter next(li);
2757 TFriendElement* T_friend=0;
2758 ofile->cd();
2759 while( (T_friend = (TFriendElement*)next()) ){
2760 // cout<<T_friend->IsA()->GetName()<<" "<<T_friend->GetName()<<hex << T_friend->GetTree() << dec<<endl;
2761 // cout<<T_friend->GetTree()->GetName()<< endl;
2762 pam_tree_clone[i] = T_friend->GetTree()->CloneTree(0);
2763 pam_tree_clone[i]->SetAutoSave(1000000);
2764 name = pam_tree_clone[i]->GetName();
2765 name.Append("_clone");
2766 // pam_tree_clone[i]->SetName(name.Data());
2767 cout << pam_tree_clone[i]->GetName() << endl;
2768 i++;
2769 }
2770
2771 delete li;
2772
2773 cout << "+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+" <<endl;
2774
2775 }
2776
2777 /**
2778 * Create clone-trees
2779 */
2780 void PamLevel2::CreateCloneTrees(TFile *ofile){
2781
2782
2783 // if the pointer is null, create a default file
2784
2785 if(!ofile){
2786 cout << "void PamLevel2::CreateCloneTrees(TFile*) -- WARNING -- Creating file: clone-tree.root "<<endl;
2787 ofile = new TFile("clone-tree.root","recreate");
2788 }
2789
2790 ofile->cd();
2791
2792 cout << "+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+" <<endl;
2793 cout << "Create new PAMELA trees "<<endl;
2794
2795
2796 run_tree_clone = new TTree("Run","PAMELA Level2 data from the GL_RUN table ");
2797 run_tree_clone->Branch("RunInfo","GL_RUN",GetPointerTo("RunInfo"));
2798 cout << "Run : branch RunInfo"<<endl;
2799 run_tree_clone->Branch("SoftInfo","SoftInfo",GetPointerTo("SoftInfo"));
2800 cout << "Run : branch SoftInfo"<<endl;
2801 // ------------------
2802 // replicate run tree
2803 // ------------------
2804 cout << "----------------------------------------------------"<<endl;
2805 cout << "irun\t | RUN\t NEVENTS\t absolute time"<<endl;
2806 for (Int_t i=0; i<run_tree->GetEntries(); i++){
2807 run_tree->GetEntry(i);
2808 cout << i<< "\t | "<<GetRunInfo()->ID<<"\t "<<GetRunInfo()->NEVENTS<< "\t "<<GetRunInfo()->RUNHEADER_TIME<<" <---> "<<GetRunInfo()->RUNTRAILER_TIME<<endl;
2809 run_tree_clone->Fill();
2810 }
2811 cout << "----------------------------------------------------"<<endl;
2812
2813
2814 sel_tree_clone = new TTree("SelectionList","List of selected events ");
2815 sel_tree_clone->Branch("RunEntry",&irun,"runentry/L");
2816 sel_tree_clone->Branch("EventEntry",&irunentry,"eventry/L");
2817
2818
2819 Int_t i=0;
2820 if(TRK1||TRK2||TRKh){
2821 pam_tree_clone[i] = new TTree("Tracker","PAMELA tracker level2 data ");
2822 if(TRK1) {
2823 pam_tree_clone[i]->Branch("TrkLevel1","TrkLevel1", GetPointerTo("TrkLevel1"));
2824 pam_tree_clone[i]->BranchRef();
2825 cout << "Tracker : branch TrkLevel1"<<endl;
2826 // cout << "CreateCloneTrees " << GetTrkLevel1()<<endl;
2827 };
2828 if(TRK2) {
2829 pam_tree_clone[i]->Branch("TrkLevel2", "TrkLevel2",GetPointerTo("TrkLevel2"));
2830 cout << "Tracker : branch TrkLevel2"<<endl;
2831 };
2832 if(TRKh) {
2833 pam_tree_clone[i]->Branch("TrkHough","TrkHough", GetPointerTo("TrkHough"));
2834 cout << "Tracker : branch TrkHough"<<endl;
2835 };
2836 i++;
2837 }
2838
2839 // Calorimeter
2840 if(CAL1||CAL2){
2841 pam_tree_clone[i] = new TTree("Calorimeter","PAMELA calorimeter level2 data ");
2842 if(CAL1) {
2843 pam_tree_clone[i]->Branch("CaloLevel1", "CaloLevel1", GetPointerTo("CaloLevel1"));
2844 cout << "Calorimeter : branch CaloLevel1"<<endl;
2845 };
2846 if(CAL2) {
2847 pam_tree_clone[i]->Branch("CaloLevel2","CaloLevel2", GetPointerTo("CaloLevel2"));
2848 cout << "Calorimeter : branch CaloLevel2"<<endl;
2849 };
2850 i++;
2851 }
2852
2853 // ToF
2854 if(TOF) {
2855 pam_tree_clone[i] = new TTree("ToF","PAMELA ToF level2 data ");
2856 pam_tree_clone[i]->Branch("ToFLevel2","ToFLevel2", GetPointerTo("ToFLevel2"));
2857 cout << "ToF : branch ToFLevel2"<<endl;
2858 i++;
2859 };
2860 // Trigger
2861 if(TRG) {
2862 pam_tree_clone[i] = new TTree("Trigger","PAMELA trigger level2 data ");
2863 pam_tree_clone[i]->Branch("TrigLevel2","TrigLevel2", GetPointerTo("TrigLevel2"));
2864 cout << "Trigger : branch TrigLevel2"<<endl;
2865 i++;
2866 };
2867 // S4
2868 if(S4) {
2869 pam_tree_clone[i] = new TTree("S4","PAMELA S4 level2 data ");
2870 pam_tree_clone[i]->Branch("S4Level2","S4Level2", GetPointerTo("S4Level2"));
2871 cout << "S4 : branch S4Level2"<<endl;
2872 i++;
2873 };
2874 // Neutron Detector
2875 if(ND) {
2876 pam_tree_clone[i] = new TTree("NeutronD","PAMELA neutron detector level2 data ");
2877 pam_tree_clone[i]->Branch("NDLevel2","NDLevel2", GetPointerTo("NDLevel2"));
2878 cout << "NeutronD : branch NDLevel2"<<endl;
2879 i++;
2880 };
2881 // Anticounters
2882 if(AC) {
2883 pam_tree_clone[i] = new TTree("Anticounter","PAMELA anticounter detector level2 data ");
2884 pam_tree_clone[i]->Branch("AcLevel2","AcLevel2", GetPointerTo("AcLevel2"));
2885 cout << "Anticounter : branch AcLevel2"<<endl;
2886 i++;
2887 };
2888 // OrbitalInfo
2889 if(ORB) {
2890 pam_tree_clone[i] = new TTree("OrbitalInfo","PAMELA oribital info ");
2891 pam_tree_clone[i]->Branch("OrbitalInfo","OrbitalInfo", GetPointerTo("OrbitalInfo"));
2892 cout << "OrbitalInfo : branch OrbitalInfo"<<endl;
2893 i++;
2894 };
2895 cout << "+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+" <<endl;
2896
2897 }
2898
2899
2900 /**
2901 * Fill tree (created with CreateCloneTrees)
2902 *
2903 */
2904 //void PamLevel2::FillNewPamTree(TTree *T){
2905 void PamLevel2::FillCloneTrees(){
2906
2907 // cout << "PamLevel2::FillCloneTrees()" << irunentry << endl;
2908
2909 for(Int_t i=0; i<NCLONES; i++){
2910 if(pam_tree_clone[i])pam_tree_clone[i]->Fill();
2911 }
2912 if(sel_tree_clone)sel_tree_clone->Fill();
2913
2914 }
2915
2916
2917 TTree* PamLevel2::GetCloneTree(TString name){
2918
2919 for(Int_t i=0; i<NCLONES; i++){
2920 if(pam_tree_clone[i]){
2921 TString na = pam_tree_clone[i]->GetName();
2922 if(!name.CompareTo(na))return pam_tree_clone[i];
2923 };
2924 }
2925 if(run_tree_clone){
2926 TString na = run_tree_clone->GetName();
2927 if(!name.CompareTo(na))return run_tree_clone;
2928 }
2929 if(sel_tree_clone){
2930 TString na = sel_tree_clone->GetName();
2931 if(!name.CompareTo(na))return sel_tree_clone;
2932 }
2933 return NULL;
2934
2935 }
2936 void PamLevel2::WriteCloneTrees(){
2937 cout << "+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+" <<endl;
2938 cout << "Write clones of PAMELA trees "<<endl;
2939 cout << run_tree_clone->GetName()<<endl;
2940 run_tree_clone->Write();
2941 cout << sel_tree_clone->GetName()<<endl;
2942 sel_tree_clone->Write();
2943 for(Int_t i=0; i<NCLONES; i++){
2944 if(pam_tree_clone[i]){
2945 cout << pam_tree_clone[i]->GetName()<<endl;
2946 pam_tree_clone[i]->Write();
2947 };
2948 }
2949 cout << "+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+" <<endl;
2950
2951 }
2952
2953 /**
2954 * Method to get level2-trees entry, the corresponding run entry and (if required) the level0 entry.
2955 */
2956 //Int_t PamLevel2::GetEntry(Int_t iee){
2957 Int_t PamLevel2::GetEntry(Long64_t iee){
2958
2959 // cout << "-------------------------------------"<<endl;
2960 // cout << "Int_t PamLevel2::GetEntry("<<iee<<")"<<endl;
2961
2962 if(!pam_tree){
2963 cout << " Int_t PamLevel2::GetEntry(Int_t) -- ERROR -- level2 trees not loaded"<<endl;
2964 return 0;
2965 }
2966
2967
2968 //
2969 // This is a sort of bug: if you don't have the run tree you don't want to exit here you want to have loaded the event anyway...
2970 //
2971 // if(!run_tree ){
2972 // cout << " Int_t PamLevel2::GetEntry(Int_t) -- ERROR -- run tree not loeaded"<<endl;
2973 // return 0;
2974 // }
2975
2976 Long64_t ii=0;
2977 //-------------------------------
2978 ii = iee;
2979 if( !pam_tree->GetEntry(ii) ){
2980 cout << " Int_t PamLevel2::GetEntry(Int_t) -- ERROR -- error reading pam tree"<<endl;
2981 return 0;
2982 }
2983 //
2984 // ... that's way I put it here. Notice that nothing change in the code (is backward compatible) since in any case you return with 0.
2985 // in theory one would like to return 1 if run is not loaded but now I don't have the will to add that 2 lines of code and it is not
2986 // a problem if you don't check the return code of getentry.
2987 //
2988 if(!run_tree ){
2989 if ( TRK0 || CAL0 || TOF0 || RUN ) { //forse cosi` va bene per tornare 1?
2990 cout << " Int_t PamLevel2::GetEntry(Int_t) -- ERROR -- run tree not loaded"<<endl;
2991 return 0;
2992 } else {
2993 return 1; //cosi` se non c'e` run esce qua...
2994 }
2995 }
2996
2997 //-------------------------------
2998 ii = iee;
2999 // Bool_t UPDATED = UpdateRunInfo(run_tree,ii);
3000 // Bool_t UPDATED = UpdateRunInfo(ii);
3001 UpdateRunInfo(ii);
3002 if(SELLI==0)irunentry = iee-runfirstentry;
3003 // if(UPDATED && run_tree_clone)run_tree_clone->Fill();
3004
3005 // cout << "PamLevel2::GetEntry("<<iee<<") "<<irun<<" "<<runfirstentry<<" "<<irunentry<<endl;
3006
3007 // cout << "irunentry "<<irunentry << endl;
3008 // cout << "runfirstentry "<<runfirstentry << endl;
3009 // cout << "nevents "<<GetRunInfo()->NEVENTS<< endl;
3010
3011 // if( TRK0 || CAL0 || TOF0 ){
3012 // if( !GetYodaEntry( ) ){
3013 // cout << " Int_t PamLevel2::GetEntry(Int_t) -- ERROR -- error reading level0 tree"<<endl;
3014 // return 0;
3015 // }
3016 // }
3017
3018
3019 return 1;
3020
3021 }
3022
3023 TrkLevel0 *PamLevel2::GetTrkLevel0(){
3024 if( !TRK0 )return NULL;
3025 if( !GetYodaEntry( ) ){
3026 cout << " Int_t PamLevel2::GetTrkLevel0() -- ERROR -- error reading level0 tree"<<endl;
3027 return 0;
3028 }
3029 return trk0_obj;
3030 };
3031 CaloLevel0 *PamLevel2::GetCaloLevel0(){
3032 if( !CAL0 )return NULL;
3033 if( !GetYodaEntry( ) ){
3034 cout << " Int_t PamLevel2::GetCaloLevel0() -- ERROR -- error reading level0 tree"<<endl;
3035 return 0;
3036 }
3037 return calo0_obj;
3038 };
3039
3040
3041 /**
3042 * Method to retrieve the level0 tree (YODA tree) that contains the current event.
3043 * Given the run ID (...), if needed it query the DB and load the proper file.
3044 * @return Pointer to the tree
3045 */
3046
3047
3048 TTree* PamLevel2::GetYodaTree( ){
3049
3050 // cout << "TTree* PamLevel2::GetYodaTree( )"<<endl;
3051 //===================================
3052 // check if iroot has changed
3053 //===================================
3054 if( irun<0 ){
3055 cout << "PamLevel2::GetYodaTree() -- ERROR -- irun = "<<irun<<endl;
3056 // cout << "In order to use this method you have to first load the RunInfo tree "<<endl;
3057 return NULL;
3058 }
3059 Int_t irootnew = run_obj->ID_ROOT_L0;
3060 // cout << "iroot "<<iroot<<endl;
3061 // cout << "irootnew "<<irootnew<<endl;
3062
3063 //===================================
3064 // load the level0 file
3065 // (if not already loaded)
3066 //===================================
3067 if( iroot != irootnew || !l0_tree){
3068 iroot = irootnew;
3069 //===================================
3070 // open the DB connection
3071 // (if not already opened)
3072 //===================================
3073 if(!dbc || (dbc && !dbc->IsConnected())){
3074 cout<<"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;
3075 cout<<"Connecting to DB"<<endl;
3076 cout<<"HOST "<<host<<endl;
3077 cout<<"USER "<<user<<endl;
3078 cout<<"PSW "<<psw<<endl;
3079 dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
3080 if( !dbc )return NULL;
3081 if( !dbc->IsConnected() )return NULL;
3082 // cout<<"...done"<<endl;
3083 cout<<"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;
3084 }else{
3085 // cout<<"DB already connected"<<endl;
3086 }
3087 GL_ROOT glroot = GL_ROOT();
3088 if( glroot.Query_GL_ROOT(iroot,dbc) ){
3089 cout << "TTree* PamLevel2::GetYodaTree( ) -- ERROR -- level0 file iroot = "<<iroot<< " does not exists"<<endl;
3090 return NULL;
3091 };
3092 TString filename = glroot.PATH + glroot.NAME;
3093 if(l0_file){
3094 l0_file->Close();
3095 l0_file->Delete();
3096 }
3097 cout << "Opening LEVEL0 file: "<< filename << endl;
3098 FileStat_t t;
3099 if( gSystem->GetPathInfo(filename.Data(),t) ){
3100 cout << " PamLevel2::GetYodaTree() -- ERROR opening file "<<endl;
3101 return NULL;
3102 }
3103 l0_file = new TFile(filename);
3104 if( !l0_file )return NULL;
3105 l0_tree = (TTree*)l0_file->Get("Physics");
3106 if(!h0_obj)h0_obj = new EventHeader();
3107 l0_tree->SetBranchAddress("Header" ,&h0_obj);
3108 prevshift = 0;
3109 //---------------------------------------------------
3110 // TRACKER:
3111 if(TRK0){
3112 if(!trk0_obj){
3113 trk0_obj = new TrkLevel0();
3114 trk0_obj->Set();
3115 };
3116 l0_tree->SetBranchAddress("Tracker" ,trk0_obj->GetPointerToTrackerEvent());
3117 }
3118 //--------------------------------------------------
3119 // CALORIMETER:
3120 if(CAL0){
3121 if(!calo0_obj){
3122 calo0_obj = new CaloLevel0();
3123 calo0_obj->Set();
3124 };
3125 l0_tree->SetBranchAddress("Calorimeter" ,calo0_obj->GetPointerToCalorimeterEvent());
3126 // cout << "PamLevel2::GetYodaTree() --- level0 calorimeter not implemented "<<endl;
3127 }
3128 //---------------------------------------------------
3129 // TOF:
3130 if(TOF0){
3131 cout << "PamLevel2::GetYodaTree() --- level0 TOF not implemented "<<endl;
3132 }
3133
3134 };
3135
3136 if(!dbc || (dbc && !dbc->IsConnected())){
3137 cout << " TTree* PamLevel2::GetYodaTree( ) -- no DB connected... hai fatto qualche cazzata "<<endl;
3138 }
3139
3140 if ( TRK0 ){
3141 TrkParams::Load(6);
3142 if( !TrkParams::IsLoaded(6) ){
3143 cout << " TTree* PamLevel2::GetYodaTree( ) -- WARNING -- VK-mask not loaded"<<endl;
3144 };
3145 TrkParams::SetCalib(run_obj,dbc);
3146 TrkParams::LoadCalib( );
3147 if( !TrkParams::CalibIsLoaded() ){
3148 cout << " TTree* PamLevel2::GetYodaTree( ) -- WARNING -- Calibration not loaded"<<endl;
3149 };
3150 }
3151
3152 // cout << l0_tree << endl;
3153
3154 return l0_tree;
3155
3156 }
3157
3158 /**
3159 * Method to retrieve the level0 tree (YODA tree) that contains the current event.
3160 */
3161 Int_t PamLevel2::GetYodaEntry(){
3162
3163 // cout << "Int_t PamLevel2::GetYodaEntry()"<<endl;
3164 if(!GetYodaTree())return 0;
3165
3166 // patch
3167 if( irunentry < 0){
3168 // cout << "Int_t PamLevel2::GetYodaEntry() -- ATTENZIONE -- irunentry negativo?!?! "<<(Int_t)irunentry<<endl;
3169 irunentry=0LL;
3170 }
3171 // ---------------------------------
3172 // if file is NOT a preselected file
3173 // ---------------------------------
3174 Long64_t quellagiusta = irunentry + (Long64_t)(run_obj->EV_FROM);
3175 // cout << " irunentry "<<irunentry<<endl;
3176 // cout << " EV_FROM "<<run_obj->EV_FROM<<endl;
3177 // cout << " quellagiusta = irunentry + EV_FROM "<< quellagiusta << endl;
3178
3179 // cout << " irun "<< irun << " irunentry "<< irunentry<<" run_obj->EV_FROM "<<run_obj->EV_FROM <<" quella giusta "<<quellagiusta << endl;
3180 // cout << " iroot "<<iroot<<" run_obj->ID_ROOT_L0 "<<run_obj->ID_ROOT_L0<<endl;
3181 // cout << " time "<< GetOrbitalInfo()->absTime << endl;
3182 // cout << " trk_calib_used "<<run_obj->TRK_CALIB_USED<< endl;
3183
3184 if( !GetOrbitalInfo() ){
3185 cout << "Int_t PamLevel2::GetYodaEntry() -- ERROR -- missing OrbitalInfo "<<endl;
3186 return 0;
3187 }
3188 if( GetOrbitalInfo()->OBT==0 && GetOrbitalInfo()->pkt_num==0 ){
3189 cout << "Int_t PamLevel2::GetYodaEntry() -- ERROR -- level2 event corrupted ?? "<<endl;
3190 return 0;
3191 }
3192
3193 // ---------------------------------------------------------------------
3194 // ATTENTION!!!
3195 // If data are not pre-processed with cleaner, the level0 tree may contain
3196 // spurious nested physics packets.
3197 // The GL_RUN variables EV_FROM, EV_TO and NEVENTS counts also these entries
3198 // while level2 tree DOES NOT!!
3199 // This means that "quellagiusta" in these cases is not correct.
3200 // In order to retrieve the correct level0 event, I implemented a check
3201 // of the OBT and pkt-number. In case of mismatch, the level0 entry number
3202 // is shift forward until when the packets match.
3203 // ---------------------------------------------------------------------
3204 Int_t answer = 0;
3205 Int_t shift =0;
3206 // printf(" siamo qui %i %i \n",shift,prevshift);
3207 Int_t maxshift = 10;
3208 do{
3209 if(shift>0){
3210 cout << " PKTNUM L2 --- "<< GetOrbitalInfo()->pkt_num << " --- L0 --- "<< GetEventHeader()->GetPscuHeader()->GetCounter()<<endl;
3211 cout << " RUN: ID "<< GetRunInfo()->ID << " ID_ROOT_L0 "<<run_obj->ID_ROOT_L0<<" ID_RUN_FRAG "<<GetRunInfo()->ID_RUN_FRAG << " EV_FROM "<<GetRunInfo()->EV_FROM <<endl;
3212 cout << " L2 <--> L0 mismatch ( irun "<<irun<<" irunentry "<<irunentry<<" shift "<<shift<<" prevshift "<<prevshift<<" )"<<endl;
3213 }
3214 answer = l0_tree->GetEntry(quellagiusta+(Long64_t)shift+(Long64_t)prevshift);
3215 shift++;
3216 if( !GetEventHeader() ){
3217 cout << "Int_t PamLevel2::GetYodaEntry() -- ERROR -- missing EventHeader "<<endl;
3218 return 0;
3219 }
3220 // cout << "PKTNUM "<<shift<<" == L2 --- "<< GetOrbitalInfo()->pkt_num << " --- L0 --- "<< GetEventHeader()->GetPscuHeader()->GetCounter()<<endl;
3221 // cout << " L2 --- "<< GetOrbitalInfo()->OBT << " --- L0 --- "<< GetEventHeader()->GetPscuHeader()->GetOrbitalTime()<<endl;
3222 // if( (quellagiusta+shift) == l0_tree->GetEntries()+1 )cout << ">>> end of level0 tree <<<"<<endl;
3223 // cout << " GetOrbitalInfo()->OBT "<< GetOrbitalInfo()->OBT << endl;
3224 // cout << " GetEventHeader()->GetPscuHeader()->GetOrbitalTime() "<< GetEventHeader()->GetPscuHeader()->GetOrbitalTime() << endl;
3225 // cout << " GetOrbitalInfo()->pkt_num "<< GetOrbitalInfo()->pkt_num << endl;
3226 // cout << " GetEventHeader()->GetPscuHeader()->GetCounter() "<< GetEventHeader()->GetPscuHeader()->GetCounter() << endl;
3227 // printf(" IDRUN %u \n",GetRunInfo()->ID);
3228 //
3229 if ( prevshift != 0 && (quellagiusta+(Long64_t)shift) == GetYodaTree()->GetEntries() ){
3230 prevshift = 0;
3231 shift = -1;
3232 };
3233
3234 }while( ( GetOrbitalInfo()->OBT != (UInt_t)(GetEventHeader()->GetPscuHeader()->GetOrbitalTime()) || GetOrbitalInfo()->pkt_num != (UInt_t)(GetEventHeader()->GetPscuHeader()->GetCounter())) && (quellagiusta+(Long64_t)shift) < GetYodaTree()->GetEntries() && shift < maxshift);
3235
3236 if ( (quellagiusta+(Long64_t)shift+(Long64_t)prevshift) > GetYodaTree()->GetEntries() || shift == maxshift ) {
3237 cout << " Big trouble here, no such event in Level0 data! " <<endl;
3238 return 0;
3239 }
3240 // cout << "LA ENTRY GIUSTA E`: "<<quellagiusta<<" (spero...)"<<endl;
3241 // return GetYodaTree()->GetEntry(quellagiusta);
3242 if ( shift > 1 ) prevshift += (shift-1);
3243
3244 return answer;
3245
3246 }
3247 /**
3248 * \Brief Set DB connection
3249 */
3250 Bool_t PamLevel2::SetDBConnection(){
3251
3252 cout<<"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;
3253 cout<<"Connecting to DB"<<endl;
3254 cout<<"HOST "<<host<<endl;
3255 cout<<"USER "<<user<<endl;
3256 cout<<"PSW "<<psw<<endl;
3257 dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
3258 if( !dbc )return false;
3259 if( !dbc->IsConnected() )return false;
3260 cout<<"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;
3261 return true;
3262
3263 }

  ViewVC Help
Powered by ViewVC 1.1.23